CHELSA-TraCE21k v1.0. Downscaled transient temperature and precipitation data since the last glacial maximum

High resolution, downscaled climate model data are used in a wide variety of applications across environmental sciences. Here we introduce the CHELSA-TraCE21k downscaling algorithm to create global monthly climatologies for temperature and precipitation at 30 arcsec spatial resolution in 100-year time steps for the last 21,000 years. Paleo orography 10 at high spatial resolution and at each timestep is created by combining high resolution information on glacial cover from current and Last Glacial Maximum (LGM) glacier databases with the interpolation of a dynamic ice sheet model (ICE6G) and a coupling to mean annual temperatures from CCSM3-TraCE21k. Based on the reconstructed paleo orography, mean annual temperature and precipitation was downscaled using the CHELSA V1.2 algorithm. The data were validated by comparisons with the glacial extent of the Laurentide ice shield based on expert delineations, proxy data from Greenland ice cores, historical 15 climate data from meteorological stations, and a dynamic simulation of species distribution throughout the Holocene. Validations show that CHELSA TraCE21k output creates a reasonable representation of the distribution of temperature and precipitation through time at a high spatial resolution, and simulations based on the data are capable of detecting effective LGM refugia of species.


Introduction 20
Since the Last Glacial Maximum (LGM), variation in climate has caused multiple changes of the Earth surface, including the rearrangement of species distributions, or even extinctions (Adams and Faure, 1997;Binney et al., 2017;Prentice et al., 1991;Velichko et al., 1997;Williams et al., 2004;Yu et al., 2010). Yet we have not fully evaluated the historical underpinnings of these changes as we have often lacked the climate data at the necessary spatial resolution. Biological entities such as species usually encounter climatic conditions at spatial resolutions (<1km 2 ) (Seo et al., 2009) that are not captured by numerical global 25 circulation models (GCMs) which run at much coarser grains (e.g. >0.5°). For many applications such as inference of ecological niches (Hutchinson, 1957), determination of growing seasons (McMaster and Wilhelm, 1997), identification of species migrations (Engler and Guisan, 2009), or modelling of high resolution species distributions (Guisan and Thuiller, 2005;Guisan and Zimmermann, 2000), temporal and spatial variability in temperature and precipitation is highly important. For such analyses, errors in the underlying climate data can strongly deteriorate the analytical power (Soria-Auza et al., 2010). For the recent past the gap between the coarse GCM resolution and the high resolution needed has be bridged using satellite data (Biasutti et al., 2011;Funk et al., 2015;Huffman et al., 2007), statistical downscaling (Karger et al., 2017aMaraun et al., 2010;Schmidli et al., 2006;Wilby et al., 1998;Wood et al., 2004), dynamical downscaling (Skamarock et al., 2019), or interpolation of meteorological station data (Daly et al., 1997;Harris et al., 2020;Hijmans et al., 2005;Meyer-Christoffer et al., 2015). While all of these methods work comparably well for current climatic conditions, satellite or station data are not 35 available before the 19 th Century (end of the 20 th century for satellite data), hampering an application of said methods to paleoclimatic models. Most paleoclimatic data at high spatial resolution is therefore based on climatologically aided interpolation of GCM output (Brown et al., 2018). This process uses the high-resolution information of current day climatologies, and adds an interpolated anomaly derived from a coarser-resolution GCM (Hunter and Meentemeyer, 2005;Willmott and Robeson, 1995). While this approach works rather well for short term time series where topography is relatively stable (Daly et al., 40 1997), it becomes difficult for longer time series where the dependence structure between variables (e.g. topography and climate) is dynamic (Maraun, 2013). This phenomenon is of concern especially in the last 21'000 years, as the topography in many regions on earth has changed drastically due to the retreating ice shields and glaciers along the poles and in high mountain areas (Scotese, 2001). While numerical climate models are able to simulate paleo-environmental conditions comparably well (Sepulchre et al., 2020), they are however, computationally very demanding, and therefore they have not been applied on 45 ecologically relevant spatial resolutions of ~1km yet. Current global kilometer-scale models only show a simulation throughput of 0.043 SYPD (simulated years per day) , which is 25-fold lower than computater-efficient simulations of 1 SYPD (Schär et al., 2019;Schulthess et al., 2018). Even with state-of-the-art supercomputers and climate models this gap can only be minimized by a factor of 20 (Neumann et al., 2019).
Climate impact studies however, often only use a reduced set of climate variables compared to those available from the 50 numerical output of climate models (Frieler et al., 2017). Such studies do not need therefore a complete representation of all climate processes at high spatial resolution. In ecological studies, for instance, precipitation is often used along with minimum and maximum temperatures for analyses of species occurrences (Woodward et al., 1990). Also, it is common practice to describe species ranges by their climate envelopes and building species distribution models (SDMs) with a relatively small set of climate predictors based on monthly minimum and maximum temperature and precipitation (Guisan and Thuiller, 2005;55 Guisan and Zimmermann, 2000).
Here we present paleo-climatic data, downscaled from the CCSM3_TraCE21k model to a 30 arcsec resolution using the CHELSA V1.2 algorithm (Karger et al., 2017a), which covers time steps of 100 years from 21k-BP to 1990, for minimum and maximum temperatures, surface precipitation, and paleo-orography. The CCSM3 TraCE-21k climate model (Carlson et al., 2012;He, 2011;Liu et al., 2009;Marcott et al., 2011) provides information on climate change over the last 21,000 years, i.e. from Last Glacial Maximum (LGM) to present. The TraCE-21k simulations are capable to reproduce many main features of post-glacial climate dynamics in various parts of the world form low to high latitudes and including abrupt climate changes (He, 2011;Liu et al., 2009). CCSM3 TraCE-21k uses the 65 Community Climate System Model version 3 (CCSM3)of the National Center for Atmospheric Research (NCAR), which is a global climate model with coupled ocean, atmosphere, sea ice, and land surface components (Collins et al., 2006). CCSM3 TraCE-21k was calculated at T31_gx3v5 resolution (Otto-Bliesner et al., 2006)

