Climate of the Past Climate change between the mid and late Holocene in northern high latitudes – Part 2 : Model-data comparisons

The climate response over northern high latitudes to the mid-Holocene orbital forcing has been investigated in three types of PMIP (Paleoclimate Modelling Intercomparison Project) simulations with different complexity of the modelled climate system. By first undertaking model-data comparison, an objective selection method has been applied to evaluate the capability of the climate models to reproduce the spatial response pattern seen in proxy data. The possible feedback mechanisms behind the climate response have been explored based on the selected model simulations. Subsequent model-model comparisons indicate the importance of including the different physical feedbacks in the climate models. The comparisons between the proxy-based reconstructions and the best fit selected simulations show that over the northern high latitudes, summer temperature change follows closely the insolation change and shows a common feature with strong warming over land and relatively weak warming over ocean at 6 ka compared to 0 ka. Furthermore, the sea-ice-albedo positive feedback enhances this response. The reconstructions of temperature show a stronger response to enhanced insolation in the annual mean temperature than winter and summer temperature. This is verified in the model simulations and the behaviour is attributed to the larger contribution from the large response in autumn. Despite a smaller insolation during winter at 6 ka, a pronounced warming centre is found over Barents Sea in winter in the simulations, which is also supported by the nearby northern Eurasian continental and Fennoscandian reconstructions. This indicates that in the Arctic region, the response of the ocean and the sea ice to the enhanced summer insolation is more important for the winter temperature than the synchronous decrease of the insolation. Correspondence to: Q. Zhang (qiong@misu.su.se)


