Articles | Volume 19, issue 2
Research article
20 Feb 2023
Research article |  | 20 Feb 2023

CHELSA-TraCE21k – high-resolution (1 km) downscaled transient temperature and precipitation data since the Last Glacial Maximum

Dirk Nikolaus Karger, Michael P. Nobis, Signe Normand, Catherine H. Graham, and Niklaus E. Zimmermann

High-resolution, downscaled climate model data are used in a wide variety of applications across environmental sciences. Here we introduce a new, high-resolution dataset, CHELSA-TraCE21k. It is obtained by downscaling TraCE-21k data, using the “Climatologies at high resolution for the earth's land surface areas” (CHELSA) V1.2 algorithm with the objective 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 at high spatial resolution and for each time step is created by combining high-resolution information on glacial cover from current and Last Glacial Maximum (LGM) glacier databases and interpolations using data from a global model of glacial isostasy (ICE-6G_C) and a coupling to mean annual temperatures from TraCE21k (Transient Climate Evolution of the last 21 000 years) based on the Community Climate System Model version 3 (CCSM3). Based on the reconstructed paleo-orography, mean annual temperature and precipitation were downscaled using the CHELSA V1.2 algorithm. The data were validated by comparisons with the glacial extent of the Laurentide ice sheet based on expert delineations, proxy data from Greenland ice cores, historical climate data from meteorological stations, and a dynamic simulation of species distributions throughout the Holocene. Validations show that the CHELSA-TraCE21k V1.0 dataset reasonably represents the distribution of temperature and precipitation through time at an unprecedented 1 km spatial resolution, and simulations based on the data are capable of detecting known LGM refugia of species.

1 Introduction