Transient glacial extent simulations: ICE6G_C
The ICE6G_C model is a refinement of the ICE5G (VM2) model (Peltier, 2004) which has been widely used to model the distribution and dynamics of major ice shields through time. The ICE6G_C explicitly models changes in ice thickness of major 80 ice sheets (e.g. the Laurentide ice sheet) from Last Glacial Maximum (LGM) till today (Argus et al., 2014;Peltier et al., 2015) at 500 year time steps.

Observational glacial extent at Last Glacial Maximum (LGM)
As the extent of the glaciers during LGM, we use data from Ehlers et al. (2011) that presents an up-to-date, detailed overview 85 of Quaternary glaciations all over the world, not only with regard to stratigraphy but also with regard to major glacial landforms and the extent of the respective ice sheets.

Observational current glacial extent: GLIMS
The GLIMS (Global Land Ice Measurements from Space) project (Raup et al., 2007) at the NSIDC (National Snow and Ice Data Center) provides data on global glacial extent and other information about glaciers including metadata on how those 90 outlines were derived. Here we use this database to delineate the current extent of the glaciers at high resolution globally.

Methods:
Downscaling is based on a derived version of the CHELSA V1.2 algorithm (Karger et al., 2017a) and forcings from CCSM3 TraCE-21k (He, 2011;Liu et al., 2009) and involved four main steps. (1) To apply the CHELSA algorithm on the TraCE-21k 95 data, we first approximated the paleo-orography using a combination of data from the digital elevation model GMTED2010 (Danielson and Gesch, 2011), large scale ice sheet configurations from ICE6G (Peltier et al., 2015), high resolution glacier extents from GLIMS (Raup et al., 2007) and (Ehlers et al., 2011), as well as sea level change data from (Miller et al., 2005).