Introduction
The development of scenarios for possible future climate change using simulations with climate models is a major scientific challenge.The models need to be validated to assure that they are able to simulate observed climate changes.Simulations of past climates allow the evaluation of how models respond to changes in external forcing such as orbital forcing and greenhouse gases.To undertake such evaluations, the Paleoclimate Modelling Intercomparison Project (PMIP) was launched (Joussaume and Taylor, 1995;Harrison et al., 2002), with a focus on two periods, the Last Glacial Maximum (LGM), 21 000 years ago, and the mid-Holocene (MH), 6000 years ago.
The mid-Holocene climate is reasonably well documented by proxy data.The main forcing compared to the late preindustrial period is the insolation change due to the Earth's slowly changing orbit around the sun (Hewitt and Mitchell, 1996;Kutzbach et al., 1996;Vettoretti et al., 1998;Joussaume et al., 1999).During the mid-Holocene, the seasonal cycle of insolation forcing in the Northern Hemisphere was larger than today, with on average 5% more solar radiation in summer, and 5% less in winter (Berger, 1978).As a result of this change of insolation, a warmer summer and a colder winter could be expected over the northern high latitudes at the mid-Holocene compared to the late pre-industrial period.However, both the proxy reconstructions and climate model simulations reveal that the surface temperature not only responded to the solar forcing but also complex processes within the climate system, such as ocean and landsurface feedbacks, were involved (Cheddadi et al., 1997;Wohlfahrt et al., 2004).Climate model predictions of the response to anthropogenic changes in atmospheric composition suggest that the northern high latitudes are particularly sensitive to the radiative forcing, mainly because of two positive feedbacks; changes in the extent of sea-ice cover, and Published by Copernicus Publications on behalf of the European Geosciences Union.
Q. Zhang et al.: Mid-Holocene climate change: model-data comparison changes in the albedo of the land surface as a consequence of changes in snow cover and the extent of forests (IPCC, 2007).Also climate model simulations for the past have shown that sea-ice-albedo feedback is a key factor, which amplifies climate change in northern high latitude Arctic region (Wohlfahrt et al., 2004;Renssen et al., 2005;Braconnot et al., 2007b).Investigation under the double CO 2 scenario in the Coupled Model Intercomparison Projects (CMIP2) shows that the range of simulated polar warming in the Arctic is between 1.5 to 4.5 times the global mean warming (Holland and Bitz, 2003).The observational evidence from the last century is consistent with climate model simulations that include changes in greenhouse gas concentrations (IPCC, 2007).Observations also show an amplified warming in the Arctic, but its magnitude varies depending on the time period analysed.Jones and Moberg (2003) found about a 2.2 times larger warming in the Arctic compared to the global warming when the periods 1861-2000 and 1977-2001 were considered but only a factor of 1.2 over the period 1901-2000.A variety of proxy records provide temporal and spatial information of climate change during the current interglacial, the Holocene (Cheddadi et al., 1997;Prentice and Jolly, 2000;Kim and Schneider, 2004).Compilations for northern high latitudes suggest that Arctic warm events may have occurred in past climate (D'Arrigo and Jacoby, 1993;Jennings and Weiner, 2003;Jiang et al., 2005;Massé et al., 2008).The mid-Holocene is one of the noticed warm periods that is well documented by proxy data and is often used to evaluate how models respond to a change primarily in insolation.Most of these data-model comparisons focus on proxy data rich areas such as southern and mid-latitude Europe (Liao et al., 1994;Harrison et al., 1998;Masson et al., 1999;Prentice et al., 1998;Guiot et al., 1999;Joussaume et al., 1999;Bonfils et al., 2004;Gladstone et al., 2005;Masson-Delmotte et al., 2006;Brewer et al., 2007).The high-latitude region is less often addressed due to the relative lack of proxy data compared to some other regions.The past and present polar amplification of climate change have, however, been recently reviewed in both CMIP and PMIP simulations by comparison with polar ice core data (Masson-Delmotte et al., 2006).Results indicate that PMIP1 atmosphere-only simulations show no consistent temperature response for the polar regions to the mid-Holocene forcing, whereas the PMIP2 coupled atmosphere-ocean climate models systematically simulate a significant mid-Holocene warming both for Greenland and Antarctica, consistent with the polar ice core data based temperature estimates.
Here, we present a model-data comparison over the northern high latitudes (60 • N-90 • N), which is more comprehensive than previous comparisons; we use the latest PMIP database and a new extensive collection of published proxybased reconstructions.In a companion study (Sundqvist et al., 2010) we analyse this data collection in detail, which includes temperature and precipitation reconstructions covering both the mid-Holocene (6000 year BP; 6 ka) and pre-industrial (AD 1500; 0 ka) periods.The purpose of that study is to obtain an overview knowledge of the climate change documented in the paleoclimate records, as well as to analyse these changes seen in the reconstructions.In the present paper, we use both the PMIP1 (http://pmip.lsce.ipsl.fr/)and PMIP2 (http://pmip2.lsce.ipsl.fr/pmip2/,version of 09-30-2008) data bases to perform a model-model and a modeldata comparison.An optimal selection method, using a cost function approach (Goosse et al., 2006), is applied to measure the discrepancy between model results and reconstructions.This model-data comparison method aims to evaluate the capability of the climate models to reproduce the spatial climate response pattern seen in the proxy data, and to identify a few model simulations that most closely resembles the reconstructions.A further goal is to explore some feedback mechanisms responsible for mid-to late Holocene climate changes.
The paper is structured as follows.In Sect.2, the temperature reconstructions for the northern hemisphere high latitudes from our companion paper (Sundqvist et al., 2010) are summarized.The PMIP simulations used in the current study are described in Sect.3, while their temperature change during the mid-Holocene is presented in Sect. 4. In Sect.5, the model data and the reconstructions are compared with the cost function approach, and the analysis for the related feedbacks is presented.Finally, a summary and discussion is given in Sect.6.

Evidence of the mid-to late Holocene climate change in reconstructions
The collection of proxy data that has been used here for comparison with climate model output is discussed in detail in the companion paper (Sundqvist et al., 2010).Here, we summarize some of the main features of the proxy data and point out some main findings in the companion paper.This set of proxy-based reconstructions is obtained through a systematic scan for published calibrated temperature and precipitation reconstructions from the region north of 60 • N, with data for both the 6 ka and 0 ka periods.Altogether, the selected proxy records are more numerous and cover a larger fraction of the entire area over northern high latitudes than in previous model-data comparisons (Wohlfahrt et al., 2004;Renssen et al., 2005;Masson-Delmotte et al., 2006).Our model-data comparison, however, solely deals with temperature because available precipitation records are too few and too uncertain to permit any conclusive analysis.
Figure 1 shows the spatial distribution of the temperature reconstructions used in this study.In total there are 72 proxy reconstructions, including 46 July and August temperature reconstructions, 6 January temperature and 20 annual mean temperature reconstructions.These reconstructions were obtained from 61 sites, mostly over the land areas neighbouring the North Atlantic sector with the highest density over Fennoscandia.There are some proxy sites in Greenland and a few in the Northern Eurasian and North American continents.Four sites over Siberia provide both January, July and annual mean temperature reconstructions.Over Scandinavia there are two sites that provide both January and July temperature reconstructions and one site that provides July and annual mean temperature reconstructions.The majority of reconstructions are from terrestrial archives with a dominance of pollen, diatoms, and chironomids as temperature proxies.A few borehole, ice-core, tree ring and speleothem records are also included.The number of marine proxies is considerably smaller, only five including marine diatoms, alkenones and foraminifera from the North Atlantic Ocean.
Our intention in this model-data comparison is to use as many data as possible from the published proxy archives and available model data.This will help to distinguish a spatial climate pattern in the proxy data that can be directly compared with the model patterns of climate response to forcing.This comparison offers an opportunity to evaluate the physical reliability of the reconstructions from different proxy data types which provide a more rigorous estimate than a single proxy approach.The multi-proxy approach may help to identify possible seasonal biases in the reconstruction estimates.Pollen records are the most abundant temperature proxy type for the time period studied here.Reconstructions from other proxies provide complementary information and hence serve as additional validation data.
When compiling the reconstructions, the uncertainty from the calibration into temperature estimates done by the respective original authors and the uncertainty in temperature change due to internal variability were considered.These uncertainties were combined into a single standard deviation of the uncertainty in each proxy series, which is here used in the cost function.For practical reasons, July and August reconstructions are used here to represent northern summer (JJA) mean temperature and, January reconstructions are used to represent northern winter (DJF) mean temperature in the model-data comparisons.The climate change in mid-to late Holocene, as recorded in the proxy data, is defined as the difference between the 500-year averages for two periods, that is, 6 ka minus 0 ka.More specifically, for the proxy data, this means the difference between the 500-year periods centred on 6000 years BP and 1500 AD.Sundqvist et al. (2010) provide a detailed motivation for this chosen definition of estimated climate change along with a discussion of the uncertainties in this quantity.
According to an unweighed average of all the temperature proxies in Fig. 1, the northern high-latitude region was about 1 • C warmer in summer, 1.8 • C warmer in winter and 2.1 • C warmer in the annual mean temperature at 6 ka, in comparison to 0 ka.The uncertainties of the overall temperature change largely depend on the uncertainty in individual reconstructions and the number of reconstructions.Uncertainties in reconstructed changes at the individual site level range from about 3.2 • C to only 0.2 • C (Sundqvist et al., 2010).The sparsity of the data distribution and the large uncertainty in winter temperature make the estimated change (warmer at 6 ka than at 0 ka) in winter temperature less reliable than the estimated changes in summer and annual mean temperatures.It has previously been noted, though, that in Europe warmer winter temperatures at 6 ka compared to 0 ka have been found in a majority of PMIP1 models and in proxy-based values (Masson et al., 1999).Such higher winter temperatures at 6 ka cannot be explained by the insolation change directly, because the winter insolation in the high latitudes was lower in the mid-Holocene than late Holocene.
Another feature pointed out in the companion paper, is the observed larger change in the annual mean temperature change compared to winter and summer temperature changes in the reconstructions.One reason for this behaviour could be that the proxy data used to reconstruct the annual temperature are not from the same archives as those for the summer and winter temperature, resulting in different locations and numbers of the data.Moreover, we cannot exclude the possibility that the response in seasons for which we have no explicit data (spring and autumn) may have larger contributions to the change in annual mean temperature than the change in winter and summer.Furthermore, one can also not exclude the possibility that the transfer functions used to estimate the temperatures from proxy data are not sufficiently robust to allow a quantitative comparison between estimates for different seasons.Usually the reconstructions are not well Table 1.PMIP models used in this study.The model name is listed as specified in the PMIP database.Type A refers to an atmosphereonly model with fixed SST, OA refers to a coupled ocean-atmosphere and OAV refers to a coupled ocean-atmosphere-vegetation model.OA + OAV means the model has both simulations.The resolution of the spectral models is given by the type of truncation, the highest total wavenumber, and the number of vertical levels, while for the grid models it is given in longitude (degree) × latitude (degree) × vertical levels.Yukimoto et al. (2006) resolved seasonally.Therefore, the problem needs to be investigated further with the aid of climate model simulations that have reasonable physical constraints and appropriately prescribed external forcing corresponding to the particular past periods.In the following section first simulations from the PMIP database to obtain an estimate of the simulated climate response to the change in forcing will be analysed and then these will be compared the reconstructions.

The PMIP mid-Holocene simulations
The PMIP simulations used in this study include nineteen atmosphere-only models with fixed SST (SSTf) from the PMIP1 database, thirteen coupled ocean-atmosphere (OA) models and six coupled ocean-atmosphere-vegetation (OAV) models from the PMIP2 database.The models are presented with their names as specified in the PMIP1 and PMIP2 databases.Their spatial resolutions and references to papers that describe the models are listed in Table 1.For most of the modelling groups, the version of the coupled OA or OAV models used for PMIP2 is identical to the version used for future climate change predictions in CMIP3, but with a lower resolution.Each model has been used to perform a mid-Holocene simulation (6 ka) and a pre-industrial control simulation (0 ka), under the same external forcing as required by the PMIP protocol.The main difference in forcing for mid-Holocene compared to 0 ka, is defined by the orbital parameter, which is represented by the eccentricity, obliquity and precession.The PMIP protocol also considers a change in the atmospheric CH 4 concentration with a lower concentration at 6 ka.The concentration for other greenhouse gases and the topography are the same in the 6 ka and 0 ka simulations (Braconnot et al., 2007a) for each model.The topography can differ between the models due to their different horizontal and vertical resolution.
In PMIP1 simulations the models were integrated at least for 11 years and the last 10 years of each simulation were stored in the database.In PMIP2 simulations the models were run for at least 100 years after the coupled model reached equilibrium and the last 100 years of the simulation were stored in the database.Hence, all results presented here are computed from 10-year averages for PMIP1 simulations and 100-year averages for PMIP2 simulations.To compare with the reconstructions, rather than using single months (e.g.July or January), we calculate boreal seasonal averages for the models, i.e.December-January-February (DJF) mean for winter, March-April-May (MAM) mean for spring, June-July-August (JJA) mean for summer, and September-October-November (SON) mean for autumn.The annual mean is computed from the twelve monthly means.
In the present study the temperature from the model outputs refers to the surface air temperature at 2 m height.The analyses are focused over the high latitudes (north of 60 • N).The simulations for 0 ka conditions are regarded here to be control runs.The climate response to the change in insolation from the mid-to late Holocene is defined as the mean climate change between two time periods, i.e. 6 ka minus 0 ka.Hence, for example, when we speak about a "warming" response here, we mean that the simulation for the 6 ka period is warmer than the simulation for the 0 ka period.The three different types of models can help to identify the different responses in the climate system; the PMIP1 atmosphere only model simulations mainly show the direct atmospheric response; the PMIP2-OA simulations introduce the feedback from both ocean and sea ice; and the PMIP2-OAV simulations introduce further feedback from vegetation.The 6 ka orbital configuration leads to an increase of the amplitude of the annual cycle of the incoming solar radiation at the top of the atmosphere in the northern hemisphere and a decrease in the southern hemisphere.The insolation averaged over northern high latitudes (north of 60 • N) shows a slight increase of about 2.9 W/m 2 for the annual mean, but a large increase by 23.5 W/m 2 for the JJA mean with a maximum increase of about 32 W/m 2 in July (Fig. 2).The DJF mean insolation shows a slight decrease by about −2.3 W/m 2 over the region.In the simulations the modern calendar has been used for both periods instead of using a celestial based calendar as suggested by Joussaume and Braconnot (1997).The consequence is that a simulated climate response that occurs in autumn will be slightly underestimated in the northern hemisphere (Braconnot et al., 2007a).The reduction in CH 4 results in 0.07 W/m 2 decrease in radiative forcing (Otto-Bliesner et al., 2006), therefore the mid-Holocene climate response is mainly due to the changed orbital configuration.

Temperature change in PMIP simulations
The surface air temperature response to the insolation is illustrated in Fig. 3 for the three types of PMIP models.In summer, all the PMIP1 and PMIP2 simulations show a warming at high latitudes in response to the enhanced summer insolation (Fig. 3c).In the PMIP1 atmosphere-only simulations, the atmospheric response produces a summer warming between 0.3 to 1.6 • C; the average for nineteen simulations is 0.84 • C (Table 2).In the PMIP2 OA and OAV simulations, the response is more consistent (the spread is between 0.8 and 1.6 • C) across the models than for PMIP1, except the OAV simulation with MRI-CGCM2.3.4nfa which has an anomalous 3.0 • C response.This suggests that a dynamic response of the ocean and sea-ice in the PMIP2-OA narrows the inter-model spread of the summer warming compared to the PMIP1 simulations, for which the modern SST and sea ice fraction are prescribed.The averaged summer temperature change for the thirteen PMIP2-OA simulations is 1.10  2).Since most of the Arctic region is covered by sea ice, this difference between PMIP2-OA and PMIP1 simulations is primarily due to the reduced sea-ice cover, which induces the well known sea-ice-albedo positive feedback.This positive feedback appears to be robust across all the thirteen PMIP2-OA simulations as indicated by the small spread between models (Fig. 3c, middle panel).The average for the five PMIP2-OAV simulations shows a response which is stronger by 0.20 • C than the response in the PMIP2-OA simulations (Table 2), reflecting the amplification due to the vegetation feedback.
In accordance with the decreased insolation at 6 ka in the other three seasons (Fig. 2), the direct atmospheric response shows cooling in spring, autumn and winter in most of the PMIP1-SSTf simulations with averaged values of −0.51 • C, −0.12 • C and −0.33 • C, respectively (Table 2).However, one or a few models with warming are seen in all these three seasons.As a result, there is no obvious change in annual mean temperature in the PMIP1-SSTf simulations.The feedback from the sea-ice does not appear to have any importance in spring, for which the atmospheric response to the insolation appears to be dominant.No distinct differences in spring temperature response are found between the three types of PMIP simulations (Fig. 3b).Most of the PMIP1 and the PMIP2-OA simulations show cooling and the PMIP2-OAV simulations show a slight cooling or no change at all in spring.An exception is for the MRI-CGCM2.3.4nfa-OAV which shows a large warming.
The influence of ocean and vegetation feedbacks on the surface air temperature appears to be more robust in autumn and winter.The average autumn warming in the PMIP2-OA simulations is 1.35 • C, which is about 0.25 • C larger than the summer on average, indicating a lagged temperature response of the ocean to the enhanced summer insolation (Fig. 3d).The additional feedback from the vegetation in the PMIP2-OAV simulations further increases the autumn warming to 2.0 • C (Table 2, the anomalous MRI-CGCM2.3.4nfasimulation is not included in the average).
In winter, ten of the thirteen PMIP2-OA simulations and all the six PMIP2-OAV simulations show a warming (Fig. 3a), whereas most of the PMIP1 simulations show a cooling (on average −0.33 • C).The winter warming is apparently caused by the ocean and the sea ice feedback is 0.88 • C, and the vegetation feedback yields another 0.67 • C warming.Compared to the 0 ka control simulation the average winter warming reaches an average of 1.22 • C in the PMIP2-OAV simulations.We noticed that the flux-adjusted and nonfluxadjusted version of MRI-CGCM2.3.4 show quite different response in winter.In Fig. 3a MRI-CGCM2.3.4fashows about 0.5 • C warm, but MRI-CGCM2.3.4nfashows about −0.6 • C cold.The cold response also shows much more sea ice cover in winter in MRI-CGCM2.3.4nfa-OA, as shown in Fig. 10 in Braconnot et al. (2007b).It implies the fluxcorrection is crucial in this model at least for Northern highlatitudes.From Table 2 and the discussion above, it can be concluded that the combined effects of orbital forcing, ocean, sea-ice feedback and vegetation feedback produce about 1.3 • C warming in summer, 2.0 • C warming in autumn and 1.2 • C warming in winter, but no clear change in spring.The important feedbacks from ocean, sea-ice and vegetation take effect from summer to winter, but their effect is most important in autumn and winter.Together, they eventually translate into an annual mean temperature response of about 1.1 • C. Strictly, one cannot separate the relative contribution from ocean or from vegetation simply by calculating the difference between the PMIP2-OAV and PMIP2-OA, because for a given model, the OA and OAV simulations do not share the same control simulation.The possible feedbacks between ocean and vegetation might not be separated as well.
In the above overall comparisons the response varies from model to model in autumn and winter, i.e. when the ocean and vegetation feedbacks take a notable effect (Fig. 3a and  d).This can be explained by the different ocean, sea-ice and vegetation physics that are applied in the various PMIP2 models.
The most consistent response, across the models, is seen for the summer temperatures in the PMIP2-OA simulations (Fig. 3c).All thirteen simulations show a warming between 0.8 and 1.6 • C, with an average of about 1.0 • C. We selected the six simulations with summer temperature responses that Corresponding to the direct response to the insolation in summer the warming over the continents is more pronounced than over the ocean.This feature is also the main pattern seen in temperature change in the thirteen simulation ensemble (not shown).However, the warming centres differ from model to model.For example, CCSM and CSIRO-Mk3L-1.1 have their maximum warming over north America and north Eurasia, as well as Greenland, whereas in MRI-CGCM2.3.4fa the warming over Greenland is small.In the GISSmodeE and MIROC3.2simulations, there is a small cooling area in the north Pacific and north Atlantic, respectively.

Model-data comparisons
Since the available reconstructions do not cover the entire high latitude region, a model-data comparison could be obtained merely over the areas covered by data.We present such a comparison for the averaged temperature change over the locations where the proxy reconstructions are available.
To this end, the corresponding model grid points for 46 summer temperature, 6 winter temperature and 20 annual mean temperature reconstructions have been selected in the ensemble means of PMIP1-SSTf, PMIP2-OA and PMIP-OAV data.
Table 3 shows the average of the temperature change in three types of PMIP simulations and in the proxy reconstructions.Generally, the simulated responses in summer, winter and annual mean temperatures at the proxy locations (Table 3) are quite similar to the responses for the entire high latitudes (Table 2); that is, the more complex models have a stronger response.In the OA and OAV models, the magnitude of the change in the simulated annual and winter mean temperatures at the proxy locations are smaller than the average over the entire region.This may indicate that the annual and winter mean response over the Arctic Ocean is larger than over land, which is where the most reconstructions are located.
In summer, the response is about the same in proxy covered areas as in the entire region.The averaged winter and annual temperature change in the PMIP2-OAV ensemble is closer to the reconstructions than PMIP2-OA and PMIP1-SSTf.In summer, the results are less conclusive in this respect.This may partly be due to the different number of simulations included in each ensemble.However, the ensemble mean for five PMIP2-OAV simulations also has a stronger response than their PMIP2-OA counterpart five-member ensemble mean.Furthermore, the five-member PMIP2-OAV ensemble mean is also closer to the reconstructions in winter and annual mean temperature change, though not in summer.
A map of the ensemble-mean temperature changes may provide important information that reflects the inter-model consistency, but could also miss some regional features by averaging across the different model simulations.To identify simulated regional features of the climate response that are likely to be realistic, it seems preferable to select one (or a few) simulations that most closely resemble the pattern seen in the proxy-based reconstructions.To achieve this goal, we applied a variant of the model-data comparison technique that was developed by Goosse et al. (2006).By using this technique, we can select among the PMIP simulations the ones that most closely resemble the reconstructions.

Selection of an optimal simulation
The principle of the method used by Goosse et al. (2006) is to select among a relatively large ensemble of simulations from one climate model, the one that is the closest to the observed climate.We have slightly adapted this technique to suit our purposes.Our goal is to find the simulation that objectively shows a best fit with the available proxy data.The selec-tion is performed by comparing each simulation to the available reconstructions in a consistent manner.This method is similar to data assimilation or the identification of analogues in meteorology and is expected to result in a good fit between model simulations and reconstructions at a reasonable cost.An important difference between our application of the method and the approach taken by Goosse et al. (2006), is the fact that we applied the objective selection method to an ensemble of model simulation with different models as opposed to an ensemble of simulations with a single model.
To select the best-fit simulation from the PMIP database, we applied a cost function where CF k is the value of the cost function for each PMIP simulation k.In our study, CF k is calculated for summer, winter and annual mean temperature.The quantity n in the summation, is the number of reconstructions used in the model-data comparison.T rec,i , is the temperature change for reconstruction i, at a particular location.T k mod,i is the value of the corresponding temperature change in the PMIP simulation k for the model grid box that contains the location of the proxy-record i. w i is a weight factor.We have compared results for two definitions of w i , to test the influence of the choice of the weight on the cost function.First we give the same value to all the weights w i in the computation of CF.In this case, w i is thus equal to 1 n , where n is the number of reconstructions.Secondly we allow w i to represent the uncertainties of the reconstructions to ensure that the reconstructions with larger uncertainties contribute less to the cost function.Here we define the weights as w i = 1 σ 2 i + 1 where σ i is the estimated uncertainty for reconstruction i, that is, the combined uncertainty due to calibration uncertainty and internal variability (Sundqvist et al., 2010).To avoid that a few reconstructions with unrealistically small estimated uncertainties dominate the value of cost function, a constant error of 1 is added in the denominator.
We also estimate the model uncertainty by analysing data from MRI-CGCM2.3.4fa-OA, which provides one 50-year simulation and one 100-year simulation.We find that the variance for the annual mean temperature between the 50year simulation and the 100-year simulation is 0.004 • K. Compared to the uncertainties in the reconstructions the uncertainties in this model can be neglected in the calculation of the cost function.In the absence of the possibility to undertake similar comparisons also for the other models, we assume that the model uncertainty is in general small compared to the reconstruction uncertainty.
The relative magnitude of CF represents the goodness of each simulation compared to other simulations.A low value of CF means that the simulated climate response pattern is close to the temperature change seen in the reconstructions, and the corresponding model is regarded to reproduce the reconstructed climate change better than the other models with higher CF values.We calculated the cost function for the zero change in models.It means that the model has no sensitivity to orbital forcing at all.One can speculate that the cost function would be very high since the model has no response.The results are not out of our expectation.In the case of T k mod,i = 0, the cost functions are high in all the seasons for all the models.
To evaluate the overall goodness-of-fit between the reconstructions and the three types of PMIP simulations, the CF values of the ensemble means for nineteen PMIP1-SSTf simulations, thirteen PMIP2-OA simulations and five PMIP2-OAV simulations are shown in Table 4. Here, we have calculated the ensemble mean temperature change for each type of PMIP simulations and then for each ensemble mean calculated the cost function based on 46 summer temperature, six winter temperature and 20 annual mean temperature reconstructions respectively.CF is determined separately for the equal weight w 1 and the weight w 2 including the uncertainty.Table 4 shows that w 2 gives smaller CF values than w 1 , indicating the importance of taking the reconstruction uncertainty into account.For winter and annual mean data, the CF values for the PMIP2-OAV ensemble are smaller than those for PMIP2-OA, which are in turn smaller than those for PMIP1-SSTf.This indicates that by including more complexity in the models, a better agreement with the proxy records, in terms of winter and annual mean temperatures, can be achieved.For summer, however, the CF values are essentially equal for all three types of models.Nevertheless, if we consider the overall performance across all seasons studied, the PMIP2 models are in closer agreement with the reconstructions than the PMIP1 models, indicating a larger degree of realism in the PMIP2 simulations.Therefore, we will in the following calculate the CF values for individual PMIP2 simulations (but not PMIP1) and select the models that give the smallest CF.
The CF values for each PMIP2 simulation are shown in Fig. 5, separately for the summer, winter and annual mean temperature changes.The OAV simulation by MRI-CGCM2.3.4nfa is excluded since its CF value for summer temperature is three times larger than the summer CF value average.We have tested for both the equal weight w 1 and uncertainty considered weight w 2 .Similar to the CF values in Table 4, we consistently find smaller CF values for all models and all seasons, when taking the uncertainties into account.This effect is largest and most clearly observable for winter data.Moreover, uncertainty weighted CF values have smaller spread among the models.When using w 2 , the CF values for the large majority of PMIP2 models are higher for the annual mean temperature changes than those for summer and winter temperatures (Fig. 5).The larger change seen in the annual mean reconstructions compared to those for summer and winter might indicate inconsistencies in the reconstructed changes of annual mean temperatures (Sundqvist et al., 2010).Consistent with the overall comparisons, most PMIP2 models have smaller CF values in summer than in winter and annual mean.The CF values for summer also show less variability across the models.We interpret this as an indication that the PMIP2 models consistently simulate well the climate response to summer insolation forcing, particularly for the Atlantic region where the most reconstructions are located.The winter CF values in these simulations are more variable.For example, the CF for GISSmodelE-OA, MIROC3.2-OA and UBRIS-HadCM3M2-OA have notably high CF values in winter.This behaviour, however, probably reflects uncertainty also due to very few winter proxy records being used, rather than enhanced uncertainty in the three models mentioned.
The aim of the above described objective method for model-data comparison is not to rank the models, but to identify few models that are the closest to the available reconstructions, while accounting for the estimated uncertainty in the latter.A good fit (i.e.small CF) between a model (few models) and proxy evidence is interpreted to indicate similar responses to external forcing.Hence, the "best-fit" models should be more suitable for studying the dynamical processes than the multi-model ensemble mean or single models with large CF values.According to the CF values, however, the goodness-of-fit for the selected simulations are different in summer, winter and in the annual mean.It is therefore not easy to identify a single model simulation that shows a best fit with the proxy data.A possible way to select the few best-fit simulations is thus to consider the ones that have the relatively smaller CF value in winter, summer and annual mean.Based on CF values in Fig. 5, two OAV models (FOAM-OAV and MRI-CGCM2.3.4fa-OAV)show overall lowest CF values.Their OA counterparts, FOAM-OA and MRI-CGCM2.3.4fa-OA, also show relatively low CF values compared to other PMIP2-OA models.In order to focus on some important feedback mechanisms, we select FOAM-OA and MRI-CGCM2.3.4fa-OAsimulation for our mechanism study which is described in the next section.In the following section, MRI-CGCM2.3.4fa-OA is referred in a short name Fig. 5. Value of the cost function for the 18 PMIP2 models for summer temperature (red), winter temperature (green), annual mean temperature (blue).The value of the cost function is normalized by the number of the reconstructions.The horizonal X-axis indicates the 18 PMIP2 models used in the comparison.Number 1 to 13 are PMIP2-OA models: 1, CCSM; 2, CSIRO-Mk3L-1.0;3, CSIRO-Mk3L-1.1;4, ECBILTCLIOVECODE; 5, ECHAM5-MPIOM1; 6, ECHAM53-LPJ; 7, FGOALS-1.0g;8, FOAM; 9, GISSmodelE; 10, IPSL-CM4-V1-MR; 11, MIROC3.2;12, MRI-CGCM2.3.4fa; 13, UBRIS-HadCM3M2.Number 14 to 18 are PMIP2-OAV models: 14, ECBILTCLIOVECODE; 15, ECHAM53-LPJ; 16, FOAM; 17, MRI-CGCM2.3.4fa; 18, UBRIS-HadCM3M2.
as MRI-OA.The reason why we prefer the OA versions, and not the OAV versions, is that we will focus on the ocean and sea-ice-albedo feedbacks without any couplings to feedbacks from the vegetation models.

Large-scale change in surface temperature and related feedbacks
Figure 6 illustrates the change in surface temperature in the two selected PMIP2-OA simulations and in the reconstructions.For the simulated change in summer temperature, the main feature in both simulations is the remarkable land-sea contrast; the warming over the Eurasian continent, North America and most of Scandinavia is much higher than over the ocean.Over the continents, the magnitude of the summer warming in FOAM-OA is much stronger than in MRI-OA, but it is mostly closer to the magnitude of the reconstructions, especially those over Siberia and Scandinavia.The reconstructions in the data-rich Fennoscandian region display a spread with majority showing a warming and only two records showing cooling.One record over north Finland shows minor cooling (less than −0.1 • C).One chironomid-based temperature reconstruction from northern central Sweden rise abruptly in late-Holocene, probably as a result of local limnological changes in lake Gilltjärnen and its catchment (Antonsson et al., 2006).
Despite similar low CF values, the simulated changes in winter in the two simulations differ rather much from each other.This is not unexpected given that all winter temperature proxy sites are located over only a small region including northern Scandinavia and northern Eurasia.This geographical spread of proxy sites is too small to constrain the model temperature response over much of the rest of the Arctic region.FOAM-OA shows a winter warming over most of the land areas, except southern Greenland, while MRI-OA shows notable cooling over much of North America and, to a smaller degree, over the southern parts of the Eurasian Arctic region.A striking feature in the simulated winter temperature response is a strong warming region centred over Barents Sea in both simulations.The geographical extent of the strongest warming is larger in FOAM-OA compared to the other model; including in FOAM-OA also northernmost Fennoscandia and the north-western Eurasian continent.There are no winter temperature reconstructions available for the ocean, but the warming centred over Barents Sea seen in the simulations is geographically connected with the warming seen over Siberia and northern Fennoscandia.The latter is in agreement with five of the six reconstructions, and is, in fact, the reason why these two simulations have low CF values for winter temperatures (Fig. 5).
The simulated change in annual mean temperatures in MRI-OA is characterized by a warming along the eastern sector of the Arctic Ocean and Eurasian continent, with the warm centre located over Barents Sea, indicating a major contribution of the winter response to the annual mean.FOAM-OA also shows a warming over the same region, but in this model the warming is more widespread over northeast Eurasia, the Arctic Ocean and even North America.In contrast, MRI-OA shows cooling in much of north America and only a small warming over the Arctic Ocean.The distribution of annual proxy records with the additional coverage obtained over Greenland is only slightly better than in winter.Except two data points over North Fennoscandia, most annual reconstructions show warming, in agreement with the same locations in MRI-OA and FOAM-OA.Because of the large areas without any proxy data, i.e. in North America and the marine areas, the simulated spatial pattern of the annual mean temperature change is not well constrained, but slightly better than for winter data.
To explore possible mechanisms behind the simulated temperature changes in the FOAM-OA and MRI-OA models, we show the simulated seasonal variations of the sea ice fraction, the sea ice thickness, the snow cover fraction over land, the surface albedo, the downward ocean surface heat flux and the upward ocean heat flux at the ice base (Fig. 7).It can be seen that, following the enhancement of the insolation in summer, the sea ice over the northern high latitudes is reduced all year around and has a maximum decrease by about 25% in MRI-OA and about 8% in FOAM-OA in August and September (Fig. 7a).As shown in Fig. 10 in Braconnot et al. (2007b), most of the PMIP2-OA   models show considerable decrease in sea ice cover during summer.In addition to the decrease in sea ice cover, the sea ice thickness also decreased all year round, with a maximum decrease of about 1.5 m in September in MRI-OA and 0.8 m in FOAM-OA (Fig. 7b).The maximum decrease in snow covered area over land also occurs in September, with a magnitude of 4% in FOAM-OA and 1.2% in MRI-OA (Fig. 7c).The combined effect from the reduction of the sea ice cover and snow cover leads to a decrease of the surface albedo all year round, with a maximum decreases by about 0.03 and 0.04 for FOAM-OA and MRI-OA, in August and September, respectively (Fig. 7d).The surface ocean receives more heat during May to August directly from the enhanced summer insolation (Fig. 7e).However, the two models behave quite differently in absolute values, as MRI-OA shows a July increase by about 20 W/m 2 and FOAM-OA only by 5 W/m 2 .Due to the reduced sea ice and surface albedo, ocean heat storage during summer is amplified.In September, following the decreased insolation over the northern high latitudes (Fig. 2), the ocean starts to release the heat to the atmosphere.Therefore, the insolation as well as the positive seaice albedo feedback result in more heat being released from the ocean to the atmosphere in autumn and winter.Figure 7 shows that this effect is more pronounced in MRI-OA than in the FOAM-OA.In summer, enhanced insolation and decreased ice thickness cause an increased downward heat flux transfer to the ocean through the ice in MRI-OA over the ice covered regions.In winter when the ocean releases heat through the ice to the atmosphere, an increased upward heat transfer is observed in MRI-OA over the ice covered regions (Fig. 7f).These simulated changes of the sea-ice, surface albedo and ocean surface heat flux, which are ultimately forced by the increased summer insolation, contribute to the significant warming during autumn and winter in the Arctic region.The feedbacks are broadly consistent in the two simulations, but the magnitude of the feedbacks resulting from the ocean response are more pronounced in MRI-OA than in FOAM-OA.This is probably due to the polar filtering employed in the FOAM ocean model OM3, which yields a poor performance at high latitudes (Pierrehumbert and Liu, personal communication, 2009).This problem also affects the high latitude ocean heat transport in FOAM-OA, which has no significant ocean heat transport polewards of 60 • N. In MRI-OA, on the other hand, the enhanced summer insolation yields a significant increase in the annual-mean northward heat transport in the North Atlantic Ocean, about 0.01 PW at 60 • N.
Figure 8 shows the spatial pattern of the summer sea ice area fraction at 6 ka, 0 ka and for the change between the two periods in the two selected simulations.In MRI-OA, a large decrease in sea ice fraction occurs over the Eurasian sector of the Arctic Ocean, centred over the Barents Sea.At 6 ka, around 50% of the Eurasian Arctic Ocean becomes icefree during summer (Fig. 8).In FOAM-OA, both at 6 ka and 0 ka, the simulated sea ice area is relatively large and covers almost the entire region north of 70 • N. The reduction of sea ice in FOAM-OA occurs over the entire region, but most strongly along the land-sea boundary, with a maximum reduction nearby the Barents Sea.The anomaly pattern in sea ice thickness, ocean surface heat flux and surface albedo also essentially follows the anomaly pattern of sea ice (not shown).The significant reduction in sea ice during summer together with the well-known sea-ice-albedo positive feedback amplifies the summer insolation forcing, and causes a larger temperature response in summer and autumn (Fig. 3 and Table 2).Simultaneously, the ocean receives more insolation during summer over the regions with reduced ice cover, which further warms the ocean subsurface.More heat is thus stored in the upper ocean in the following season (Fig. 9).The maximum warming response in the upper ocean occurs in September, in consistency with the much enhanced high-latitude temperature response in PMIP2-OA compared to PMIP1-SSTf models (Fig. 3 and Table 2).The warmer upper ocean leads to less formation of ice in winter and reduced sea ice thickness (Fig. 7b).This, in turn promotes the ocean to release more heat to the atmosphere over the regions with decreased ice cover (Fig. 7e).The increased upward heat flux through thinner ice also contributes to the enhanced winter time surface heat flux to the atmosphere over the Arctic Ocean (Fig. 7f).The pattern of the change in simulated DJF mean sea level pressure (SLP) in FOAM-OA (Fig. 10a) and in MRI-OA (Fig. 10b) is consistent with the winter warm conditions that shown in Fig. 6.A similar feature in the two simulations is the decreased pressure over the polar region and the increased pressure over part of southern Europe and central Asia.In FOAM-OA, the mid-latitude Atlantic and North America also exhibit increased SLP, indicating a mean shift towards a more positive phase of the North Atlantic Oscillation (NAO) and even the Northern Annular Mode (NAM) in this model.This feature has previously been addressed by Gladstone et al. (2005) in a comparison of the NAO in PMIP2-OA simulations.Most models show no significant SLP change in terms of its mean and its interannual variability.There is, however, a tendency towards an increased meridional SLP gradient in PMIP2-OA as well as in PMIP2-OAV models.The stronger pressure gradients over the north Atlantic region may have driven the storm track further north at 6 ka, resulting in warm conditions over the north Atlantic and northern Fennoscandia (Thompson and Wallace, 1998).However, in MRI-OA, a decreased pressure over the midlatitude Atlantic does not lead to such a positive shift in NAO.The warm centre over the Barents Sea in both simulations is more related to the strengthened Siberian high and the weakened Arctic low.Corresponding strengthened SLP gradients between the continent and the Arctic Ocean should have enhanced the northward inflow of warm air into the Arctic region.

Summary and discussion
We have performed model-model and model-data comparisons to examine the climate response to the change in insolation between the mid-and late Holocene.Between 6 ka and 0 ka, the altered orbital forcing leads to an increase by 23.5 W/m 2 of the top of the atmosphere insolation over northern high latitude in summer, and a slight decrease by −2.3 W/m 2 in winter over the same region.The PMIP1 simulations with fixed SST show that the atmospheric response to this orbital forcing produces on the average a 0.8 • C warm- ing in summer and a cooling in the rest of the year.In PMIP2 ocean-atmosphere coupled simulations, the sea-icealbedo feedback enhances the summer warming to 1.1 • C, and the thermal inertia of the ocean leads to a 1.4 • C warming in autumn and a 0.6 • C warming in winter, while the cooling in spring remains the same as in the PMIP1 simulations.The PMIP2-OAV simulations that include vegetation feedbacks also show warming in summer, autumn and winter; the changes being beyond 1.0 • C in winter and summer, and reaching 2.0 • C in autumn.When comparing these results with the temperature changes seen in reconstructions from all available temperature proxy data (for summer, winter and annual mean), the results from the PMIP2-OAV simulations most closely resemble the reconstructions.This indicates that when feedbacks from the ocean, sea ice, and vegetation are included in the models, the simulated climate response is generally in better agreement with the proxy-based reconstructions.
Based on an objective cost function approach we selected two PMIP2-OA models that showed the best overall agreement with the available proxy data.The two models show a spatial response pattern that is largely consistent with the reconstructions in summer.The simulated summer temperature change follows closely from the insolation change and shows strong warming over land and relative weak warming over ocean.Hence, simulated summer warming is seen everywhere in the high-latitude regions.The large majority of proxy-based reconstructions also shows a warming in summer, but a few records suggest a cooling.The consistent simulated warming and dominance of warming in the summer temperature proxies together suggest that the few proxies that indicate cooler summers at 6 ka compared to 0 ka should be questioned in terms of their long-term changes.For winter and the annual mean fields, however, the geographical distribution of the proxy sites is too small to adequately constrain the models across the entire study region.Because of this, we find rather large differences between the response pattern seen in the two selected models for winter and annual mean data.A pronounced warming center is found over Barents Sea in winter in both models, that is supported by the few available nearby reconstructions from northern Fennoscandia and northern Eurasia.Winter warmth over the Barents Sea is also found in warm periods in the last millennium, such as the 15th century Arctic warming and early 20th century warming (Goosse et al., 2003;Crespin et al., 2009).
Furthermore, analysis of the two selected best-fit OA models shows that in summer and early autumn, the reduction in sea ice coverage and snow coverage lead to a decrease in surface albedo, which together with ocean feedbacks, leads to about 1.2 • C increase of summer temperature.During summer, a reduced sea ice extent and a lower surface albedo as well as a reduction in sea ice thickness enhance the warming over the ocean.Furthermore, there is an increased transfer of heat from the atmosphere to the ocean, which in turn results in more heat being stored in the upper ocean.During autumn winter more heat is released from the ocean to the atmosphere that causes a local decrease in sea level pressure and increases the pressure gradient over the region.This enhances the northward flow of warm air into the region, and moreover results in a significant winter warming despite the weaker insolation.This result indicates that in the Arctic region, the response of the ocean and the sea ice to the enhanced summer insolation and the associated change in atmospheric pressure patterns and heat transport are more important for the winter temperature than the synchronous decrease of the insolation.
It is evident that the few proxy-data reconstructions for winter and annual mean temperatures cause large uncertainties in the model-data comparisons.Therefore, caution must be taken when attempting to attribute causes (in reality) for the warming response seen in the models.However, if this is also the true climate response, the exploration for more proxy data is needed to test the reliability of the measured response and for the examination of the possible mechanisms behind the response.
In many previous model-data comparison studies (Cheddadi et al., 1997;Masson et al., 1999;Prentice et al., 1998;Guiot et al., 1999;Bonfils et al., 2004;Masson-Delmotte et al., 2006;Brewer et al., 2007) the climate variables considered are not summer, winter and annual mean temperature, but bioclimatic variables such as growing degree days, temperature of the coldest month and a humidity index that have been identified to better reflect the plant physiology.Most of the temperature reconstructions we used in this paper are indeed from bioclimatic proxies such as pollen, and this type of proxy mainly provides information of climate in summer.Even if an individual proxy type may theoretically be more affected by indices such as growing degree days, it is a fact that the target for calibration chosen by the original proxy data investigators has often been seasonal mean temperatures.We have chosen to study seasonal mean temperatures to make it possible to obtain and compare a large number of proxy series.If we had chosen to study different target climate indices for different proxies, then it would not have been possible to include as many proxy series in the analysis as now.It may be that studies of more "optimal" calibration targets for different proxies could increase signal-to-noise ratios, but such a study has been beyond our scope.Our choice instead has the advantage of making model-data comparisons of temperatures possible for both summer and winter.However, since few data outside the summer season are available, especially winter, it is hard to get a complete picture of the climate response.For example, the change in annual mean temperature in proxy reconstructions is not equivalent to the average for summer and winter temperature; it is even larger than the summer and winter temperature change in our case.The monthly mean results from climate models can help us to explain why this happens.For instance, the strongest warming response in the models are found in autumn.Because of feedbacks in the climate system, such a response could explain the large response seen in annual mean temperature proxy data.Hence, model studies can help to clarify the nature of the climate signals seen in proxy data.

Fig. 1 .
Fig.1.Locations of the temperature reconstructions from proxy data from high latitudes (60 • N north).There are 46 summer temperature, 6 winter temperature and 20 annual mean temperature reconstructions, from in total 61 sites.Red circles denote proxies for summer temperature only.Green circles: proxies for the annual mean temperature.Blue circles: sites with proxies for both summer, winter and annual temperature.Orange circles: summer and winter temperature.Purple circles: summer and annual mean temperature.

Fig. 2 .
Fig. 2.Incoming solar radiation (W/m 2 ) at the top of the atmosphere averaged over 60 • N-90 • N for the mid-Holocene (6 ka, red), pre-industrial (0 ka, blue) and the difference between the two (6 ka-0 ka, black).The left-hand vertical axis refers to the absolute values and the right-hand scale to the difference values.

Fig. 3 .
Fig. 3. Seasonal change in surface air temperature ( • C) averaged over northern high latitudes (60 • N-90 • N) for the PMIP simulations.(a) DJF mean, (b) MAM mean, (c) JJA mean and (d) SON mean.The red circles are for PMIP1, the blue ones for PMIP2-OA and the purple ones for PMIP2-OAV.Note that the changes are computed with respect to the control simulation of the respective model, that is, 6 ka minus 0 ka.

Fig. 6 .
Fig. 6.Large scale pattern in surface temperature( • C) change in FOAM-OA (left column), MRI-OA (middle column), and reconstructions (right column).Top row is for summer temperature, represented by JJA mean for model data and July and August temperature for reconstructions; middle row is for winter temperature, represented by DJF mean for data and January temperature for reconstructions; bottom row is for annual mean temperature.

Fig. 6 .
Fig. 6.Large scale pattern in surface temperature ( • C) change in FOAM-OA (left column), MRI-OA (middle column), and reconstructions (right column).Top row is for summer temperature, represented by JJA mean for model data and July and August temperature for reconstructions; middle row is for winter temperature, represented by DJF mean for model data and January temperature for reconstructions; bottom row is for annual mean temperature.

Fig. 7 .
Fig. 7. Change in seasonal variation of (a) Sea ice area fraction (%), (b) Sea thickness (m), (c) Snow area fractions over land (%), (d) Surface albedo, (e) Downward ocean surface heat flux (W/m 2 , positive means ocean receive heat from atmosphere), and (f) Upward ocean heat flux at the ice base (W/m 2 , positive means ocean release heat to atmosphere).The area average is calculated for the 60 • N-90 • N latitudes in FOAM-OA (red curve) and MRI-OA (blue curve).

Table 3 .
Summary of seasonal changes in temperature ( • C) averaged over the locations of available reconstructions.Ensemble means of nineteen PMIP1-SSTf simulations, thirteen PMIP2-OA simulations and five PMIP2-OAV simulations (The MRI-CGCM2.3.4nfa-OAV is excluded) are given, as well as an ensemble mean for five PMIP2-OA simulations to compare with their OAV counterparts.