Since the Last Glacial Maximum (LGM), variation in climate has caused multiple changes in the Earth surface, including the rearrangement of species distributions or even species extinctions (Prentice et al., 1991; Velichko et al., 1997; Adams and Faure, 1997; Williams et al., 2004; Yu et al., 2010; Binney et al., 2017). 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 < 1 km (Seo et al., 2009) that are beyond the spatial resolution of numerical global circulation models (GCMs), which run at much coarser spatial resolution (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 modeling of high-resolution species distributions (Guisan and Zimmermann, 2000; Guisan and Thuiller, 2005), temporal and spatial variability in temperature and precipitation is of utmost importance. For such analyses, imprecisions 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 for many ecological applications has been bridged using statistical downscaling (Wilby et al., 1998; Wood et al., 2004; Schmidli et al., 2006; Maraun et al., 2010; Karger et al., 2017a, 2020), dynamical downscaling (Skamarock et al., 2021), or interpolation of meteorological station data (Daly et al., 1997; Hijmans et al., 2005; Meyer-Christoffer et al., 2015; Harris et al., 2020). While all of these methods work comparably well for current climatic conditions, station data are not available before the 19th century (end of the 20th century for satellite data), hampering an application of said methods to paleo-climatic models. Most paleo-climatic data at high spatial resolution are therefore based on climatologically aided interpolation (or the change factor method) 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 (Willmott and Robeson, 1995; Hunter and Meentemeyer, 2005). While this approach works rather well for short-term time series where topography is relatively stable (Daly et al., 1997), it becomes impractical 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 sheets and glaciers in polar regions 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 computationally very demanding, and therefore they have not been applied at ecologically relevant spatial resolutions of < 1 km yet. Current global kilometer-scale models only show a simulation throughput of 0.043 SYPD (simulated years per day) (Fuhrer et al., 2018), which is 25-fold lower than desired computationally efficient simulations of 1 SYPD (Schulthess et al., 2018; Schär et al., 2019). 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 output of numerical climate models (Frieler et al., 2017). Such studies therefore do not need 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; thus species distribution models (SDMs) are often built using a relatively small set of climate predictors based on monthly minimum and maximum temperature and precipitation (Guisan and Zimmermann, 2000; Guisan and Thuiller, 2005).

Here we present paleo-climatic data, downscaled from the CCSM3_TraCE21k (Transient climate evolution of the last 21 000 years using the Community Climate System Model Version 3) model output (hereafter: TraCE21k) to a 30 arcsec resolution using the CHELSA V1.2 (Climatologies at high resolution for the earth's land surface areas) algorithm (Karger et al., 2017a), which covers time steps of 100 years from 21 ka to 1950 plus four additional time steps until 1990 (TraCE21k), for minimum and maximum temperatures, surface precipitation, and paleo-orography.

2 Input data

2.1 TraCE-21k transient climate simulations

The TraCE-21k (Transient climate evolution of the last 21 000 years) simulation using the CCSM3 (Community Climate System Model Version 3) climate model (Liu et al., 2009; He, 2011; Marcott et al., 2011; Carlson et al., 2012) provides information on climate change over the last 21 000 years, i.e., from the Last Glacial Maximum (LGM, hereafter defined as 21 ka similar to Ehlers et al., 2011) to present. The TraCE-21k simulation reproduces many main features of post-glacial climate dynamics in various parts of the world from low to high latitudes and includes abrupt climate changes (Liu et al., 2009; He, 2011). The TraCE-21k simulation output has a T31_gx3v5 resolution (Otto-Bliesner et al., 2006). It uses a coarse-resolution dynamic global vegetation model (DGVM). The coupled atmosphere–ocean model in CCSM3 is based on the Community Atmospheric Model 3 (CAM3), on 26 vertical hybrid coordinate levels. The land and atmosphere components in CCSM3 in the TraCE-21k simulations use the same resolution. The parameterizations of the DGVM are largely based on the Lund–Potsdam–Jena (LPJ) DGVM. The ocean model in CCSM3 uses the NCAR (National Center for Atmospheric Research) version of the Parallel Ocean Program (POP) with 25 vertical levels, and the sea ice model is the NCAR Community Sea Ice Model (CSIM).

2.2 Observational climatology: CHELSA V1.2

CHELSA (Climatologies at high resolution for the earth's land surface areas) V1.2 is a high-resolution (30 arcsec) climate dataset for Earth's land surface areas (Karger et al., 2017b). It includes monthly means of daily 2 m mean, minimum, and maximum temperature and monthly precipitation rates at 30 arcsec resolution for the time period 1979–2013. CHELSA V1.2. is calculated with the CHELSA V1.2 topographic downscaling algorithm (Karger et al., 2017a), using the ERA-Interim (ECMWF Re-Analysis-Interim) reanalysis (Berrisford et al., 2009) as forcing data and GPCC (Global Precipitation Climatology Center) data (Meyer-Christoffer et al., 2015) for its bias correction.

2.3 Global model of glacial isostasy: ICE-6G_C (VM5a)

We used the output data of the ICE-6G_C (VM5a) (hereafter ICE6G) model as a basis for the extent of the major ice sheets at 1 resolution. ICE6G is a refinement of the ICE-5G (VM2) (hereafter ICE5G) global model of glacial isostasy (Peltier, 2004) which has been widely used to model the distribution of major ice sheets through time. ICE6G improves ICE5G by applying all available global positioning system (GPS) measurements of vertical motion of the crust that constrain the thickness of local ice cover as well as the timing of its removal. ICE6G explicitly outputs changes in ice thickness of major ice sheets (e.g., the Laurentide ice sheet) from the LGM till today (Argus et al., 2014; Peltier et al., 2015) at 500-year time steps.

2.4 Observational glacial extent at Last Glacial Maximum (LGM)

As the extent of the glaciers during the LGM, we use data from Ehlers et al. (2011) that present a detailed overview of Quaternary glaciations all over the world, with regards not only to stratigraphy but also to major glacial landforms and the extent of the respective ice sheets.

2.5 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 outlines were derived. Here we use this database to delineate the current extent of the glaciers at high resolution globally.

2.6 Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010)

The Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010) (Danielson and Gesch, 2011) dataset contains elevation data for the globe collected from various sources. Here we use the 30 arcsec version of the data that represents the mean elevation of all 7.5 arcsec grid cells that represent the highest available resolution of the data.

2.7 Bathymetric DEM

We use the General Bathymetric Chart of the Oceans (GEBCO) 2014 (Weatherall et al., 2015) as bathymetry. Although GEBCO also includes land surface altitude, we only use it for the oceans, and we keep as land altimetric data those of the CHELSA V1.2 algorithm (that being GMTED2010) to maintain comparable topography at the land surface.

2.8 Global sea level change

We used data from Miller et al. (2005) for the estimation of global sea level change from 21 ka to 1990. The data provide global estimates of sea level change over the last 100 million years. The entire time series of sea level change is based on a variety of proxy data, with the data used here dating back to the LGM, mainly based on tropical reef proxies (Miller et al., 2005).

3 Methods

Downscaling is based on the CHELSA V1.2 algorithm (Karger et al., 2017a) using forcing from TraCE-21k simulations (Liu et al., 2009; He, 2011) and involving several processing steps (Fig. 1). The CHELSA V1.2 algorithm needs a dynamic forcing in the form of GCM output (Karger et al., 2020) or gridded reanalysis data (Karger et al., 2017a, 2021b), as well as a surface orography (i.e., topography above sea level) to run a suite of downscaling algorithms for key climatic variables such as air temperature and precipitation. As the orography at different time steps between 21 ka and current times is not available at the high resolution required for the CHELSA V1.2 algorithm, we approximated it 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 for current conditions (Raup et al., 2007) and LGM conditions (Ehlers et al., 2011), and sea level change data from Miller et al. (2005) (Fig. 1). We then ran the CHELSA V1.2 algorithm on the paleo-orography using a bias-corrected version of the TraCE-21k simulations as a forcing. Details on these steps are described in the following sub-sections.

Figure 1Graphical representation of the different steps employed in downscaling TraCE-21k simulations using the CHELSA V1.2 algorithm. Input datasets are indicated by a blue box; output data are indicated by a green box. Rhombi indicate processing steps; t indicates discrete time steps, with t=0 being the LGM. Red lines indicate processing steps that are run iteratively over all time steps; black lines indicate computations that were run only once.


3.1 Paleo-orography

The first step in estimating the paleo-orography was carried out for the LGM (21 ka). For this time point, both estimates of glacial extents from Ehlers et al. (2011) and estimates of glacier thickness from ICE6G exist. We first combined the topographic information from GMTED2010 on land and that of GEBCO into a bedrock topography that provides the current bedrock topography ectopo (including current-day glaciers; see ff.). To create a bedrock orography etbed (i.e., topography adjusted for sea level without glaciers except for currently glaciated areas), we used the information on past sea level changes and set all elevations to 0 so that

(1) e t bed = e c topo - sl t .

To include the orography of the glaciers we first converted the polygons of the glacial extents from Ehlers et al. (2011) into point locations (Fig. 2a, black dots) and extracted their elevation from etbed (Fig. 2a) at time t=0 (LGM, 21 ka), resulting in the surface elevation of the glacial boundaries (gb) egbtbed. To combine the high-resolution estimates from Ehlers et al. (2011), with the coarser (1) resolution of ICE6G, we randomly sampled 100 point locations per 1 grid cell from ICE6G and again extracted the surface elevation of the glaciers from etbed at the ICE6G time step that is nearest to the time step t=0 (LGM, 21 ka), resulting in eGt (Fig. 2b). All points that did not fall within the high-resolution glacial extent were omitted so that only points within the high-resolution estimate of glacial extent from Ehlers et al. (2011) remained (Fig. 2b). Then both point datasets egbtbed and eICE6Gtorog were combined to represent a point sample of the surface elevation eGt within the high-resolution glacial extent of Ehlers et al. (2011) (Fig. 2c). Next, this point sample was spatially interpolated to a grid of 30 arcsec resolution applying multilevel B-spline interpolations. By this, we achieved an interpolated gap-free high-resolution estimate of glacial surface elevation etint at t=0 (LGM, 21 ka) (Fig. 2d). The multilevel B-splines use a B-spline approximation to eGt and start using the coarsest grid ϕ0 from an overall set of grids ϕ0, ϕ1, ϕn, with n=14 generated using optimized B-spline refinements (Lee et al., 1997). The resulting B-spline function f0(eGt) then gives the first approximation of etint=f0(eGt) and leaves a deviation

(2) Δ 1 e t int = e t int - f 0 ( x c , y c )

at each grid cell c location xc,yc,etint. Then the next control lattice ϕ1 is used to approximate f1(Δ1etintc). This approximation is then repeated on the sum of

(3) f 0 + f 1 = e t int - f 0 x c , y c - f 1 x c , y c

at each grid cell cxc,yc,etintcn times, resulting, in our case, in the gap-free interpolated glacial surface etint. The interpolated glacial surface was then combined with etbed to the topography et (Fig. 2e) using

(4) e t = e t topo , e t topo e t int e t int , otherwise .

The final orography etoro at time step t=0 (i.e., topography above sea level) (Fig. 2f) is then generated using

(5) e t oro = 0 , e t sl t e t , e t > sl t ,

with slt being the sea level at time step t. Although this approach includes changes in the glacial surface and sea level rise, it ignores changes in bedrock elevation due to upwelling after glacier melt.

3.2 Interpolation of glacier extent and thickness between LGM and current

As high-resolution estimates of glacial surface elevation are not available for time steps t other than the LGM and current day, we use a combination of mean annual 2 m air temperature data together with sea surface elevation slt and ICE6G orography eICE6Gtorog data to estimate etoro at each time step t≠0 and t=221. The rationale behind this approach is that temperature and glacier extents are interdependent, and a change in temperature will translate into a change in glacial extent or thickness. The procedure to generate high-resolution glacial surfaces is explained in the following sections.

3.2.1 Bias correction of air temperature

In a first step, the orography etoro with t=0 (LGM, 21 ka) was used to downscale mean annual temperature tast. GCMs such as the CCSM3 normally exhibit a large bias in temperatures or precipitation (Cannon et al., 2015; Maraun, 2013). We therefore applied a change factor bias correction based on the bias observed between current annual-mean 2 m air temperatures tascurobs from CHELSA V1.2 normals resampled to a 0.5 grid resolution and that of TraCE21k simulated for the same time period tascurmod, spline-interpolated using the same multilevel B-spline interpolation method as described in Sect. 3.1 to 0.5 grid resolution. The resolution of 0.5 follows the same procedure as used in the CHELSA V1.2 algorithm (Karger et al., 2017a). We used the time period 1980–1990 (= cur) to calculate this bias, as it is the only time period for which CHELSA V1.2 data and TraCE21k overlap. The change factor was then calculated as

(6) Δ tas = tas cur obs - tas cur mod .

This effectively preserves the trends observed in temperature but simultaneously assumes that the bias has also been conserved over time (Maraun, 2016). The bias-corrected temperatures tastcor are then given by

(7) tas t cor = tas cur obs - Δ tas .

3.2.2 Downscaling mean annual air temperature

To achieve a high-resolution approximation of near-surface air temperatures (Fig. 2g), 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 taz at altitudes az and vertical levels z 26 (992.5 hPa) to 20 (600.5 hPa) of the T31_gx3v5 grid that contain all surface elevations so that

(8) Γ = n ( a z t a z ) - ( a z ) ( t a z ) n ( a z 2 ) - ( a z ) 2 .

Temperature at the surface at a high spatial resolution (tas) was then calculated by

(9) tas t = Γ t e t oro + tas t cor .

3.2.3 Glacier extent approximation using mean annual air temperature

We assume that air temperature is correlated to glacier extent and use this relationship to estimate the boundaries of glaciers for each time step separately. To do so, we use mean annual 2 m air temperature at the boundary of the interpolated high-resolution glacier orography. We then transformed the glacier elevations etint to a polygon Gt and then transformed the outline of this polygon to a point sampling of the glacier boundaries at time t=0. Mean annual 2 m air temperatures at this glacier boundary were then extracted for this point sample, which gives the local annual-mean air temperature tast=0gb under which a glacier had a boundary at the LGM. To set this in relation to current mean annual air temperatures at current glacier boundaries tast=221gb (with t=221 being the year 1990), we calculated the difference between current and LGM boundary temperatures. The resulting point locations for both tast=0gb and tast=221gb were then spatially interpolated using a multilevel B-spline (as described in Sect. 3.1) to result in a gap-free surface and then subtracted, resulting in Δtascurobs (Fig. 2g).

As the orography for the next time step is not known yet, we estimated the near-surface air temperature tast+1est for the reduction in glacier extent similarly to the time step before so that

(10) tas t + 1 est = Γ t + 1 e t oro + tas t + 1 cor .

The binary glacial extent Gt+1 at t+1 is then approximated as

(11) G t + 1 = 1 , tas t + 1 est < tas t gb + Δ tas cur obs - 1 t - t G 0 = 1 0 , otherwise .

This glacial extent is then used again to estimate the combined topography at t+1 in the same way as described in Sect. 3.1. As ICE6G has a 500-year resolution we used the ICE6G orography that is closest to each time step. As ICE6G only includes information on the major ice sheets, smaller ice sheets in the Alps do not include a sample of eICE6Gtorog. In the case of smaller ice sheets, the surface orography from ICE6G is replaced by a point sample of the elevation of the glacier boundary under current conditions egbcorog (Fig. 2h). The glacier orography in this case is then created by using a spline interpolation between egbtbed and egbcorog. In Eq. (11) the second term in the condition for Gt+1 linearly scales Δtascurobs 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 unrealistically strong shift in the glacial extent. We then repeated the transformation of the glacial extent Gt+1 to all point locations and repeated the procedure for the temperature coupling to estimate the orography et+1oro. Near-surface air temperatures for t+1 have then been approximated using

(12) tas t + 1 = Γ t + 1 e t + 1 oro + tas t + 1 cor .

3.3 Downscaling mean monthly precipitation rates

3.3.1 Orographic wind effects

The estimation of high-resolution precipitation follows a variant of the CHELSA V1.2 algorithm (Karger et al., 2017a). The CHELSA V1.2 algorithm assumes that orography is one of the main drivers of precipitation (Spreen, 1947; Basist et al., 1994; Daly et al., 1997; Sevruk, 1997; Böhner, 2006; Gao et al., 2006; Böhner and Antonic, 2009; Karger et al., 2017a). In tropical convective regimes, precipitation typically increases up to the condensation level around 1000–1500 m above surface, while the exponentially decreasing moisture content in the mid- to upper troposphere results in a drying above the condensation level and 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 almost linear (Weischet and Endlicher, 2008). To approximate the effects of orographic precipitation we used the CHELSA V1.2 algorithm, which is explained in more detail below.

We used 10 m u-wind and v-wind components of TraCE-21k to calculate wind direction. Both wind components were projected to a world Mercator projection at a 4 km grid resolution using a multilevel B-spline interpolation similar to the one described in Sect. 3.1. Windward and leeward effects are assumed to be best represented at resolutions larger than 1 km (Daly et al., 1994); we therefore chose a grid resolution of 4 km for the underlying digital elevation model. The wind effect H was then calculated using


where dWHi and dLHi refer to the horizontal distances between the focal 4 km grid cell in the windward and leeward direction, and dWZi and dLZi are the corresponding vertical distances compared with the focal 4 km cell following the wind trajectory. The second summand in the equation for HW, where dLHi < 0, accounts for the leeward impact of previously traversed mountain chains. The horizontal distances in the equation for HL, where dLHi≥0, lead to a longer-distance impact of leeward rain shadow. The final wind effect parameter is calculated as H=HLHW. Both equations were applied to each grid cell at the 30 arcsec resolution in a world Mercator projection. Orographic precipitation effects are less pronounced just above the surface, as well as in the free atmosphere above the planetary boundary layer (Daly et al., 1997; Oke, 2002; Stull, 1988; Karger et al., 2020). 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 an 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 near-surface relative humidity (hurs) using

(15) LCL = 20 + ( tas / 5 ) ( 100 - hurs )

(Lawrence, 2005). The LCL has been interpolated to a 30 arcsec resolution using a B-spline interpolation. To create a boundary-layer-height-corrected wind effect HB, the wind effect grid H containing LCL was then proportionally distributed to all grid cells falling within a respective T31 grid cell using

(16) H B = H 1 - z - LCL z - z max h ,

with zmax being the maximum distance between the LCL at elevation z and all grid cells at a 30 arcsec resolution falling within a respective T31 grid cell. In Eq. (16), h is a constant of 9000 m, and z is the respective elevation from GMTED2010 (Danielson and Gesch, 2011) with

(17) LCL z = LCL + z GCM + f ,

zGCM being the elevation of the TraCE21k grid cell and f 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 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 (Liu et al., 2013). To account for these effects, we used a variant of the windward–leeward equations with a linear search distance of 300 km in circular steps of 5 from 0 to 355 for each grid cell. The calculated leeward index was then scaled towards higher elevations using

(18) E = i = 1 n 1 d WHi tan - 1 d LZi d WHi 0.5 i = 1 n 1 d LHi z h .

The h value has been set to 9000 m. EHB will give the first approximation of the orographic precipitation intensity pI.

3.3.2 Bias-correcting precipitation and downscaling

Precipitation, similar to temperature, exhibits a rather large bias in TraCE-21k (Figs. 3, 4). To remove this bias, we applied a change factor bias correction similar to the one described in Sect. 3.2.1 with the reference period 1980–1990. Here, we used a multiplicative change factor to avoid precipitation rates < 0. Additionally, we included a constant of c= 0.0001 kg m−2 per month to avoid division by zero so that

(19) Δ pr m = ( pr cur m mod + c ) / ( pr cur m obs + c ) ,

with m being the respective month of the year. The bias-corrected precipitation rate for pmcor is then calculated by

(20) pr m cor = pr cur m obs / Δ pr m .

To achieve the distribution of monthly precipitation pro given the approximated orographic precipitation intensity pIc at each grid location (xc,yc), we used a linear relationship between prmcor and prIc using

(21) pr o = p I c 1 n i = 1 n p I c i pr m cor ,

where n equals the number of 30 arcsec grid cells of pI that fall within a 0.5 grid cell of pmcor.

3.4 Downscaling mean monthly near-surface air temperatures

The downscaling of monthly near-surface air temperatures (tas, tasmax, tasmin) follows the methods described in Sect. 3.2.2, with the only difference being that instead of mean annual temperature, tasmax and tasmin are used, where tas = (tasmax + tasmin)/2. The temperatures have again first been bias-corrected using




with m being the respective month of the year in Eqs. (22)–(25).

Figure 2Illustration of several steps performed to estimate the surface orography and the temperature and precipitation fields in CHELSA-TraCE21k during the Last Glacial Maximum (21 ka). The upper row gives an example of the interpolation of the European ice sheets; the lower row shows an example of the resulting orography and environmental variables in the western part of the European Alps. (a) Topographic information at t=0 (LGM) is combined with information on past sea levels and the boundary of ice sheets (black dots) for which the surface elevation is extracted. (b) Within the extent of ice sheets, surface elevation is extracted from the ICE6G orography for a random sample of points for t=0. (c) Both point samples from (a) and (b) are combined and interpolated (d) to estimate the orography of the glaciers. (e) The interpolated glacier orography and the sea-level-adjusted topography are then combined. (f) The high-resolution (30 arcsec) orography (shown here for the western Alps) is then used as a basis at t=0 for (g) a lapse-rate-based downscaling of air temperature. (h) From the high-resolution temperatures, information on the glacier boundaries during the LGM (black) and current times (red) is extracted, and the difference is interpolated to correct the temperature-based shrinking and expansion of the glaciers. (i) Based on the orography the windward–leeward index is calculated (shown for July 21 ka), which builds the basis for the (j) precipitation approximation.

4 Output validation

Direct validation of the temperature (Fig. 1) and precipitation (Fig. 2) output at high resolution for paleo-time series 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 1 km paleo-climatic dataset. Therefore, to validate 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.

4.1 Validation using current (historical) observations

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 are available. To do so we calculated monthly climatologies for each month for tasmax, tasmin, and pr from both TraCE-21k and CHELSA-TraCE21k. We then compared the values measured at each station to those simulated in both TraCE-21k and CHELSA-TraCE21k.

The original TraCE-21k data show large deviations and root mean square errors (RMSEs) from the observed data (Fig. 3). This is expected as a climate model running for such long time periods needs to have coarse resolution, as well as a large degree of generality and realism, which decreases the accuracy of a model when compared to observations. The temperature variables perform well in TraCE-21k, with r∼0.8 for all months, but have deviations and RMSE similar to those of precipitation, which most likely can be attributed to the coarse resolution of the climate model. TraCE-21k also seems to overestimate temperature extremes for both tasmax and tasmin (Fig. 4).

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 CHELSA V1.2 algorithm improves the correlation between observed and modeled data 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 CHELSA-TraCE21k increases, which is reflected in an r value of ∼0.7 and a lower standard deviation and RMSE (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.

Figure 3Taylor diagrams comparing the relationship between TraCE-21k (blue) and CHELSA-TraCE21k (red). Data are shown for the 20th-century time period with average monthly observational data from the Global Historical Climate Network (GHCN) for the time period 1950–1990. Each dot represents a specific month.


Figure 4Scatterplots comparing precipitation, maximum, and minimum temperature. Data are aggregated from TraCE-21k and CHELSA-TraCE21k for the 20th-century time period with observational data from the Global Historical Climate Network (GHCN) for the time period 1950–1990.


4.2 Comparison with temperature proxies from ice core data

We compared the downscaled temperatures with the Greenland ice core reconstructions of Buizert et al. (2014, 2018) to check the performance of the downscaling at eight ice core locations on the Greenland Ice Sheet (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 21 ka to 1990 (Buizert et al., 2014, 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 TraCE-21k temperature data.

Compared to the temperature reconstructions from ice cores, 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, RMSE, and mean absolute error (MAE) at the remaining four sites at the center of the GIS (NEEM, NGRIP, GISP2, Dye 3) (Fig. 5). Overall, both CHELSA-TraCE21k and TraCE21k show a warm bias before the Heinrich 1 event (i.e., the break-off of large groups of icebergs from Greenland into the North Atlantic, 16.8 ka) and roughly after the 8.2 ka event at four of the sites (ReCAP, Agassiz, Hans Tausen Iskappe, Camp Century). At four sites (ReCAP, Agassiz, Hans Tausen Iskappe, Camp Century) a cold bias is present after the Younger Dryas (Fig. 5). At the four other sites, CHELSA-TraCE21k usually shows a warm bias before the H1 and TraCE-21k a cold bias before the H1 (Fig. 5). After the Younger Dryas (12.9–11.6 ka), both models show a cold bias at these sites. At the Camp Century site, the TraCE-21k data are close to the δ15N-based temperature reconstructions before the H1 event, and CHELSA-TraCE21k shows a warm bias, while after the Younger Dryas the situation is reversed (Fig. 5).

Figure 5Comparison of temperature anomalies from current times (1950–1990) for the CHELSA-TraCE21k time series data (blue), the TraCE21k data (red), and the temperature reconstructions from ice cores (black) for eight sites across Greenland. The gray horizontal line indicates the current observed temperature during the period 1981–2010 from CHELSA V2.1. data. Temperatures are plotted as anomalies from the current temperature recorded at the respective location of the ice cores.


The bias observed after downscaling might be related to biases in all the different input sources, such as the TraCE-21k bias being amplified, a bias in the ice core proxy data itself, or the bias correction using the simple change factor method. With the available data, these potential causes cannot be clearly disentangled but should be kept in mind for applications of the data.

4.3 Validation of glacier extent

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 exist.

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 data from Dyke (2004). The data consist of expertly delineated glacier maps based on a chronological database of radiocarbon dates and contain > 4000 dates located in North America (Dyke, 2004). To compare both datasets, we first calculated the glacial extent from CHELSA-TraCE21k by assigning a binary value to each 1 km grid cell in a Lambert conformal conic projection so that each data point compared equals 1 km2 either being covered by a glacier [1] or being free of a glacier [0]. We assigned a 1 if the simulated glacier height was above the paleo-terrain elevation and a 0 if it was lower than 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 day and that of the respective paleo-time step. 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 expertly delineated extent, we rasterized the polygons provided by Dyke (2004) for the years 18–1 ka to the 1 km resolution, extent, and projection of the simulated glacial cover and assign a 1 where the polygon intersects with a 1 km 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 balanced accuracy, which is defined as (sensitivity + specificity)/2. Additionally we report Cohen's kappa and the true-skill statistic (Allouche et al., 2006).

The test validations of the glacial extent show a good performance over most time steps (Fig. 6), but with a notable drop in accuracy at 8 ka, where all validation metrics drop significantly. Aside from the drop at 8 ka, the performance of the glacial-extent simulations performed well. The marked drop in performance around 8 ka might be due to the 8.2 ka 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 forcing data in CHELSA-TraCE21k. Additionally, we used the data from the extent of the ice sheets over Fennoscandia from 22 to 10 ka (Stroeven et al., 2016) for all time steps for which ICE6G data and data from Stroeven et al. (2016) were available. The results (Supplement Fig. S1) show, similar to the Laurentide ice sheet, that the accuracy is relatively high until 10.5 ka, with a drop in accuracy at 10 ka. Therefore, we assume that the temperature coupling does introduce errors in the time between 10 and 6 ka, as is evident from the comparison with the ice sheets of North America and over Fennoscandia.

Figure 6(a) Comparisons of estimated glacial extents of the Laurentide ice sheet from 16 to 11 ka. Blue delineates the interpolated ice sheet extent from CHELSA-TraCE21k, and red shows the estimated extent from Dyke (2004). While the retreat of the main Laurentide ice sheet is similar in both estimations, the Cordilleran ice sheet covering the Rocky Mountains retreats faster in the estimations by Dyke (2004) compared to CHELSA-TraCE21. (b) Performance comparison using three different metrics (balanced accuracy, Cohen's kappa, and true-skill statistic) from a comparison of CHELSA-TraCE21k and Dyke (2004).

5 Plausibility test using dynamic simulation of effective plant refugia

Transient long-term climatic data have a wide range of possible applications, ranging from population genetics (Leugger et al., 2022; Yannic et al., 2020), community ecology (Staples et al., 2022), and biodiversity buildup (Garcés-Pastor et al., 2022; Alsos et al., 2022) to evolutionary biology (Cerezer et al., 2022), just to name a few. Here we use one application in paleoecology as a plausibility test to additionally check if the transient CHELSA-TraCE21k data can reliably detect known LGM refugia of plant species. Climatic changes during the last glacial cycle since the LGM have had a significant influence on the distribution of ecosystems (Williams and Jackson, 2007), species (Hewitt, 1999; Hampe and Jump, 2011), and as a result on intraspecific genetic structures and speciation (Alsos et al., 2012; Yannic et al., 2014, 2020; Pellissier et al., 2015).

Tracing the distribution of species through time is, however, challenging as the spatio-temporal distributions of species strongly depend on environmental suitability (Guisan and Zimmermann, 2000), spatial accessibility of a given location (Svenning and Skov, 2004; Normand et al., 2011), and species dispersal abilities (Engler and Guisan, 2009). A dynamic simulation of species distributions can integrate all these aspects and therefore provides a valuable test bed for climatic data (Nobis and Normand, 2014). However, the spatio-temporal resolution of climate data needed for such simulations has 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 from CHELSA-TraCE21 since 17 ka (the coldest recorded temperatures in the CHELSA-TraCE21k model for Europe) to reconstruct refugia of the deciduous tree gray alder (Alnus incana) in Europe before post-glacial climate warming. Similar to Nobis and Normand (2014), we first calibrated a generalized linear model (GLM) (Nelder and Wedderburn, 1972) using current presences and absences of gray alder within polygons of the Atlas Florae Europaeae (AFE) (Jalas and Suominen, 1976) as the response variable and current annual-mean temperature and precipitation from CHELSA-TraCE21k as predictors calculated as zonal-mean values of 5×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 area-under-the-receiver-operating-characteristic-curve (AUC) value of 0.89.

Then, the GLM model was used to predict the suitability of gray alder from 17 ka till today in 500-year steps with 5 km 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 (Keep it simple stupid migration) model (Nobis and Normand, 2014), which iteratively uses a simple 3×3 cell algorithm to calculate the spatial spread from a given origin from 17 ka 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, last access: 30 December 2019).

We tested for each AFE polygon of the current gray alder distribution all 25×25 km areas across Europe as potential refugia. All 5×5 km grid cells of those areas suitable at 17 ka 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 gray alder was unknown a priori, KISSMig simulations used 1 to 10 iterations for each 500-year step, corresponding to a maximum migration rate of 10 to 100 m a−1. For each iteration number, the combined spread pattern 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 gray alder a maximum migration rate of 50 m a−1. 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 modeling of A. incana distributions at 17 ka shows that simulations based on CHELSA-TraCE21k successfully 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.

Figure 7Distribution of Alnus incana in Europe (based on the Atlas Flora Europaea; Jalas and Suominen, 1976) in current times (line shaded/hatched) and reconstructed effective refugia at 17 ka (dark-green alpha hull polygons) using dynamic species distribution modeling based on KISSMig and CHELSA-TraCE21k. The entire suitable habitat for A. incana at 17 ka is indicated as light green. Although almost all of south-central Europe was suitable for A. incana at 17 ka, it might not have occurred at all locations due to dispersal constraints, which are considered in the dark-green KISSMig reconstructed distribution. Colored circles indicate the population genetic structure of A. incana, taken from Dering et al. (2016), where each color represents a genetic cluster. Lines indicate the most likely effective refugia a genetic cluster can be associated with, given dispersal and climatic constraints. Current genetic clustering of populations indicates that the modeling of A. incana distributions at 17 ka successfully detected glacial effective refugia in the southern Alps (dark blue), southern Norway (light blue), northern Norway (pink), the Balkans (dark red), eastern Romania (turquoise), and the Black Sea (dark red and violet). As we only use the current distribution of A. incana within the AFE extent the results might be biased outside of it.

6 Conclusions

Although both the original TraCE-21k and 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). 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 themselves (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 TraCE-21k. If, however, these trends are already underestimated by the TraCE-21k data, they will also be present in the downscaled data.

The estimation of glacial extents shows an accuracy of > 80 % compared to expert delineations of the glacial extent of the Laurentide ice sheet. There is, however, a clear drop in accuracy at the 8 ka 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 that are 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 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 involve errors from 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 a spatial autocorrelation artifact of the interpolation approach used (e.g., see Fig. S1; 13 ka). Another source of error is that changes in bedrock due to the release of pressure from the melting ice sheets are not yet included in the algorithm. This can potentially result in several hundred meters of bias in affected areas that have not been taken into account in the current version of the algorithm.

The CHELSA-Trace21k data seem to be able to recreate the distribution of temperature and precipitation in a meaningful 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 with a dynamic distribution model was able to accurately detect refugia, even those of a few kilometers in extent (Parducci et al., 2012), which cannot be detected using coarse climate data.

Code availability

Downscaling codes are based on Karger et al. (2017a), and all modules used are open source and integrated into SAGA-GIS, available here: (Conrad and Wichmann, 2015). The code unique to this study is written in R and creates the paleo-orography and glacier interpolations and is also available on Zenodo (, greenmind1980, 2021).

Data availability

All post-processed data and additional input files other than those provided by TraCE21k can be accessed at (, Karger et al., 2021a). The data are published under a Creative Commons Attribution 2.0 Generic (CC BY 2.0) license.


The supplement related to this article is available online at:

Author contributions

DNK, MN, and NZ developed the idea. DNK developed the model and implemented the code. DNK, MN, and SN validated the data. NZ, SN, and CHG funded the project. DNK wrote the first version of the manuscript, and all authors contributed to subsequent revisions.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Financial support

This research has been supported by the Swiss Federal Institute for Forest, Snow and Landscape Research (ex-CHELSA and ClimEx grants); BiodivERsA COFUND (FeedBaCks grant); the Swiss National Science Foundation (grant nos. 20BD21_193907, 20BD21_184131); the ERA-NET BiodivERsA–Belmont Forum (FutureWeb grant); the Swiss Data Science Center (SPEEDMIND and COMECO grants); and the Aarhus University Research Foundation.

Review statement

This paper was edited by Irina Rogozhina and reviewed by two anonymous referees.


Adams, J. M. and Faure, H.: Preliminary Vegetation Maps of the World since the Last Glacial Maximum: An Aid to Archaeological Understanding, J. Archaeol. Sci., 24, 623–647,, 1997. 

Allouche, O., Tsoar, A., and Kadmon, R.: Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS), J. Appl. Ecol., 43, 1223–1232,, 2006. 

Alsos, I. G., Ehrich, D., Thuiller, W., Eidesen, P. B., Tribsch, A., Schönswetter, P., Lagaye, C., Taberlet, P., and Brochmann, C.: Genetic consequences of climate change for northern plants, P. R. Soc. B, 279, 2042–2051,, 2012. 

Alsos, I. G., Rijal, D. P., Ehrich, D., Karger, D. N., Yoccoz, N. G., Heintzman, P. D., Brown, A. G., Lammers, Y., Pellissier, L., Alm, T., Bråthen, K. A., Coissac, E., Merkel, M. K. F., Alberti, A., Denoeud, F., Bakke, J., and PHYLONORWAY CONSORTIUM: Postglacial species arrival and diversity buildup of northern ecosystems took millennia, Sci. Adv., 8, eabo7434,, 2022. 

Argus, D. F., Peltier, W. R., Drummond, R., and Moore, A. W.: The Antarctica component of postglacial rebound model ICE-6G_C (VM5a) based on GPS positioning, exposure age dating of ice thicknesses, and relative sea level histories, Geophys. J. Int., 198, 537–563,, 2014. 

Basist, A., Bell, G. D., and Meentemeyer, V.: Statistical Relationships between Topography and Precipitation Patterns, J. Climate, 7, 1305–1315,<1305:SRBTAP>2.0.CO;2, 1994. 

Berrisford, P., Dee, D., Fielding, K., Fuentes, M., Kallberg, P., Kobayashi, S., and Uppala, S.: The ERA-interim archive, ERA Rep. Ser., 1, 1–16, 2009. 

Binney, H., Edwards, M., Macias-Fauria, M., Lozhkin, A., Anderson, P., Kaplan, J. O., Andreev, A., Bezrukova, E., Blyakharchuk, T., Jankovska, V., Khazina, I., Krivonogov, S., Kremenetski, K., Nield, J., Novenko, E., Ryabogina, N., Solovieva, N., Willis, K., and Zernitskaya, V.: Vegetation of Eurasia from the last glacial maximum to present: Key biogeographic patterns, Quaternary Sci. Rev., 157, 80–97,, 2017. 

Böhner, J.: General climatic controls and topoclimatic variations in Central and High Asia, Boreas, 35, 279–295,, 2006. 

Böhner, J. and Antonic, O.: Land-Surface Parameters Specific to Topo-Climatology, in: GEOMORPHOMETRY: CONCEPTS, SOFTWARE, APPLICATIONS, edited by: Hengl, T. and Reuter, H. I., Geomorphometry: Concepts, Software, Applications, Elsevier Science, 195–226,, 2009. 

Brown, J. L., Hill, D. J., Dolan, A. M., Carnaval, A. C., and Haywood, A. M.: PaleoClim, high spatial resolution paleoclimate surfaces for global land areas, Sci. Data, 5, 1–9,, 2018. 

Buizert, C., Gkinis, V., Severinghaus, J. P., He, F., Lecavalier, B. S., Kindler, P., Leuenberger, M., Carlson, A. E., Vinther, B., Masson-Delmotte, V., White, J. W. C., Liu, Z., Otto-Bliesner, B., and Brook, E. J.: Greenland temperature response to climate forcing during the last deglaciation, Science, 345, 1177–1180,, 2014. 

Buizert, C., Keisling, B. A., Box, J. E., He, F., Carlson, A. E., Sinclair, G., and DeConto, R. M.: Greenland-Wide Seasonal Temperatures During the Last Deglaciation, Geophys. Res. Lett., 45, 1905–1914,, 2018. 

Cannon, A. J., Sobie, S. R., and Murdock, T. Q.: Bias Correction of GCM Precipitation by Quantile Mapping: How Well Do Methods Preserve Changes in Quantiles and Extremes?, J. Climate, 28, 6938–6959,, 2015. 

Carlson, A. E., Ullman, D. J., Anslow, F. S., He, F., Clark, P. U., Liu, Z., and Otto-Bliesner, B. L.: Modeling the surface mass-balance response of the Laurentide Ice Sheet to Bølling warming and its contribution to Meltwater Pulse 1A, Earth Planet. Sc. Lett., 315–316, 24–29,, 2012. 

Cerezer, F. O., Machac, A., Rangel, T. F., and Dambros, C. S.: Exceptions to the rule: Relative roles of time, diversification rates and regional energy in shaping the inverse latitudinal diversity gradient, Glob. Ecol. Biogeogr., 31, 1794–1809,, 2022. 

Conrad, O. and Wichmann, V.: SAGA GIS, (last access: 16 September 2018), 2015. 

Daly, C., Neilson, R. P., and Phillips, D. L.: A Statistical-Topographic Model for Mapping Climatological Precipitation over Mountainous Terrain, J. Appl. Meteorol., 33, 140–158,<0140:ASTMFM>2.0.CO;2, 1994. 

Daly, C., Taylor, G. H., and Gibson, W. P.: The PRISM approach to mapping precipitation and temperature, Proc. 10th AMS Conf Appl. Climatol., 20–23, (last access: 29 August 2018), 1997. 

Danielson, J. J. and Gesch, D. B.: Global multi-resolution terrain elevation data 2010 (GMTED2010): U.S. Geo-logical Survey Open-File Report 2011–1073, 26 pp., (last access: 29 August 2018), 2011. 

Dering, M., Latałowa, M., Boratyńska, K., Kosiński, P., and Boratyński, A.: Could clonality contribute to the northern survival of grey alder [Alnus incana (L.) Moench] during the Last Glacial Maximum?, Acta Soc. Bot. Pol., 86, 1–14,, 2016. 

Dyke, A. S.: An outline of North American deglaciation with emphasis on central and northern Canada, in: Developments in Quaternary Sciences, vol. 2, edited by: Ehlers, J. and Gibbard, P. L., Elsevier, 373–424,, 2004. 

Ehlers, J., Gibbard, P. L., and Hughes, P. D.: Quaternary Glaciations – Extent and Chronology, Volume 15, 1st Edition, ISBN 9780444534477, 2011. 

Engler, R. and Guisan, A.: MigClim: Predicting plant distribution and dispersal in a changing climate, Divers. Distrib., 15, 590–601,, 2009. 

Erb, M. P., Jackson, C. S., Broccoli, A. J., Lea, D. W., Valdes, P. J., Crucifix, M., and DiNezio, P. N.: Model evidence for a seasonal bias in Antarctic ice cores, Nat. Commun., 9, 1361,, 2018. 

Frei, C. and Schär, C.: A precipitation climatology of the Alps from high-resolution rain-gauge observations, Int. J. Climatol., 18, 873–900,<873::AID-JOC255>3.0.CO;2-9, 1998. 

Frieler, K., Lange, S., Piontek, F., Reyer, C. P. O., Schewe, J., Warszawski, L., Zhao, F., Chini, L., Denvil, S., Emanuel, K., Geiger, T., Halladay, K., Hurtt, G., Mengel, M., Murakami, D., Ostberg, S., Popp, A., Riva, R., Stevanovic, M., Suzuki, T., Volkholz, J., Burke, E., Ciais, P., Ebi, K., Eddy, T. D., Elliott, J., Galbraith, E., Gosling, S. N., Hattermann, F., Hickler, T., Hinkel, J., Hof, C., Huber, V., Jägermeyr, J., Krysanova, V., Marcé, R., Müller Schmied, H., Mouratiadou, I., Pierson, D., Tittensor, D. P., Vautard, R., van Vliet, M., Biber, M. F., Betts, R. A., Bodirsky, B. L., Deryng, D., Frolking, S., Jones, C. D., Lotze, H. K., Lotze-Campen, H., Sahajpal, R., Thonicke, K., Tian, H., and Yamagata, Y.: Assessing the impacts of 1.5  C global warming – simulation protocol of the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP2b), Geosci. Model Dev., 10, 4321–4345,, 2017. 

Fuhrer, O., Chadha, T., Hoefler, T., Kwasniewski, G., Lapillonne, X., Leutwyler, D., Lüthi, D., Osuna, C., Schär, C., Schulthess, T. C., and Vogt, H.: Near-global climate simulation at 1 km resolution: establishing a performance baseline on 4888 GPUs with COSMO 5.0, Geosci. Model Dev., 11, 1665–1681,, 2018. 

Gao, X., Xu, Y., Zhao, Z., Pal, J. S., and Giorgi, F.: On the role of resolution and topography in the simulation of East Asia precipitation, Theor. Appl. Climatol., 86, 173–185,, 2006. 

Garcés-Pastor, S., Coissac, E., Lavergne, S., Schwörer, C., Theurillat, J.-P., Heintzman, P. D., Wangensteen, O. S., Tinner, W., Rey, F., Heer, M., Rutzer, A., Walsh, K., Lammers, Y., Brown, A. G., Goslar, T., Rijal, D. P., Karger, D. N., Pellissier, L., Heiri, O., and Alsos, I. G.: High resolution ancient sedimentary DNA shows that alpine plant diversity is associated with human land use and climate change, Nat. Commun., 13, 6559,, 2022. 

Gherghel, I. and Martin, R. A.: Postglacial recolonization of North America by spadefoot toads: integrating niche and corridor modeling to study species' range dynamics over geologic time, Ecography, 43, 1499–1509,, 2020. 

greenmind1980: greenmind1980/CHELSA_TraCE21k: Version 1.0 (V1.0), Zenodo [code],, 2021. 

Guisan, A. and Thuiller, W.: Predicting species distribution: offering more than simple habitat models, Ecol. Lett., 8, 993–1009, 2005. 

Guisan, A. and Zimmermann, N. E.: Predictive habitat distribution models in ecology, Ecol. Model., 135, 147–186, 2000. 

Hampe, A. and Jump, A. S.: Climate Relicts: Past, Present, Future, Annu. Rev. Ecol. Evol. S., 42, 313–333,, 2011. 

Harris, I., Osborn, T. J., Jones, P., and Lister, D.: Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset, Sci. Data, 7, 1–18,, 2020. 

He, F.: Simulating Transient Climate Evolution of the Last Deglaciation with CCSM3, PhD Thesis, University of Wisconsin Madison, Madison, WC, USA, 171 pp., (last access: 2 February 2017), 2011. 

Hempel, S., Frieler, K., Warszawski, L., Schewe, J., and Piontek, F.: A trend-preserving bias correction – the ISI-MIP approach, Earth Syst. Dynam., 4, 219–236,, 2013. 

Hewitt, G. M.: Post-glacial re-colonization of European biota, Biol. J. Linn. Soc., 68, 87–112,, 1999. 

Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G., and Jarvis, A.: Very high resolution interpolated climate surfaces for global land areas, Int. J. Climatol., 25, 1965–1978,, 2005. 

Hunter, R. D. and Meentemeyer, R. K.: Climatologically Aided Mapping of Daily Precipitation and Temperature, J. Appl. Meteorol., 44, 1501–1510,, 2005. 

Hutchinson, G. E.: Population Studies: Animal Ecology and Demography – Concluding Remarks, Cold Spring Harb. Sym., 22, 415–427,, 1957. 

Jalas, J. and Suominen, J. (Eds.): Atlas Florae Europaeae. Distribution of Vascular Plants in Europe. 3. Salicaceae to Balanophoraceae. – The Committee for Mapping the Flora of Europe & Societas Biologica Fennica Vanamo, Helsinki, 128 pp., ISBN 951-9108-02-5, 1976. 

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., Zimmermann, N. E., Linder, H. P., and Kessler, M.: Climatologies at high resolution for the earth's land surface areas, Sci. Data, 4, 170122,, 2017a. 

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., Zimmermann, N. E., Linder, H. P., and Kessler, M.: Climatologies at high resolution for the earth's land surface areas, Dryad Digital Repository [data set],, 2017b. 

Karger, D. N., Schmatz, D. R., Dettling, G., and Zimmermann, N. E.: High resolution monthly precipitation and temperature timeseries for the period 2006–2100, Sci. Data, 7, 248,, 2020. 

Karger, D. N., Nobis, M., Normand, S., Graham, C. H., and Zimmermann, N. E.: CHELSA-TraCE21k: Downscaled transient temperature and precipitation data since the last glacial maximum – EnviDat, envidat [data set],, 2021a. 

Karger, D. N., Wilson, A. M., Mahony, C., Zimmermann, N. E., and Jetz, W.: Global daily 1 km land surface precipitation based on cloud cover-informed downscaling, Sci. Data, 8, 307,, 2021b. 

Kobashi, T., Severinghaus, J. P., Brook, E. J., Barnola, J.-M., and Grachev, A. M.: Precise timing and characterization of abrupt climate change 8200 years ago from air trapped in polar ice, Quaternary Sci. Rev., 26, 1212–1222,, 2007. 

Körner, C.: The use of “altitude” in ecological research, Trends Ecol. Evol., 22, 569–574,, 2007. 

Lawrence, M. G.: The Relationship between Relative Humidity and the Dewpoint Temperature in Moist Air: A Simple Conversion and Applications, B. Am. Meteorol. Soc., 86, 225–234,, 2005. 

Lawrimore, J. H., Menne, M. J., Gleason, B. E., Williams, C. N., Wuertz, D. B., Vose, R. S., and Rennie, J.: An overview of the Global Historical Climatology Network monthly mean temperature data set, version 3, J. Geophys. Res.-Atmos., 116, D19121,, 2011. 

Leugger, F., Broquet, T., Karger, D. N., Rioux, D., Buzan, E., Corlatti, L., Crestanello, B., Curt-Grand-Gaudin, N., Hauffe, H. C., Rolečková, B., Šprem, N., Tissot, N., Tissot, S., Valterová, R., Yannic, G., and Pellissier, L.: Dispersal and habitat dynamics shape the genetic structure of the Northern chamois in the Alps, J. Biogeogr., 49, 1848–1861,, 2022. 

Liu, M., Bárdossy, A., and Zehe, E.: Interaction of valleys and circulation patterns (CPs) on spatial precipitation patterns in southern Germany, Hydrol. Earth Syst. Sci., 17, 4685–4699,, 2013. 

Liu, Z., Otto-Bliesner, B. L., He, F., Brady, E. C., Tomas, R., Clark, P. U., Carlson, A. E., Lynch-Stieglitz, J., Curry, W., Brook, E., Erickson, D., Jacob, R., Kutzbach, J., and Cheng, J.: Transient Simulation of Last Deglaciation with a New Mechanism for Bølling-Allerød Warming, Science, 325, 310–314,, 2009. 

Maraun, D.: Bias correction, quantile mapping, and downscaling: Revisiting the inflation issue, J. Climate, 26, 2137–2143, 2013. 

Maraun, D.: Bias Correcting Climate Change Simulations – a Critical Review, Curr. Clim. Change Rep., 2, 211–220,, 2016. 

Maraun, D., Wetterhall, F., Ireson, A. M., Chandler, R. E., Kendon, E. J., Widmann, M., Brienen, S., Rust, H. W., Sauter, T., Themeßl, M., Venema, V. K. C., Chun, K. P., Goodess, C. M., Jones, R. G., Onof, C., Vrac, M., and Thiele-Eich, I.: Precipitation downscaling under climate change: Recent developments to bridge the gap between dynamical models and the end user, Rev. Geophys., 48, RG3003,, 2010. 

Marcott, S. A., Clark, P. U., Padman, L., Klinkhammer, G. P., Springer, S. R., Liu, Z., Otto-Bliesner, B. L., Carlson, A. E., Ungerer, A., Padman, J., He, F., Cheng, J., and Schmittner, A.: Ice-shelf collapse from subsurface warming as a trigger for Heinrich events, P. Natl. Acad. Sci. USA, 108, 13415–13419,, 2011. 

McMaster, G. S. and Wilhelm, W. W.: Growing degree-days: one equation, two interpretations, Agr. Forest Meteorol., 87, 291–300,, 1997. 

Meehl, G. A., Stocker, T. F., Collins, W. D., Friedlingstein, P., Gaye, A. T., Gregory, J. M., Kitoh, A., Knutti, R., Murphy, J. M., Noda, A., Raper, S. C. B., Watterson, I. G., Weaver, A. J., and Zhao, Z.-C.: Global climate projections, in: Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, UK, ISBN 978 0521 88009-1, 2007. 

Meyer-Christoffer, A., Becker, A., Finger, P., Rudolf, B., Schneider, U., and Ziese, M.: GPCC Climatology Version 2015 at 0.25: Monthly Land-Surface Precipitation Climatology for Every Month and the Total Year from Rain-Gauges built on GTS-based and Historic Data., Glob. Precip. Climatol. Cent. GPCC,, 2015. 

Miller, K. G., Kominz, M. A., Browning, J. V., Wright, J. D., Mountain, G. S., Katz, M. E., Sugarman, P. J., Cramer, B. S., Christie-Blick, N., and Pekar, S. F.: The Phanerozoic Record of Global Sea-Level Change, Science, 310, 1293–1298,, 2005. 

Nelder, J. A. and Wedderburn, R. W. M.: Generalized Linear Models, J. R. Stat. Soc. Ser. A-Gen., 135, 370–384,, 1972. 

Neumann, P., Düben, P., Adamidis, P., Bauer, P., Brück, M., Kornblueh, L., Klocke, D., Stevens, B., Wedi, N., and Biercamp, J.: Assessing the scales in numerical weather and climate predictions: will exascale be the rescue?, Philos. T. R. Soc. A, 377, 20180148,, 2019. 

Nobis, M. P. and Normand, S.: KISSMig – a simple model for R to account for limited migration in analyses of species distributions, Ecography, 37, 1282–1287,, 2014. 

Normand, S., Ricklefs, R. E., Skov, F., Bladt, J., Tackenberg, O., and Svenning, J.-C.: Postglacial migration supplements climate in determining plant species ranges in Europe, Philos. T. R. Soc. B, 278, 3644–3653,, 2011. 

Oke, T. R.: Boundary layer climates, Routledge, 464 pp., ISBN 9780415043199, 2002. 

Otto-Bliesner, B. L., Brady, E. C., Clauzet, G., Tomas, R., Levis, S., and Kothavala, Z.: Last Glacial Maximum and Holocene Climate in CCSM3, J. Climate, 19, 2526–2544,, 2006. 

Parducci, L., Jørgensen, T., Tollefsrud, M. M., Elverland, E., Alm, T., Fontana, S. L., Bennett, K. D., Haile, J., Matetovici, I., Suyama, Y., Edwards, M. E., Andersen, K., Rasmussen, M., Boessenkool, S., Coissac, E., Brochmann, C., Taberlet, P., Houmark-Nielsen, M., Larsen, N. K., Orlando, L., Gilbert, M. T. P., Kjær, K. H., Alsos, I. G., and Willerslev, E.: Glacial Survival of Boreal Trees in Northern Scandinavia, Science, 335, 1083–1086,, 2012. 

Pellissier, L., Eidesen, P. B., Ehrich, D., Descombes, P., Schönswetter, P., Tribsch, A., Westergaard, K. B., Alvarez, N., Guisan, A., Zimmermann, N. E., Normand, S., Vittoz, P., Luoto, M., Damgaard, C., Brochmann, C., Wisz, M. S., and Alsos, I. G.: Past climate-driven range shifts and population genetic diversity in arctic plants, J. Biogeogr., 43, 461–470,, 2015. 

Peltier, W. R.: Global glacial isostasy and the surface of the ice-age earth: The ICE-5G (CM2) Model and GRACE, Annu. Rev. Earth Pl. Sc., 32, 111–149,, 2004. 

Peltier, W. R., Argus, D. F., and Drummond, R.: Space geodesy constrains ice age terminal deglaciation: The global ICE-6G_C (VM5a) model, J. Geophys. Res.-Sol. Ea., 120, 450–487,, 2015. 

Prentice, I. C., Bartlein, P. J., and Webb, T.: Vegetation and Climate Change in Eastern North America Since the Last Glacial Maximum, Ecology, 72, 2038–2056,, 1991. 

Raup, B., Racoviteanu, A., Khalsa, S. J. S., Helm, C., Armstrong, R., and Arnaud, Y.: The GLIMS geospatial glacier database: A new tool for studying glacier change, Global Planet. Change, 56, 101–110,, 2007. 

Rotunno, R. and Houze, R. A.: Lessons on orographic precipitation from the Mesoscale Alpine Programme, Q. J. R. Meteor. Soc., 133, 811–830,, 2007. 

Schär, C., Fuhrer, O., Arteaga, A., Ban, N., Charpilloz, C., Di Girolamo, S., Hentgen, L., Hoefler, T., Lapillonne, X., Leutwyler, D., Osterried, K., Panosetti, D., Rüdisühli, S., Schlemmer, L., Schulthess, T., Sprenger, M., Ubbiali, S., and Wernli, H.: Kilometer-scale climate models: Prospects and challenges, B. Am. Meteorol. Soc., 101,, 2019. 

Schmidli, J., Frei, C., and Vidale, P. L.: Downscaling from GCM precipitation: a benchmark for dynamical and statistical downscaling methods, Int. J. Climatol., 26, 679–689,, 2006. 

Schulthess, T. C., Bauer, P., Wedi, N., Fuhrer, O., Hoefler, T., and Schär, C.: Reflecting on the goal and baseline for exascale computing: a roadmap based on weather and climate simulations, Comput. Sci. Eng., 21, 30–41, 2018. 

Scotese, C. R.: Atlas of earth history, PALEOMAP project, (last access: 16 September 2018), 2001. 

Seo, C., Thorne, J. H., Hannah, L., and Thuiller, W.: Scale effects in species distribution models: implications for conservation planning under climate change, Biol. Lett., 5, 39–43,, 2009. 

Sepulchre, P., Caubel, A., Ladant, J.-B., Bopp, L., Boucher, O., Braconnot, P., Brockmann, P., Cozic, A., Donnadieu, Y., Dufresne, J.-L., Estella-Perez, V., Ethé, C., Fluteau, F., Foujols, M.-A., Gastineau, G., Ghattas, J., Hauglustaine, D., Hourdin, F., Kageyama, M., Khodri, M., Marti, O., Meurdesoif, Y., Mignot, J., Sarr, A.-C., Servonnat, J., Swingedouw, D., Szopa, S., and Tardif, D.: IPSL-CM5A2 – an Earth system model designed for multi-millennial climate simulations, Geosci. Model Dev., 13, 3011–3053,, 2020. 

Sevruk, B.: Regional Dependency of Precipitation-Altitude Relationship in the Swiss Alps, in: Climatic Change at High Elevation Sites, edited by: Diaz, H. F., Beniston, M., and Bradley, R. S., Springer, the Netherlands, 123–137,, 1997. 

Skamarock, C., Klemp, B., Dudhia, J., Gill, O., Liu, Z., Berner, J., Wang, W., Powers, G., Duda, G., Barker, D., and Huang, X.: A Description of the Advanced Research WRF Model Version 4.3, No. NCAR/TN-556+STR,, 2021. 

Soria-Auza, R. W., Kessler, M., Bach, K., Barajas-Barbosa, P. M., Lehnert, M., Herzog, S. K., and Bohner, J.: Impact of the quality of climate models for modelling species occurrences in countries with poor climatic documentation: a case study from Bolivia, Ecol. Model., 221, 1221–1229, 2010. 

Spreen, W. C.: A determination of the effect of topography upon precipitation, Eos T. Am. Geophys. Un., 28, 285–290,, 1947. 

Staples, T. L., Kiessling, W., and Pandolfi, J. M.: Emergence patterns of locally novel plant communities driven by past climate change and modern anthropogenic impacts, Ecol. Lett., 25, 1497–1509,, 2022. 

Stroeven, A. P., Hättestrand, C., Kleman, J., Heyman, J., Fabel, D., Fredin, O., Goodfellow, B. W., Harbor, J. M., Jansen, J. D., Olsen, L., Caffee, M. W., Fink, D., Lundqvist, J., Rosqvist, G. C., Strömberg, B., and Jansson, K. N.: Deglaciation of Fennoscandia, Quaternary Sci. Rev., 147, 91–121,, 2016. 

Stull, R. B. (Ed.): An Introduction to Boundary Layer Meteorology, Springer Netherlands, Dordrecht,, 1988. 

Svenning, J.-C. and Skov, F.: Limited filling of the potential range in European tree species, Ecol. Lett., 7, 565–573,, 2004. 

Velichko, A. A., Andreev, A. A., and Klimanov, V. A.: Climate and vegetation dynamics in the tundra and forest zone during the late glacial and holocene, Quatern. Int., 41–42, 71–96,, 1997. 

Weatherall, P., Marks, K. M., Jakobsson, M., Schmitt, T., Tani, S., Arndt, J. E., Rovere, M., Chayes, D., Ferrini, V., and Wigley, R.: A new digital bathymetric model of the world's oceans, Earth Space Sci., 2, 331–345,, 2015. 

Weischet, W. and Endlicher, W.: Einführung in die Allgemeine Klimatologie, Schweizerbart Science Publishers, Stuttgart, Germany, 342 pp., ISBN 978-3-443-07155-4, 2008. 

Wilby, R. L., Wigley, T. M. L., Conway, D., Jones, P. D., Hewitson, B. C., Main, J., and Wilks, D. S.: Statistical downscaling of general circulation model output: A comparison of methods, Water Resour. Res., 34, 2995–3008,, 1998. 

Williams, J. W. and Jackson, S. T.: Novel climates, no-analog communities, and ecological surprises, Front. Ecol. Environ., 5, 475–482, 2007. 

Williams, J. W., Shuman, B. N., III, T. W., Bartlein, P. J., and Leduc, P. L.: Late-Quaternary Vegetation Dynamics in North America: Scaling from Taxa to Biomes, Ecol. Monogr., 74, 309–334, 2004. 

Willmott, C. J. and Robeson, S. M.: Climatologically aided interpolation (CAI) of terrestrial air temperature, Int. J. Climatol., 15, 221–229,, 1995. 

Wood, A. W., Leung, L. R., Sridhar, V., and Lettenmaier, D. P.: Hydrologic Implications of Dynamical and Statistical Approaches to Downscaling Climate Model Outputs, Clim. Change, 62, 189–216,, 2004. 

Woodward, F. I., Fogg, G. E., Heber, U., Laws, R. M., and Franks, F.: The impact of low temperatures in controlling the geographical distribution of plants, Philos. T. R. Soc. B, 326, 585–593,, 1990. 

Yannic, G., Pellissier, L., Ortego, J., Lecomte, N., Couturier, S., Cuyler, C., Dussault, C., Hundertmark, K. J., Irvine, R. J., Jenkins, D. A., Kolpashikov, L., Mager, K., Musiani, M., Parker, K. L., Røed, K. H., Sipko, T., Þórisson, S. G., Weckworth, B. V., Guisan, A., Bernatchez, L., and Côté, S. D.: Genetic diversity in caribou linked to past and future climate change, Nat. Clim. Change, 4, 132–137,, 2014.  

Yannic, G., Hagen, O., Leugger, F., Karger, D. N., and Pellissier, L.: Harnessing paleo-environmental modeling and genetic data to predict intraspecific genetic structure, Evol. Appl., 13, 1526–1542,, 2020. 

Yu, Z., Loisel, J., Brosseau, D. P., Beilman, D. W., and Hunt, S. J.: Global peatland dynamics since the Last Glacial Maximum, Geophys. Res. Lett., 37, L13402,, 2010. 

Short summary
Here we present global monthly climate time series for air temperature and precipitation at 1 km resolution for the last 21 000 years. The topography at all time steps is created by combining high-resolution information on glacial cover from current and Last Glacial Maximum glacier databases with the interpolation of an ice sheet model and a coupling to mean annual temperatures from a global circulation model.