Paleo-orography
The first step in the downscaling was to combine annual mean temperature data with sea surface elevation and ICE6G data to estimate combined ice + surface orography. The rationale behind this approach is that temperature and glacier extents are to some degree interdependent, and a change in temperature will translate to a change in glacial extent. This approach was implemented by first combining the high resolution glacial extents from (Ehlers et al., 2011) with the ICE6G data. To do so, 105 we randomly sampled 100 point location per 1° grid cell from ICE6G and extracted the height of the glacier plus the surface elevation at time 0 (LGM). All points that did not fall within the high-resolution glacial extent were omitted. We then converted the polygons of the glacial extents into point locations as well, and extracted their elevation from the DEM at time , where = − with being the difference between the current and past sea surface elevation at time t.
Then both point datasets were combined to containing a point sample of the surface orography. This point sample was 110 then spatially interpolated to a grid of 30 arcsec resolution applying multilevel B-spline interpolation. The multilevel B-splines use a B-spline approximation to and starts using the coarsest grid 0 from an overall set of grids 0 , 1 , … , with n = 14 generated using optimized B-spline refinements (Lee et al., 1997). The resulting B-spline function 0 ( ) then gives the first approximation of et. 0 ( ) and leaves a deviation: 115 at each location ( , , ). Then the next control lattice 1 is used to approximate 1 (∆ 1 ). This approximation is then repeated on the sum of: at each point ( , , ) times resulting, in our case, in the gap free interpolated glacial surface . The interpolated glacial surface was then combined with to the orography using:

Temperature coupling
The estimated orography at the LGM was then used to downscale mean annual temperature =0 . GCMs such as the CCSM3 normally exhibit a large bias in temperatures or precipitation (Cannon et al., 2015;Maraun, 2013). We therefore 130 applied a change factor bias correction on current mean annual temperatures from CHELSA V1.2 resampled to a 0.5° grid resolution, where the change factor was calculated from the ∆ = − . This effectively preserves the trends observed in temperature, but simultaneously assumes that the bias has been conserved over time (Maraun, 2016). For the interpolation of ∆ to the 0.5° grid resolution of the same multilevel B-spline interpolation was used as described above. The bias corrected temperatures are then given by = − ∆ . To achieve a high-resolution 135 temperature approximation, we used a lapse rate-based downscaling from atmospheric temperature data at the TraCE-21k pressure levels. The lapse rates are based on a linear approximation from average temperatures at altitudes at vertical levels 26 to 20 of the T31 grid so that: Temperature at the surface at a high spatial resolution was then calculated by: To couple the outline of the glacier to the mean annual temperature, we transformed the glacier elevations = < , to a polygon, and then transformed the outline of this polygon to a point sampling of the glacier boundaries at time t: .
Temperatures at the glacier boundary where then extracted for each boundary location coordinate ( , ) during the LGM ( = 0) which gives the local temperature under which a glacial boundary existed during the LGM ( 0 , 0 ). As this simple coupling between glacial boundaries and annual temperatures is however only a rough estimate, we additionally 150 included current temperatures and calculated the difference between current and LGM boundary temperatures. Here, the locations for the current glacial boundaries were extracted ( , ) and locations of current and LGM glacial boundary temperatures were combined. The resulting point locations were then interpolated using a multilevel-B-spline to result in a gap-free surface ∆ . As the orography for the next time step is not known yet, we estimated the near surface air temperature +1 for the glacial melt similar the orography from the time step before. So that: 155 The binary glacial extent +1 at + 1 is then approximated as: here the second term in the condition for +1 linearly scales ∆ over the entire number of time steps. This correction is necessary, as otherwise the entire bias would be added at the first time step, resulting in an unrealistic strong shift in the glacial extent. We then repeated the transformation of the glacial extent +1 to all point locations and repeated the procedure for the 165 temperature coupling to estimate the orography +1 . Near surface air temperatures for + 1 have been then approximated using:

