Climate of the Past Title : Statistical downscaling of a climate simulation of the last glacial cycle : temperature and precipitation over Northern Europe

Comment: All have a sort of “map”, the function of latitude and longitude sj(x, y). For the method to be valid, any contribution of these maps to the final result needs to be fully independent of time: it needs to be limited to the correction of a sort of “local bias” in the interpolated EMIC result. The potential to achieve this probably depends on the “diversity” in the calibration data: in the extreme case which would only use one calibration time-slice, such a map could, depending on the available degrees of freedom, represent most of the data. We are evidently not in this case, as there are 2 glacial time-slices + the recent past: the risk of exaggerating the contribution of this time-independent term (which might not be time-independent in the real world) is limited. However, this suggests that the model needs to be carefully validated, because there is a potential to obtain a very good fit for the calibration data, but much less good results when predicting other time periods even in the “calibration range” of each predictor. A related difficulty is that the model includes quite complex functions, in particular s7(x, y, TCLI). This particular contribution (s7) suggests that some locations are more “sensitive” to the EMIC’s output than others (could an interpretation be provided?). This is perhaps true, but again, there is a potential for this contribution to appear independent of time for the calibration cases, while it is not clear that it is more generally true (in the real world).


Introduction
Climate risk assessments and bioclimatic studies on a millennial timescale require regional climate data.Present state-ofthe-art comprehensive atmospheric general circulation models (GCMs), coupled with modules simulating the biosphere and sea ice, are major tools for the study of past, present and future climates.The resolution of such global models is usually 100-300 km (Flato et al., 2013).The relatively high resolution and complex calculations make GCMs suitable for modelling the climate on a 100-1000 yr timescale.For routinely simulating time periods covering full glacial cycles, the GCMs are computationally far too demanding.Nevertheless, several GCMs have been used for time-slice simulations of the distant past climate with prescribed ice sheets (Renssen et al., 2005;Otto-Bliesner et al., 2006;Brandefelt and Otto-Bliesner, 2009;Liu et al., 2009;Kjellström et al., 2010;Singarayer and Valdes, 2010;Strandberg et al., 2011;Braconnot et al., 2011Braconnot et al., , 2012)).Smith and Gregory (2012) used an AOGCM (atmosphere-ocean general circulation model) in a transient simulation of the last glacial cycle; however, in their work the ice sheet was prescribed and not coupled with the climate component of the model.
Earth system models of intermediate complexity (EMIC; Claussen et al., 2002;Petoukhov et al., 2005) provide a practical approach for simulating climate in timescales covering entire glacial cycles (∼ 100 000 yr).A reasonable computing time is achieved by reducing the problem's size (i.e.number of grid points) and model complexity.Full earth system models employ eddy-resolving general circulation models for atmosphere and ocean and are coupled with detailed submodels for bio-geochemical processes and cycles.EMICs, in contrast, are coarsely discretized and an atmosphere model is typically dynamical-statistical where mean flow is explicit and eddies are parameterized.Ocean models can be multibasin but zonally averaged.Despite these simplifications, in an EMIC intercomparison of eight different models (Petoukhov et al., 2005), the equilibrium and transient responses to a doubling of the atmospheric CO 2 concentration were within the range of corresponding GCM simulations.EMICs have been extensively utilized in simulating past glaciations (Berger et al., 1999(Berger et al., , 2003;;Wang and Mysak, 2002;Calov et al., 2005a, b) and future climate on a timescale of 100 000 yr (Berger and Loutre, 2002;Archer and Ganopolski, 2005;Cochelin et al., 2006).
The EMIC simulations provide continuous global data; however, there are uncertainties related to the model simplifications and low spatial resolution, and these make the direct output of EMICs unsuitable for regional studies.In climate science, several downscaling methods have been used.A common method is to downscale the results of a global model dynamically with a regional model.Regional models usually have sophisticated atmosphere and biosphere modules at a resolution of ∼ 50 km and are therefore computationally demanding.Hence downscaling of this kind is possible for short time-slice simulations of 50-100 yr, not for an entire 100 000 yr long simulation.A computationally less demanding method is statistical downscaling, whereby the statistical relationships between the observed small-scale variables (derived from observations) and larger-scale variables (e.g. from a global model) are derived using, for example, regression analysis.These derived statistical relations are then applied to downscale the large-scale variables of climate simulations to a smaller scale.Dynamical vs. statistical downscaling from hemispheric to regional scales involves a number of fundamental choices.It seems that the single most important factor affecting the quality of the downscaling is the simulation skill of the driving model (for a thorough discussion, see Racherla et al. (2012).Vrac et al. (2007) introduced a statistical downscaling method for palaeoclimatological purposes, based on generalized additive models (GAMs; Wood, 2006).They fitted a GAM-type regression model by finding the statistical relationships between the observed recent past climate  and the low-resolution CLIMBER-2 (CLIMate and BiosphERe Model 2) EMIC (Petoukhov et al., 2000(Petoukhov et al., , 2005) ) model simulation of recent past climate.They then used this regression model to downscale the recent past and last glacial maximum climates simulated by CLIMBER-2 over Western Europe.All these GAMs were calibrated with the recent past climate, and were used relying on the assumption that the statistical relations between the large and small scales remain unchanged for the last glacial maximum and the recent past climate.
We investigated the usability of GAMs in downscaling the large-scale variables of the CLIMBER-2 EMIC model simulations by Ganopolski et al. (2010) over Europe in climatic conditions ranging from glacial to interglacial.For this purpose we fit a GAM not only to downscale the observed recent past climate but also to downscale relatively high-resolution simulated glacial climate.For a glacial period characterized by an extensive ice sheet over Fennoscandia, about 21 kyr before present (BP), we utilized data from a CCSM4 (Community Climate System Model) GCM simulation by Brady et al. (2013) and an RCA3 (Rossby Centre regional Atmospheric climate model) regional model simulation by Strandberg et al. (2011) forced by a CCSM3 simulation by Brandefelt and Otto-Bliesner (2009).For a glacial period characterized by a small ice sheet over Fennoscandia, about 44 kyr BP, we utilized data from an RCA3 simulation by Kjellström et al. (2010).RCA3 as a regional climate model has higher spatial resolution than the global model CCSM4 and in this study RCA3 is used because of its presumably more detailed depiction of climate, whereas CCSM4 is used because its output enables downscaling to be done for any location on earth.As stated by Vrac et al. (2007), the use of statistical downscaling in palaeoclimatology is based on the assumption that the fitted regression model remains unchanged over time.We evaluated the GAMs by presenting their skills and investigate thereby which explanatory variables produce the best fit.We also tested the GAMs for downscaling the temperature of the CLIMBER-2 simulation from 10 kyr BP up to the present.The bilinearly interpolated CLIMBER-2 surface temperature and the downscaled temperature were compared to two temperature reconstructions: Laihalampi, Finland (Heikkilä and Seppä, 2003), and Gilltjärnen, Sweden (Antonsson et al., 2006).

Downscaling with generalized additive models
Statistical regression models are calibrated by establishing statistical relationships between large-scale variables (called explanatory variables or predictors, X 1 , . .., X p ) and a local variable (Y , called the response or predictand).Our largescale data are the global CLIMBER-2 output (Ganopolski et al., 2010) covering the last glacial cycle, with climate conditions ranging from glacial to interglacial.The range of the large-scale variables, X, defines the calibration range.Hence, the calibration range of the regression model should cover glacial to interglacial climates.For local data representing interglacial climate we utilize observations of the recent past climate by the Climate Research Unit, CRU.For local data representing climate for certain time steps of the last glacial cycle we utilize simulations by a GCM (CCSM4; Gent et al., 2011) and a regional model (RCA3; Samuelsson et al., 2011).These two data sets were used separately.The CCSM4 LGM (last glacial maximum) simulation (Brady et al., 2013) was used for the annual GAMs over the area we call Western Eurasia (36-70 • N, 10 • W-69 • E).The RCA3 downscaling (Strandberg et al., 2011) of the CCSM3 LGM simulation (Brandefelt and Otto-Bliesner, 2009) was used in the monthly GAMs over Northern Europe (54-70 • N, 3-35 • E).We aim at finding the relationships between the highresolution data and the large-scale data for the corresponding time periods.These relationships are then used with the large-scale data to obtain high-resolution data where they are otherwise not available.

Calibrating the statistical model
We used GAMs (Hastie and Tibshirani, 1990;Wood, 2006) in which the statistical expectation (E) of the predictand (Y ) were modelled using the sum of univariate spline functions of the p predictors (X 1 , . .., X p ), such that where the potentially non-linear spline functions f j had a non-parametric form.The calibration of GAMs is datadriven, and the shape of the response is not forced to any parametric form.GAMs are called semiparametric models, as the probability distribution of the predictand should be known.Our predictand was either the mean temperature (temperature 2 m above the surface) or the mean total precipitation.The variables used for the calibration of the monthly GAMs were monthly means of each grid point.Respectively, the annual GAMs were calibrated using annual means of each grid point.As temperature data classically satisfy the normality assumption, the response variable mean temperature was expected to follow a normal distribution.According to Cheng and Qi (2002), the cumulative precipitation data can be modelled by a lognormal distribution, hence the total precipitation response variable was log-transformed and expected to follow a normal distribution.The statistical relationships were calibrated using the data for 1 month (or annual) at a time.For fitting the GAM, all the data were to be represented at the same spatial resolution.We used a resolution of 1.5 • × 0.75 • over Northern Europe and Western Eurasia.These two regions are defined based on the available climate simulation data sets; the RCA3 simulations cover Northern Europe and the CCSM4 the area we call here Western Eurasia.As the size of the grid cells of our data varied with latitude, we associated a weight with each grid cell depending on the latitude of the data.The functions f j in Eq. (1) were defined by cubic regression splines with a restricted maximum number of degrees of freedom under wiggliness penalization (Wood, 2006).This means that although the maximum number of the degrees of freedom per predictor was restricted here to 15, the penalization allowed reducing the effective number of degrees of freedom of each spline to what was really needed by the problem.
As predictors on the right-hand side of Eq. ( 1) we used both meteorological and geographic predictors.As meteorological predictors we tested CLIMBER-2 large-scale temperature, relative humidity, precipitation, and lapse-rate data.As geographic predictors we tested terrain properties (elevation, shortest distance to the ocean, distance to nearest water body, slope angle, direction of the slope with respect to north), latitude, longitude, and shortest distance to the ice sheet margin.We tested altogether about 50 combinations of predictors and found the present to be the best.All the predictors were physically reasonable.Ice sheets have a large impact on local climate and that is why one of the predictors tested was selected to be the shortest distance to the ice sheet margin.The predictor side's temperature, relative humidity, precipitation, lapse rate and, over ice-free areas, the terrain properties were extracted from the CLIMBER-2-SICOPOLIS (SImulation COde for POLythermal Ice Sheets) simulation data by Ganopolski et al. (2010) at 0, 44 and 21 kyr BP to represent the recent past, a stadial during MIS (marine isotope stage) 3 and LGM, respectively.Over icecovered areas the terrain properties used as predictors were extracted from the data of the RCA3 and CCSM4 simulations.In earlier work GAMs have been developed using observed data only.However, if we want to have GAMs valid also for, for example, glacial periods, we need data depicting those conditions.The use of climate simulations as this work widens the data window and should hence lead to an outcome that covers various climate conditions better than the system that is developed using present data only.
The shortest distance to the nearest glacier margin used as predictor was computed from the RCA3 and CCSM4 ice sheet data.In ice-free areas the distance was expressed by positive values and in areas covered by ice the distance was expressed by negative values.In the calibration of the GAMs we implement the ice sheet extent and topography corresponding to the calibration data; e.g. when calibrating by RCA3 data we implement the RCA3 ice sheet extent and topography.Further, for predictions, the ice sheet extent and topography of SICOPOLIS is used.If the GAMs were calibrated by SICOPOLIS ice sheet data, we would have problems near the ice sheet margins as they differ in the SICOPO-LIS simulation and RCA3 and CCSM4.
The final predictors were determined by minimization of the generalized cross-validation (GCV) score, which was computed as where N stands for the number of observations, Y m for the GAM-modelled surface air temperature or log-precipitation (log meaning natural logarithm), and d stands for the effective number of degrees of freedom.The GCV score is an estimate of the expected mean square predictor error inside the calibration predictor values' range.The GCV has the same units as Y .Further, in computing the GCV, the effective degrees of freedom were also penalized by increasing their cost with a weight γ = 1.4 to prevent overfitting of the model (Kim and Gu, 2004).We allowed ourselves to combine several predictors, if that was found to improve the skill and was physically reasonable.We tested numerous combinations of predictors for fitting the GAMs in Eq. ( 1).Based on the GCV values, the best fits for downscaling the logprecipitation with GAMs were attained with the equation log(P ) (P CLI , x, y, h, d) = s 1 (P CLI ) + s 2 (x, y) (3) in which P represents the downscaled precipitation, P CLI is the precipitation of the large-scale data (CLIMBER-2), x and y are the longitude and latitude, respectively, h is the elevation, d is the direction of the steepest slope, and s 1 . . .s 4 are spline functions.
For downscaling the annual temperature and the monthly temperature of May, June, July, August and September with GAMs, the best fits were attained with the equation and for the monthly temperature of January, February, March, April, October, November and December with the equation + s 8 (h) +s 9 (d, s) , in which T is the downscaled temperature, T CLI is the temperature of the large-scale data (CLIMBER-2), the slope, s, is the angle of the steepest slope, i is the shortest distance to the ice sheet margin (positive if ice-free, negative if covered with ice) and s 5 . . .s 9 are spline functions.

Large-scale input data, CLIMBER-2-SICOPOLIS
The large-scale data to be downscaled were the output of a full last glacial cycle simulation (126 000 yr BP-present day) with the CLIMBER-2-SICOPOLIS model system by Ganopolski et al. (2010).CLIMBER-2 is an EMIC consisting of six earth system components: atmosphere, ocean, sea ice, land surface, terrestrial vegetation and ice sheets (Petoukhov et al., 2000).The model has a relatively low spatial resolution: for the atmospheric module the latitudinal resolution is 10 • and the longitudinal resolution is roughly 51 • .
The CLIMBER-2 model successfully describes the seasonal variability of a large set of characteristics of the climate system (Petoukhov et al. 2000), and simulates the climate response to changes in different types of forcing and boundary conditions within the range of corresponding GCM simulations (Ganopolski et al., 2001;Petoukhov et al., 2005).For simulating glacial climates, the CLIMBER-2 model has been coupled with the high-resolution three-dimensional thermomechanical ice sheet model SICOPOLIS (Greve, 1997).The resolution of the ice sheet model is 1.5 • × 0.75 • and the model domain extends in the Northern Hemisphere from 21 to 85.5 • N. The climate and ice sheet components are coupled bidirectionally using a physically based energy and mass balance interface (SEMI) as described in detail by Calov et al. (2005a).The simulations by Ganopolski et al. (2010) were forced by variations in the earth's orbital parameters calculated following Berger (1978), and atmospheric greenhouse gas concentrations derived from the Vostok ice core.The evolution of the radiative forcing of atmospheric dust was parameterized to alter proportionally with the simulated global ice volume as by Schneider et al. (2006).

Recent past climate
For the recent past climate predictands, we used the CRU high-resolution climate data, version 2.1.(Mitchell and Jones, 2005).For land areas the original resolution of the monthly mean surface temperature and total precipitation data was 0.5 • × 0.5 • .In practise this CRU data is the same as the often-used, reanalysed ERA-40 (ECMWF 40 Years Reanalysis) data set (Uppala et al., 2005).However, the use of ERA 40 or the other reanalysed data sets such as NCEP (National Centers for Environmental Prediction; Kistler et al., 2001) could be considered in the possible future research.To cover areas where CRU high-resolution data were not available, i.e. sea areas, we used the Jones et al. (1999) 5 • × 5 • temperature data.Similarly, where high-resolution precipitation data were not available, we used the Climate Prediction Center (CPC) Merged Analysis of Precipitation (CMAP; Xie and Arkin, 1997) derived from the years 1979-2000, with the original data resolution of 2.5 • × 2.5 • .

Glacial climate over Northern Europe
For the glacial climate predictands over Northern Europe we used the 50 yr of monthly mean temperature and precipitation values of two simulations with the regional climate model RCA3 (Kjellström et al., 2005, Samuelsson et al., 2011) forced by simulations with the GCM CCSM3 (Collins et al., 2006;Kiehl et al., 2006).The first data set was that of an RCA3 simulation by Kjellström et al. (2010), representing a stadial within the MIS 3 around 44,000 yr BP, a cold period with a relatively small Fennoscandian ice sheet as illustrated in  conditions.In particular, the CCSM3 model compares generally well in the PMIP3 (Kageyama et al., 2006) and it is within the ensemble spanned by participating models; for instance, regarding the simulated temperature of the warmest and coldest month in the latitude-longitude sectorial averages of interest.We thus conclude that, for the present study, the choice of CCSM3 is equally well justified as the selection of any other model in the PMIP ensemble.

Glacial climate over Western Eurasia
The glacial climate predictands over Western Eurasia were produced using annual mean temperature and precipitation data from two simulations with GCM CCSM4: the last 20 yr of the CCSM4 LGM simulation (Brady et al., 2013) and the period AD 1961-1990 of the CCSM4 recent past simulation (Gent et al., 2011).First we computed the differences between the CCSM4 simulated LGM and recent past climates and then we added these changes to the CRU's annual temperature and precipitation data sets (Sect.2.3.1).For temperature we applied the absolute change ( • C) and for precipitation the relative change (%).

Evaluation of the GAMs
For evaluating the fitted GAMs, we computed several statistical quantities as follows: (i) the percentage of explained deviance: where Y 0 stands for the predictand data (CRU/RCA3/CCSM4) and Y m for the GAM-modelled surface air temperature or log-precipitation; (ii) the spatial correlation between the CRU/RCA3/CCSM4 and the GAM-modelled sur-face air temperature or log-precipitation: For error estimates we computed the root-mean-square error, which is sensitive to the maximum errors, and the mean absolute difference, which is less sensitive to the maximum errors.
Further, in Supplement 1, we tested the monthly GAMs in predicting other months' temperatures or precipitation over Northern Europe.In Supplement 2, the annual GAMs for Western Eurasia (calibrated only by recent past and LGM data) were used for predicting annual mean temperature and precipitation for 44 kyr BP, and the results were compared to the output of the RCA3 regional climate model simulation by Kjellström et al. (2010).

Comparison with temperature reconstructions
The calibrated GAMs were applied to downscale the Holocene climate (10 kyr BP-present) of the simulation by Ganopolski et al. (2010).The resulting annual mean temperature was compared with two pollen-based reconstructions of temperature for Laihalampi, Finland (Heikkilä and Seppä, 2003), and Gilltjärnen, Sweden (Antonsson et al., 2006).The RMSE and MAD were computed for the bilinearly interpolated CLIMBER-2 surface temperature and reconstructions as well as for the downscaled temperature and reconstructions.

Evaluation of GAMs for precipitation
Equation (3) resulted in the skills presented in Table 1.The precipitation GAMs showed good correlation (> 0.72) with the fitted data.The MAD (RMSE) of the monthly precipitation GAMs varied between 7 mm month −1 (12 mm month −1 ) in February and 12 mm month −1 (21 mm month −1 ) in October.The GAMs were able to explain the spatial variance of the total precipitation by 59-85 %.
Figure 2 depicts the splines of the GAM downscaling January log-precipitation over Northern Europe.Figure 2a www.clim-past.net/10/1489/2014/Clim.Past, 10, 1489-1500, 2014  demonstrates the effect of the raw CLIMBER-2 January precipitation: the precipitation increased with CLIMBER-2's precipitation, which indicates that the spline function reduces the precipitation in months of low precipitation simulated with CLIMBER, and enhances the precipitation in months of high CLIMBER-simulated precipitation.This further leads to the conclusion that the spline function increases the amplitude of the precipitation annual cycle, which is a well justified assumption because of CLIMBER's coarse resolution and very small precipitation variability.Figure 2b shows the effect of the longitude and latitude: the spline increased precipitation on the west coast of Norway, as well as in Denmark and southern Sweden and decreased it over central and northern Finland.The effect of elevation on precipitation depended on the season.In the annual GAM and March-September GAMs, the precipitation increased with elevation (not shown).However, for cooler months, October-January, the elevation effect was weaker, and for elevations between 1000 and 2000 m the precipitation even decreased with elevation (Fig. 2c).To speculate, this may have a meteorological explanation: at relatively low elevations (< 750 m) topography enhances precipitation, between 750 and 2000 m elevations on glaciers the top climate is characterized by smaller precipitation amounts and, again, when we go above 2500 m (mountain tops) precipitation increases as a function of elevation.The effect of the direction of the slope is shown in Fig. 2d: southerly and westerly slopes (∼ 180 to ∼ 270 • ) tended to increase precipitation, whereas northerly and easterly slopes (∼ 0 to ∼ 90 • ) decreased it.Table S1.1 in Supplement 1 shows the skills of the monthly GAMs in predicting other months' precipitation over North-Table 1. Skills of the monthly total log-precipitation GAMs for the Northern European area, and of the annual total log-precipitation ( * ) GAM for the Western Eurasian area.Values in parentheses represent total precipitation.The GCV score is defined in Eq. ( 2).The percentage of explained deviance (%ED), the spatial correlation (Cor), the root-mean-square error (RMSE), and mean absolute difference (MAD) are defined in Eqs.(6-10), respectively.

Month
GCV ern Europe; the correlations and errors are in the same range as for the fitting data, suggesting that these GAMs do make reasonable predictions with other than just calibration data.

Comparison of downscaled precipitation with observation
Figure 3 demonstrates the raw CLIMBER-2 January precipitation (top row in Fig. 3), the CRU (first panel of the second row in Fig. 3) and RCA3 (second and third panels of the second row in Fig. 3) regional-scale precipitation, the GAM downscaled January precipitation (third row in Fig. 3), and the difference between the CRU/RCA3 data and the GAM downscaled precipitation (bottom row in Fig. 3).By comparing rows three and two in Fig. 3, we see that the GAM reproduced well the spatial distribution of precipitation over Northern Europe, the highest precipitation rates occurring on the western slopes of the Scandinavian mountains and the smallest over the central parts of the Fennoscandian ice sheet.

Evaluation of GAMs for temperature
Equations ( 4) and ( 5) resulted in the skills presented in Table 2.For the temperature GAMs the correlations were high, explaining 96-99 % of the spatial variance of mean temperature.The MAD (RMSE) of the temperature GAMs varied between 0.6 • C (0.8 • C) in April and 1.6   calibration ranged from a glacial to an interglacial climate.The CLIMBER-2 temperature was used directly, i.e. linearly (see Eqs. 4,5).With the splines for longitude and latitude some of the cold bias of CLIMBER-2 over Northern European land areas was corrected, as seen for the July temper-ature GAM in Fig. 4a.Over Western Eurasia, the longitude and latitude splines produced warming over Western Europe and cooling over continental Russia (not shown).In all the GAMs, the splines decreased the temperature with increasing altitude, as shown for July in Fig. 4b.In the May-September and annual GAMs (Eq.4) the splines also decreased the temperature with increasing distance to ice-free areas, as shown for July in Fig. 4b.In the GAMs for October-April (Eq.5) the direction of the slope also seemed to have an influence on the local temperature: the splines increased the temperature on the westerly slopes and decreased it on the easterly slopes (not shown).Table S1.2 in the Supplement 1 shows the skills of the monthly GAMs in predicting other months' temperatures over Northern Europe.For the months May-September the correlations and errors are in the same range as for the fitting data (Table 2), suggesting that these GAMs do make reasonable predictions with other than just calibration data.For the temperature GAMs of the months October-April, the monthto-month test gave higher errors than for the calibration data (Table 2), suggesting that these GAMs are not suitable for downscaling other months than the fitted one.

Comparison of downscaled temperatures with observations
In Fig. 5 the bilinearly interpolated CLIMBER-2 July temperatures are in the top row, the CRU-observed July temperature is in the first panel of the second row, the RCA3simulated July temperatures are in the second and third panels of the second row, the GAM-downscaled July temperatures are in the third row, and the differences between the CRU/RCA3 data and the GAM-downscaled temperatures are in the bottom row.By comparing the second and third rows in Fig. 5, we see that the GAM reproduced well the spatial distribution of July temperature over Northern Europe, e.g. the Scandinavian mountains and the strong temperature gradient near the ice sheet margin were well brought out by the GAM.The GAM seemed to work especially well on ice-free land areas that were not in the close vicinity of ice sheets or mountains, the residuals being typically 0-2 • C.There were, however, also some differences, e.g. the downscaling of the recent past climate over the Scandinavian mountains was somewhat too cold, and the temperature gradient on the eastern margin of the ice sheet in the downscaled GAMs was less steep than in the RCA3 simulation.
Figure 6 shows pollen-based reconstructions of annual mean temperature, the bilinear interpolation of the annual mean surface temperature simulated by CLIMBER-2, and the GAM-downscaled annual mean temperature at two locations in Northern Europe during the Holocene.The pollenbased reconstructions were produced with a weighted averaging partial least squares regression and calibration technique.Bootstrapping-based root-mean-square errors of prediction for the reconstructed values are sample-specific and vary from 1.0 to 1.5 • C (Birks et al., 2010).When proxybased reconstructions are used for comparison with model results, it must be borne in mind that the output of such re-constructions have been shown to be sensitive to the spatial extent of the calibration data sets used (Salonen et al., 2013).However, while the absolute values are highly sensitive to the climatic characteristics of the calibration data set, the shapes of the relative palaeotemperature curves seem comparatively robust, as the curve shapes mostly remain similar as the calibration data is spatially shifted (Salonen et al., 2013).Hence, it is important to note that the shapes of the reconstructed temperature curves are generally consistent with those based on GAM downscaling for the Holocene (Fig. 6).As for the absolute values, the MAD of the reconstructions and the bilinear interpolation of the CLIMBER-2 surface temperature was larger than 3.8 • C and the RMSE > 3.9 • C. The MAD of the reconstructions of the GAM output was less than 1.4 • C and the RMSE < 1.6 • C. The difference between GAM_Western Eurasia and GAM_Northern Europe is due to the difference of the driving model; that is, CCSM4 versus RCA3.
Based on Fig. 6, it seems at first that the GAM is simply a bias correction of CLIMBER simulation towards the palaeoreconstructions.The bias of the reconstructions is however unknown and thus we cannot definitely say that the GAM is an improvement over CLIMBER only because it is closer to the reconstructions (we of course hope it is, but we cannot say).Bearing in mind that the reconstructions are independent of the GAM and CLIMBER and the climate model data from CCSM4 and RCA3 are only one depiction of climate, it is nevertheless satisfying to see that the GAM is in proximity of observed climate.If we believe in the reconstructions, then we can conclude that the additional information provided to the GAM, on top of the CLIMBER simulation, is adequate and very general since the locations of the reconstructions are completely random for the GAM, and we have all reasons to further believe that the GAM would fit reconstructions in   Heikkilä and Seppä (2003) for Laihalampi, and by Antonsson et al. (2006) for Gilltjärnen.The green contours are interpolations from CLIMBER-2 simulations by Ganopolski et al. (2010).The GAM_Western Eurasia (red curves) and GAM_Northern Europe (purple curves) data were downscaled by GAM models from the CLIMBER-2 data in this paper.
other locations as well.Unfortunately, there are extremely few sites with such data available.

Discussion and conclusions
We found that with an apt selection of predictor parameters, GAMs can be useful to statistically downscale low-resolution CLIMBER-2 EMIC simulation data ranging from a glacial to an interglacial climate.For calibrating the GAMs, we utilized both physical predictors -such as the large-scale precipitation and temperature data of CLIMBER-2 -and geographical predictors -such as elevation, latitude, longitude, the direction of the steepest slope, the angle of the steepest slope, and the shortest distance to the ice sheet margin.The final selection of predictors was based on the statistical skill of each GAM.
The GAMs fitted for temperature (precipitation) were able to approximately reproduce the observed and modelled temperature (precipitation) fields.The mean errors in the mean temperature (total precipitation) were of an order of magnitude of 1-2 • C (10-20 mm month −1 ) as the climate varied from glacial to interglacial over Europe.These error estimates are at the same scale as in the GAMs by Vrac et al. (2007).However, the calibration ranges in Vrac et al. (2007) covered only an interglacial climate, whereas our fitting included a glacial climate as well.Though some of the detailed spatial features were not fully captured in regions with large spatial variation of topography, the main spatial patterns were relatively well captured by the GAMs taking into account the very coarse resolution of downscaled parameters.
In Supplement 1 the monthly GAMs were used to predict other months.The skills of precipitation GAMs and temperature GAMs for months May-September in predicting other months were in the same range as for the calibration months.For the October-April temperature GAMs the errors were high, suggesting that these GAMs are not capable of downscaling outside of their calibration range.Also, Vrac et al. (2007) found that some of their GAMs calibrated over Western Europe were not valid outside their calibration domain (over North America or Northern Europe).
One critical assumption in the use of the statistical downscaling method is that the statistical relationship of the climate parameters used should remain stationary over time (Vrac et al., 2007).This assumption can be questioned with palaeoclimatic simulations spanning tens of thousands of years and reaching back to the last glacial with markedly different climatic boundary values than at present.For the validation of the statistical downscaling method, we compared our simulated Holocene temperatures with two pollenbased quantitative annual mean temperature reconstructions from Finland and Sweden.The results of these comparisons were generally congruent and thus support the validity of the statistical downscaling approach, at least in Northern Europe.Moreover, GAM-downscaled CLIMBER-2 temperature showed better agreement with the pollen-based reconstructions than the bilinearly interpolated CLIMBER-2 surface temperature.
The downscaling method enables the generation of highresolution climate data from the low-resolution CLIMBER-2-SICOPOLIS simulations of the past.These data could be used in bioclimatic modelling and long-term climate assessments.The first comparison of the GAM-downscaled CLIMBER-2 temperature to palaeoclimatological recon-structions showed better agreement than the bilinearly interpolated CLIMBER-2 temperature.To further study the correspondence between these downscaled simulations, the next step will be to downscale the precipitation and temperature with the fitted GAMs over the whole last glacial cycle simulation, and compare the results with palaeoclimatological reconstructions.Finally, the GAMs could be further developed by adding predictors or additional simulation data such as the ensemble mean of the PMIP3 LGM simulations.
The Supplement related to this article is available online at doi:10.5194/cp-10-1489-2014-supplement.

Figure 2 .
Figure 2. Components of the GAM model in Eq. (3) fitted to downscale the CLIMBER-2 January log-precipitation in the recent past, MIS 3 stadial and LGM climates.The blue curves show the estimated effects of (a) CLIMBER-2 total precipitation, (b) longitude and latitude, (c) elevation, and (d) direction of the steepest slope.

Figure 3 .
Figure 3. 782 783 784 Figure 3. January mean total precipitation at 0 (left column), 44 (middle column) and 21 kyr BP (right column) as simulated by the global model CLIMBER-2 (top row), as observed during the period 1961-1990 by the CRU (first figure of the second row), as simulated by the regional model RCA3 (second and third figures of the second row), as predicted by the GAM model (third row), and as the difference between the statistical model and observation/simulation (bottom row).The unit is millimetres.The data of the first and the second rows have been bilinearly interpolated to a 1.5 • × 0.75 • resolution.31

Figure 4 .
Figure 4. Components of the GAM model in Eq. (4) fitted to downscale the CLIMBER-2 July mean temperature in the recent past, MIS 3 stadial and LGM climates.The blue curves show the estimated effects of (a) longitude and latitude, and (b) elevation and distance to the nearest ice sheet (or ice cap) margin.In ice-free areas the distance was expressed by positive values and in areas covered by ice the distance was expressed by negative values.

Figure 5 .
Figure 5. 790 791 Figure 5. July mean temperature at 0 (left column), 44 (middle column) and 21 kyr BP (right column) as simulated by the global model CLIMBER-2 (top row), as observed during the period 1961-1990 by the CRU (first panel of the second row), as simulated by the regional model RCA3 (second and third panels of the second row), as predicted by the GAM model (third row), and as the difference between the statistical model and observation/simulation (bottom row).The unit is degrees Celsius.The data of the first and the second rows have been bilinearly interpolated to a 1.5 • × 0.75 • resolution.33

Figure 6 .
Figure 6.Comparison of simulation data and pollen-based reconstructions of annual mean temperatures from (a) Laihalampi in Finland and (b) Gilltjärnen in Sweden.The blue curves represent reconstructions byHeikkilä and Seppä (2003) for Laihalampi, and byAntonsson et al. (2006) for Gilltjärnen.The green contours are interpolations from CLIMBER-2 simulations byGanopolski et al. (2010).The GAM_Western Eurasia (red curves) and GAM_Northern Europe (purple curves) data were downscaled by GAM models from the CLIMBER-2 data in this paper.

Table 2 .
Skills and calibration ranges of the monthly mean temperature GAMs for the Northern European area, and of the annual mean temperature ( * ) GAM for the Western Eurasian area.Same abbreviations as in Table1.