Precipitation estimation
The estimation of high resolution precipitation follows a variant of the CHELSA V1.2 algorithm (Karger et al. 2017). The CHELSA V1.2 algorithm assumes elevation to be the main driver of vertical precipitation gradients, which is however, highly idiosyncratic (Basist et al., 1994;Böhner, 2006;Böhner and Antonic, 2009;Daly et al., 1997;Gao et al., 2006;Karger et al., 2017a;Sevruk, 1997;Spreen, 1947). In tropical convective regimes, precipitation typically increase up to the condensation 175 level around 1000-1500 m above surface, while exponentially decreasing moisture content in the mid-to upper troposphere results in a drying above the condensation level resulting in non-linear precipitation lapse rates (Körner, 2007). Furthermore, negative precipitation lapse rates are common under the extremely dry polar climates. In contrast, at mid-latitudes and in the subtropics, precipitation generally increases with increasing elevation due to advection. As a consequence, summits of the Alps or other high mountain ranges exhibit high rainfall (Rotunno and Houze, 2007), and lapse rates for precipitation are 180 almost linear (Weischet and Endlicher, 2008). To approximate the effects of orographic precipitation we used a variant of the CHELSA V1.2 algorithm, which is explained in detail below.
We used 10 m u-wind and v-wind of TraCE-21k to calculate wind direction. Both wind components were projected to a world Mercator projection and then interpolated to a 4 km grid resolution using a multilevel B-spline interpolation similar to the one used for the bias correction surface. Windward and leeward effects are assumed to be best represented at resolutions larger 185 than 1km (Daly et al., 1994), we therefore chose a grid resolution of 4km for the underlying digital elevation model. The wind effect H was then calculated using: as well as in the free atmosphere above the planetary boundary layer (Daly et al., 1997;Karger et al., 2020;Oke, 2002;Stull, 2012). The highest impact of orography is considered just at the boundary layer height where the airflow interacts with the terrain. We used the lifted condensation level (LCL) as indicator of the altitude at which the wind effect exerts the highest contribution to precipitation. The LCL has been calculated using the mean air temperature (tas) and mean relative humidity (RH) using: 205 (Lawrence, 2005). The LCL has been interpolated to the CHELSA resolution using a B-spline interpolation. To create a boundary layer height corrected wind effect The wind effect grid containing was then proportionally distributed to all 210 grid cells falling within a respective T31 grid cell using: With being the being the maximum distance between the lifted condensation level at elevation and all grid cells 215 at a 30 arc sec. resolution falling within a respective T31 grid cell, ℎ being a constant of 9000 m, and being the respective elevation from GMTED2010 (Danielson and Gesch, 2011) with: 220 being the height of the monthly means of daily mean boundary layer from TraCE21k, being the elevation of the TraCE21k grid cell, and being a constant of 500 m which takes into account that the level of highest precipitation is not necessarily at the lower bound of the LCL, but slightly higher (Karger et al., 2017a).
The wind effect algorithm cannot distinguish extremely isolated valleys inside highly-elevated mountain areas (Frei and Schär, 1998). Such valleys are situated in areas where the wet air masses flow over an orographic barrier and are prevented from 225 flowing into deep valleys. These effects are mainly confined to large mountain ranges, and are not as prominent in small-to intermediate-sized mountain ranges . Such dry valleys are situated in areas where the wet air masses flow over an orographic barrier and are prevented from flowing into deep valleys . To account for these effects, we used a variant of the windward leeward equations with a linear search distance of 300 km in steps of 5° from 0° to 355° circular for each grid cell. The calculated leeward index was then scaled towards higher elevations using: 230 was set to 9000 m.. ℎ has been set to 9000 m. * will give the first approximation of precipitation intensity .
To achieve the distribution of monthly precipitation given the approximated precipitation intensity at each grid location 235 ( , ), we used a linear relationship between and using: where equals the number of 30 arc sec. grid cells that fall within a 0.25 grid cell.

Output validation
Direct validation of the temperature (Fig 1) and precipitation (Fig. 2) output at high resolution for paleo time series generally relies on proxies, as direct observations of both variables are not available. Although global temperature time series exist, they only give global means, and do not allow validation of the performance of a 1km paleo climatic dataset. Therefore, to validate 245 the CHELSA_TraCE21k dataset we complement a simple comparison of the simulated time series to proxy data and current observations, with approaches of validating derived parameters from the simulated temperature and precipitation that directly benefit from a very high horizontal resolution.

Validation using current (historical) observations 265
We used data from the Global Historical Climate Network (GHCN) monthly database V.3 (Lawrimore et al., 2011) to validate the performance of the downscaling algorithm during the last time step of the CHELSA-TraCE21k for which station data is available. To do so we calculated monthly climatologies for each month for tasmax, tasmin, and pr from both TraCE21k and CHELSA-TraCE21k. We then compared the values measured at each station to that simulated in both, TraCE21k and CHELSA-TraCE21k. 270 Although the downscaling algorithm might increase the performance of the temperature and precipitation estimates during the historical period, this does not imply that this improvement is equal during the entire transient time series. To further validate the data, we therefore compared it to more derived parameters for which time series data exists.

Comparison with temperature proxies from Ice core data
We compared the downscaled temperatures with the Greenland ice core reconstructions of Buizert et al. (Buizert et al., 2014(Buizert et al., , 275 2018 to check the performance of the downscaling at eight ice core locations on the Greenland Ice Shield (GIS). Although both temperature reconstructions and GCM generated temperatures have uncertainties connected to them (Erb et al., 2018), the ice core data are so far the best possible validation dataset that spans the entire deglaciation period from 21kBP to 1990 (Buizert et al., 2014(Buizert et al., , 2018. To assess the performance gain of the downscaling over the coarse resolution TraCE21k data, we compare the ice core annual mean near surface temperature reconstructions with both the CHELSA-TraCE21k, and the original 280 TraCE21k temperature data.

Validation of glacier extent between 18kPB and 1kBP
As the ice core temperature reconstructions have associated uncertainties, it is impossible to disentangle if potential differences between the ice core data and the model data are due to uncertainties in the reconstructions. To validate the downscaled temperature data further, we used the interpolated extent of glaciers in CHELSA_TraCE21k, and compared it to glacial extent 285 data from Dyke et al. 2003. The data consists of expert delineated glacier maps based on a chronological database of radiocarbon dates and contains >4000 dates located in North America (Dyke et al., 2003). To compare both datasets, we first calculated the glacial extent from CHELSA_TraCE21k by assigning each 1km gridcell in a Lambert Conformal Conic projection (+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +x_0=0 +y_0=0 +ellps=clrk66 +units=m +no_defs, extent= xmin : -3698534, xmax: 3596392, ymin: -1259083, ymax : 4674296). to either to be covered by a glacier [1] or being free of glacier 290 [0]. We assigned a 1 if the simulated glacier height was above the paleo terrain elevation, and a 0 if it was lower or equal to the paleo terrain elevation. The paleo terrain elevation was calculated using the current terrain elevation minus the sea level difference between current and that of the respective paleo timestep. As the current terrain elevation already includes extent glaciers, this elevation dependent procedure of assigning glacial extents would result in the current glaciers being assigned a 0. Therefore we assigned all grid cells covered by extent glaciers a 1. To compare the simulated glacial extent to the expert delineated extent, we rasterized the polygons provided by Dyke et al., 2003 for the years 18 kBP -1 kBP to the 1km resolution, extent, and projection of the simulated glacial cover and assign a 1 where the polygon intersects with a 1km raster cell and a 0 otherwise.
We then calculated three different test values to identify if the simulations correctly predict the presence and absence of a glacier. As the dataset is highly unbalanced between absences of glaciers [0] and presences of glaciers [1] through time we use 300 balanced accuracy which is defined as: (sensitivity+specificity)/2. Additionally we report Cohen's Kappa, and the True Skill Statistic (Allouche et al., 2006).

Accuracy of temperature and precipitation estimated during the current (historical) period
The original TraCE21k data shows large deviations and root mean squared errors (RSME) from the observed data (Fig. 3). 305 This is somewhat expected as a climate model running for such long time periods needs to have coarse resolution, as well as a large degree of generalism and realism which decreases the accuracy of a model when compared to observations. The temperature variables perform well in TraCE21k with r ~ 0.8 for all months, but have deviations and RSME similar to those of precipitation, that most likely can be attributed to the coarse resolution of the climate model. Trace21k also seems to overestimate temperature extremes both for tasmax and tasmin (Fig. 4). 310 The precipitation however, does not perform well in the model with r ~ 0.4 and large deviations from actual values (Fig. 3) and overall precipitation seems to be too low in the model (Fig. 4).
The downscaling algorithm improves the correlation between observed and modelled and decreases the standard deviation for all three parameters (Fig. 3). The downscaling for the temperature variables increases the correlation to r ~ 0.95 for all months and decreases the standard deviation substantially (Fig. 3). Similarly the performance of the precipitation estimation in 315 CHELSA-TraCE21k increases, which is reflected in a r ~0.7 and a lower standard deviation and RSME (Fig. 3). The underestimation of precipitation in the TraCE21k is reduced, but the downscaling algorithm still has a considerable bias (Fig.   3) during the historical period.

Performance of the downscaling algorithms for temperatures compared to ice cores
Compared to the temperature reconstructions, the downscaled CHELSA_TraCE21k model had reduced bias at four of the ice core sites located at the edges of the GIS (ReCAP, Agassiz, Hans Tausen Iskappe, Camp Century), but increased the bias, 330 RMSE and MAE at the remaining four sites at the center of the GIS (NEEM, NGRIP, GISP2, Dye 3) (Fig. 5). Overall both https://doi.org/10.5194/cp-2021-30 Preprint. Discussion started: 3 May 2021 c Author(s) 2021. CC BY 4.0 License.
CHELSA_TraCE as well as TraCE21k show a warm bias before the H1 event, and roughly after the 8.2k event at four of the sites (ReCAP, Agassiz, Hans Tausen Iskappe, Camp Century). At three sites (ReCAP, Agassiz, Hans Tausen Iskappe, Camp Century) a cold bias is present after the younger Dryas (Fig. 5). At the four other sites, the biases are idiosyncratic with CHELSA_TraCE usually showing a warm bias before the H1, and TraCE a cold bias before the H1 (Fig. 5). After the younger 335 dryas, both models show a cold bias at these sites. At the Camp Century site, the Trace data are close to the 15 N-based temperature reconstructions before the H1 event, and CHELSA_TraCE shows a warm bias, while after the younger dryas the situation is reversed (Fig.  5).

Accuracy of the glacier extent between 18kPB and 1kBP
The test validations of the glacial extent show a good performance over most time steps (Fig. 6), but with a notable drop in 345 accuracy at 8 kBP where all three validation metrics drop significantly. Aside from the drop at 8 kBP, the performance of the glacial extent simulations performed well. The marked drop in performance around 8 kBP might be due to the 8.2 kiloyear event, which marked a strong decrease in global temperatures, most likely due to meltwater fluxes from the collapsing Laurentide ice sheet. The strong coupling between temperature and glacial extent in CHELSA-TraCE21k generates an increase in glacial extent more than a sudden collapse during this time period, which seems to override the signal from the ICE6G 350 forcing data in CHELSA-TraCE21k.  cycle since the LGM, have had a significant influence on the distribution of ecosystems (Williams and Jackson, 2007), species (Hampe and Jump, 2011;Hewitt, 1999), and as a result on intraspecific genetic structures and speciation (Alsos et al., 2012;Pellissier et al., 2015;Yannic et al., 2014Yannic et al., , 2020. 365 Tracing the distribution of species through time is, however, challenging as the the spatio-temporal distributions of species strongly depend on environmental suitability (Guisan and Zimmermann, 2000), spatial accessibility of a given location (Normand et al., 2011;Svenning and Skov, 2004), and species dispersal abilities (Engler and Guisan, 2009). A dynamic simulation of species distributions can integrate all these aspects and therefore provides a valuable testbed for paleoclimatic data (Nobis and Normand, 2014). However, the spatio-temporal resolution of climate data needed for such simulations have 370 been limited to comparable coarse grain climatic data (Gherghel and Martin, 2020), which usually creates a mismatch between the climate derived from the model and the climate actually experienced by an organism (Seo et al., 2009).
Here, we use the downscaled transient temperature and precipitation data since 17kyBP (the coldest recorded temperatures in the TracCE21k model for Europe) to reconstruct refugia of the deciduous tree Grey Alder (Alnus incana) in Europe before post glacial climate warming. Similar to Nobis & Normand (2014), we first calibrated a generalized linear model (GLM) 375 (Nelder and Wedderburn, 1972) using current presences and absences of Grey Alder within polygons of the Atlas Florae Europaeae (AFE) (Kurtto et al., 2018) as the response variable and current annual mean temperature and precipitation from CHELSA-TraCE21k as predictors calculated as zonal mean values of 5 x 5 km rasterized AFE polygons. Despite the simplicity of the model it showed a fair to good model fit with a 10-fold cross-validated AUC value of 0.89.
Then, the GLM model was used to predict the suitability of Grey Alder from 17kBP till today in 500-year steps with 5 km 380 resolution and Lambert azimuthal equal-area projection. Glaciated areas were defined as unsuitable and were taken from the CHELSA-TraCE21k glacial reconstructions. We used the resulting time series of climatic suitability as input to the KISSMig model (Nobis & Normand 2014), which iteratively uses a simple 3 x 3 cell algorithm to calculate the spatial spread from a given origin from17kyBP to present. Presences and absences were weighted equally for the initial GLM calibration and KISSMig used squared suitability values to fulfill basic empirical expectations (see http://purl.oclc.org/wsl/kissmig). 385 We tested for each AFE polygon of the current Grey Alder distribution all 25 x 25 km areas across Europe as potential refugia.
All 5 x5 km grid cells of those areas suitable at 17kyBP were kept as refugia, if the respective AFE polygon was accessible and the spread pattern generated the lowest number of false positives when compared to the current AFE distribution. Because the migration ability of Grey Alder was unknown a priori, KISSMig simulations used one to 10 iterations for each 500-year step, corresponding to a maximum migration rate of 10 to 100 m/yr. For each iteration number the combined spread pattern 390 from all detected effective refugia was compared with the current distribution based on F1-scores. The optimized iteration number was identified by optimizing F1, which showed for Grey Alder a maximum migration rate of 50 m/yr. For a comparison with genetic clusters (Dering et al., 2016), the locations of that study were linked to the detected effective refugia with the shortest euclidean distance for simplicity.
Current genetic clustering of populations indicates that the modelling of A. incana distributions at 17kyBP successfully 395 detected glacial refugia in the Southern Alps, southern Norway, northern Norway, the Balkans, and eastern Romania (Fig. 7). The situation in eastern Europe is more complex, with most refugia located in Russia. However, since we only used the current distribution of A. incana in western Europe the results might be biased towards the east. Although almost the entire Europe was suitable for A. incana at 17kyBP, it might have not been at all locations due to dispersal constraints.

Conclusions 410
Although both, the original TraCE21k, as well as the downscaled CHELSA-Trace21k data track the relative temperature change well compared to ice cores, both models have relatively high temperature biases in absolute temperatures. Both the original data and the downscaled data have a warm bias before the younger Dryas, and a cold bias after it relative to the ice core proxy data. There are several reasons for this: coupled atmosphere ocean general circulation models (GCMs) such as CCSM3 cannot provide regional-scale or unbiased information on a variety of climatic processes (Meehl et al., 2007). 415 Temperatures from ice cores themselves are only based on proxy data, and the overall performance of such proxy data in estimating absolute temperatures is connected to biases in itself (Erb et al., 2018). The downscaling of the CHELSA-Trace21k data involves a trend-preserving (Hempel et al., 2013) change factor step to explicitly preserve the trends in TraCE21k. If however, these trends are already underestimated by the TraCE21k data, they will also be present in the downscaled data.
The estimation of glacial extents works comparably well when compared to expert delineations of the glacial extent of the 420 Laurentide ice shield. There is, however, a clear drop in performance at the 8.4kY event when atmospheric methane concentration decreased leading to a cooling and drying of the northern hemisphere (Kobashi et al., 2007). The strong coupling of the ice interpolations with only temperature might cause the decrease in performance as the downscaling algorithm ignores changes in precipitation which is only present in the driving ICE6G data. As the downscaling assumes an increase in glacial boundaries with cooling, this effect might not be realistic under an overall drying climate, and the fast shifts in temperatures 425 over only 150 years (Kobashi et al., 2007) might also not be well represented in a model with 100 year resolution. Another problem in the estimation of the glacial extent might be the result of the applied B-spline interpolation. The resulting ice cover from this interpolation can, in some areas, only be a few meters thick, not representing real glaciers, but rather spatial autocorrelation artefact of the interpolation approach used.
The CHELSA-Trace21k data seems to be able to recreate the distribution of temperature and precipitation in a meaningful 430 manner, so that the use of the data in subsequent analysis produces meaningful results. The reconstruction of the refugia for Alnus incana shows that the combination of high resolution climate data, together with a dynamic distribution model was able to accurately detect refugia, even those of a few kilometers in extent (Parducci et al., 2012) that cannot be detected using coarse climate data.

Code availability 435
Downscaling codes are based on (Karger et al., 2017a), and all modules used are open source and integrated into SAGA-GIS, available here: https://sourceforge.net/projects/saga-gis/. The code unique to this study is written in R and creates the paleoorography and glacier interpolations, and is also available on zenodo (DOI:10.5281/zenodo.4545753).