Species Favourability Shift in Europe due to Climate Change : A Case Study for Fagus sylvatica L . and Picea abies ( L . ) Karst . Based on an Ensemble of Climate Models

Climate is the main environmental driver determining the spatial distribution of most tree species at the continental scale. We investigated the distribution change of European beech andNorway spruce due to climate change.We applied a species distribution model (SDM), driven by an ensemble of 21 regional climate models in order to study the shift of the favourability distribution of these species. SDMs were parameterized for 1971–2000, as well as 2021–2050 and 2071–2100 using the SRES scenario A1B and three physiological meaningful climate variables. Growing degree sum and precipitation sum were calculated for the growing season on a basis of daily data. Results show a general north-eastern and altitudinal shift in climatological favourability for both species, although the shift is more marked for spruce. The gain of new favourable sites in the north or in the Alps is stronger for beech compared to spruce. Uncertainty is expressed as the variance of the averaged maps and with a density function. Uncertainty in species distribution increases over time. This study demonstrates the importance of data ensembles and shows how to deal with different outcomes in order to improve impact studies by showing uncertainty of the resulting maps.


Introduction
Climate change affects tree species because trees are sessile, long lived, and comparatively slow in reacting to a changing environment, for example, trees cannot easily colonize new habitats in a strongly fragmented landscape.Because of the slow growth rates of trees to reach maturity, the typical planning horizon in the forestry sector is, for example, for spruce 60 to 80 years in the future.Hence, forestry has to deal with the challenges of a rapidly changing environment in order to avoid breakdown of stands.The same is true for city planners or municipal authorities with regard to tree species.Therefore, forest owners and related professionals need a reliable assessment about shifts of species distribution in order to make management decisions [1].
Even though genetic resources are considered as a key factor of forest ecosystems to adapt to climate change [2], one has to keep in mind that adaptation processes are slow, that niche conservatism is described for some genera [3], and that species will not survive or at least will not have a positive population growth rate if environmental conditions at a site will not meet the physiological requirement of a species.
Species distribution models (SDMs) can be used to gain ecological insights and to predict distributions across landscapes and to extrapolate in space and time [4].SDMs relate distribution data (e.g., occurrences at known locations) with information on the environmental conditions at these sites [4].On a continental scale, tree species distribution is mainly caused by climate variables [5] and only on a finer scale by soil or terrain.For example, Coudun et al. showed that soil nutritional factors improve distribution models of species with a high demand of alkaline cations in the soil solution like Acer campestre [6].Additionally, tree species distribution is influenced by management which can result in over-or underdispersion of some species.Reference [7], for example, reviews major effects on beech distribution in north-eastern Europe on the one hand due to conversion of beech stands into agricultural land and later into coniferous forests since the medieval epoch and on the other hand due to planting and promotion of beech near or even outside the reputed distribution range in Prussian times.The SDM approach assumes equilibrium between species distribution and environmental conditions.Problems might arise in some areas where a species is not present due to human impact or when the postglacial migration is not completed (underdispersion).If those areas have comparable environmental conditions as areas with presences, a model prediction for those areas is possible.
SDMs most often use presence-absence data, the model output is a probability of occurrence.This probability can be transformed into binomial data with the help of a threshold value (e.g., [8]).According to Real et al. [9], the probabilities do not only depend on the values of the predictor variables but also on the relative proportion of presences and absences.In order to overcome this effect and to describe the environmental favourability of a site for or against a species, Real et al. [9] suggest the use of a favourability function to receive results that are independent of an unbalanced proportion of presences and absences.Especially beech has more absences than presences in the data base; therefore the favourability transformation is used in this study.Norway spruce (Picea abies (L.) Karst.) and European beech (Fagus sylvatica L.), two of the most important tree species of the humid zone of Europe referring to forestry, are examined in this study.Spruce is considered to be strongly affected by the climate change, whereas beech is considered to be a possible suitable alternative as sites become dryer and warmer in central Europe [1].Both species are relatively unspecific to soil and terrain conditions [10,11]; therefore, their distribution is mainly determined by climatic constraints.Thus, our study concentrates on the relation of climate and distribution.Plant growth is limited amongst others by low temperatures (e.g., [12]) and, therefore, precipitation sum (sumVeg) and temperature sum (sumVeg) are calculated in respect of the start and end of the growing season.Beside these two variables, minimum temperature of the dormancy (minDorm) is used to describe species distribution.
Projected changes in temperature and precipitation patterns differ within seasons in Europe [13,14].Furthermore, the annual cycle changes, that is, the vegetation period is enlarging [15,16] and has to be represented by chosen climate indices (CI).Furthermore, plant growth is limited amongst others by low temperatures (e.g., [12]) and therefore physiological meaningful indices need to be related to a period of higher temperatures.If bud burst is earlier and leaf falls later in the year, the indices forcing the SDM have to take the start and end date of the vegetation period into account.Thus, forcing values are related to length of vegetation period in this study.The terms vegetation period and growing season are used synonymously and refer to a climatological definition.
When using climate data for SDM approaches, the question of appropriate climate data or models has to be solved.
There are a variety of data for present conditions or future climate projections [17][18][19], calculated with different models, bounding conditions, and assumptions.Many studies on species distribution refer to the WorldClim database [20,21], because these worldwide data have a very high spatial resolution [17].WorldClim data are interpolated climate surfaces for global land areas.Future data are calculated with a statistical downscaling of the output of global circulation models (GCM) and these anomalies are applied to the WorldClim data.In general, it is hard to judge which data set is the best [22].Therefore, the use of an ensemble of data seems to be the state-of-the-art solution to the question of the right data set.This is especially true for future climate data because there are many different data experiments with global circulation models (GCM) which in turn are forcing different dynamic regional circulation models (RCM).The GCMs are set up under defined assumptions, the so-called scenarios (Special Report on Emissions Scenarios, SRES).In this study, only the SRES A1B scenario [23] and one SDM technique has been considered in order to concentrate on differences of circulation models and downscaling methods and to show their effects on SDM output.The A1 storyline and family describe a future world of very rapid economic growth, global population that peaks in mid-century and declines thereafter, and the rapid introduction of new and more efficient technologies [23].The technological emphasis of A1B is a balance across all sources.That means that similar improvement rates apply to all energy supply and end use technologies [23].Cumulative emissions of greenhouse gases are in a broad range and in the middle of the different scenarios and families [23].Due to this, A1B is often used in climate impact studies (e.g., [24]).
The framework from assumptions of boundary conditions for climate change projection to final merged results of spatial favourability distribution for tree species is given by a chain of separate steps: For every chain member, several different models, approaches, or assumptions can be taken into account and enlarge the number of ensemble members.Due to the usage of data and model ensembles, a corresponding number of SDM results had to be merged.This final step can be taken as an additional chain member.
In this study, we focus on the methodological framework due to the chain of steps and present an approach focusing on ensembles of different GCMs and RCMs.Investigation area is Europe, distribution models are built by means of climate data from the period 1971-2000, and maps are calculated for this period and two future timesteps (2021-2050 and 2071-2100).The study is focusing on the following issues:

Data and Data Preprocessing
The presented SDM approach is driven by downscaled regional climate model data.The data are taken from the ENSEMBLE-Project [13,19] where data experiments with different RCMs, forced by different GCMs, are available.Within the ENSEMBLES-Project 22, data experiments have been realized with regional climate over Europe, all following the assumptions of the IPCC emission scenario A1B [23].In case of spruce we left out one realization and in case of beech two (marked in Table 1).The decision to exclude these datasets is based on the shape of the single response curves.Downloaded ENSEMBLE data were preprocessed in order to receive a uniform file format.To ensure a later on merging, all datasets were remapped to a unique domain from 15.125 ∘ W to 39.875 ∘ E and 35.125 ∘ N to 70.125 ∘ N and a grid of 0.25 ∘ resolution (according to the E-OBS data grid [25]).Beside these remappings, all units were harmonised for timesteps (in UNIX epoch, seconds after 01.01.1970), daily surface temperature (in Kelvin) and precipitation (in mm/day).All preprocessing of regional climate data were realized with Climate Data Operators (CDO) version 1.5.3 [26].All in all, we used 21 climate datasets and 3 time slices (see Table 1).Due to too many missing data in nine data experiments for the year 2100, this year was excluded completely for the appropriate realizations.
Beside data describing the climate conditions, observation data of the geographical distribution of present occurrence of tree species is used in this study.The data for beech and spruce were taken from the ICP Forests large-scale forest condition monitoring (Level I) from the period 1987 to 2007.In order to receive a somewhat regular grid, some plots were left out where the density of plots was considerably more than one plot per 225 km 2 (e.g., due to shift of inventory plots within a country or duplicates).This led to more than 7500 inventory plots in a resolution of approximately one plot per 225 km 2 across Europe [27].
The available information for crown condition within the monitoring data were transformed into presence and absence data as required by the SDM approach of this study.Inventory data often comprise fallacious absences [28,29], for example, due to impact by forest management for long time periods.In case of beech there is fallacious absence in southern Scandinavia because there beech is not surveyed within the monitoring program.This means that the data include absences where a species could live successfully but it is not present due to other reasons than unfavourable climate conditions.The term fallacious absence used by [28] refers to this theoretical concept and does not mean that data are improper in a technical way.As we aim to describe the climatic favourability, we came up against fallacious absences by using a vegetation map [30] that provides an overview of vegetation units and tree species composition.We attached the information about vegetation units to the Level I plots and transformed it into presence and absence data for beech and spruce.Finally, the two columns were merged.Therefore, a zero means no presence in both of the datasets, whereas a one means a presence in at least one of the datasets.In order to receive a clear signal of the altitudinal alpine distribution limit, we added 56 pseudoabsences above the alpine tree line.The resulting final data set consisted of 7596 plots covering most of Europe with a lower spatial coverage in the eastern European parts in the chosen spatial extent.
We aimed to find independent data for validation and therefore merged data from three different sources: first, the left out Level I data; second, Level II (ICP Forests Intensive Monitoring Programme, http://icp-forests.net/page/level-ii/)data, and, third, data from the European Forest Genetic Resources Programme (http://www.euforgen.org/,resp., http://www.eufgis.org/).All data were transformed into presence and absence data, no treatment for false absences was carried out and the use of the data from the Genetic Resources Programme was restricted to stands with more than one species in order to increase the probability that a full species inventory was behind the data.All in all, a maximum of 2135 stands could be used for validation.The data were unevenly distributed but covered most parts of Europe.Presences and absences of spruce fitted better to rough distribution maps (http://www2.biologie.unihalle.de/bot/agchorologie/areale/index.php?sprache=E/) than presences and absences of beech.Validation data were combined with climate datasets.Because of missing values, for example, close to the shore line in some combinations of GCM with RCM, we left out all those validation sites that did not contain the complete climate information of all 21 ensemble members.

Methods
The overall method is to relate species occurrence data to an ensemble of climate data and to transfer this relation to climate projections for this ensemble.The distribution modeling, quality assessment, and merging of the final maps is described in this section (Figure 1).The section refers to the last to chain links of the chain mentioned in the introduction.
The SDMs are forced by climate indices which refer to a vegetation period.These indices were calculated for each grid box.A grid box is the smallest spatial unit of the climate datasets which have a resolution of 25 km.As the vegetation period is the base for the indices, the start and end day in the year of this period has to be determined first.The vegetation period, or synonymously the thermal growing season, in this study is defined as a period with daily mean temperatures ≥5 ∘ C for spruce and ≥10 ∘ C for beech.We used a running average over 11 days for the calculation of the start and end day using the daily resolved temperature data.Temperatures had to be above or below the threshold for at least five consecutive Table 1: Ensemble members for distribution calculation and the appropriate GCMs and RCMs.Detailed description of the models is reported in [13].days.Furthermore, some limitations based on a day length criterion were set as default: the vegetation period was not allowed to start before first of March (60 day) and could not end before first of July (152 day) but had to end latest after first of November (305 day).
If five consecutive days with temperature above, respectively, below 5 ∘ (10 ∘ )C were detected, the first day of the consecutive timesteps is taken as the start respective end day in the year.The following equations are describing the dedications of vegetations periods begin as well as their ending: The length of the vegetation period (VP) can be defined as The following indices are used, which are derived from daily resolved temperature and precipitation values considering the information on the length of the vegetation period.
Precipitation sum in the vegetation period: sumVeg growing degree days in the vegetation period: sumVeg total minimum temperature outside the vegetation period: with  = daily precipitation amount (mm) The selection of physiological meaningful variables was done according to our hypotheses on the climate dominated distribution of beech and spruce on the continental scale and according to preliminary work (variable importance of a boosted regression tree analysis; see [4] for details on boosted regression trees).
minDorm is detected based on a smoothed temperature line to reduce misdetection of extreme values and stored for every year.The smooth is again the running average over 11 days.Finally, a mean over the respective 30 year time period was calculated for each index and for each grid box.The climate indices were calculated for each ensemble member and time period separately and joined with the species distribution data.
In order to understand the time dependent trends of the climate data, time series of the indices were calculated as a mean over a small test region for each climate data set.We chose the Bavarian Forest as a test region 12.625 ∘ to 13.875 ∘ E and 48.875 ∘ N to 49.125 ∘ N as this landscape is suitable for growing beech and spruce.
We used generalized additive models (GAMs, [31][32][33][34]) using the R-package MGVC [33] (R-project version 2.15.0 (2012-03-30)) to calibrate the single SDMs.GAMs are semiparametric extensions of generalized linear models.They allow fitting of response curves with a nonparametric smoothing function instead of parametric terms.This improves the exploration of species responses to environmental gradients.GAMs were used without interaction terms in order to build parsimonious models.An individual GAM was fitted to each climate ensemble member of the present day time slice  for beech and spruce.Response curves were calculated for the climate indices (sumVeg, sumVeg and minDorm) and checked for physiological plausibility according to [35].We used thin plate regression splines [34]; the degree of smoothness in the package mgcv is in general estimated by a generalized cross validation (GCV) criterion within certain limits set, for example, by the dimension  of the basis used to represent the smooth term.-values were first set to 3 and adapted (maximum  = 6) in some cases in order to meet assumptions on physiological plausibility.These values affect the maximum possible degree of freedom for each term of the GAM.One climate data set was excluded due to implausible precipitation data (no decrease of favourability at low precipitation for beech and spruce in Europe, that is, no distribution limit due to drought) and another one for beech (implausible response curve of variable growing degree days) after checking the single response curves for plausibility.All in all, a separate calibration of the GAM for each data set was done.
Different quality criteria were taken from the GAM output and additional parameters like the area under the receiver operating curve (AUC), sensitivity, and specificity were calculated for learning and validation data with the R package PresenceAbsence [36].Additionally, a tenfold cross validation was run, leaving out one-tenth of the data each time and using these data for validation.The GAMs were used for the calculation of distribution maps.These maps depict the probability of occurrence for Picea abies and Fagus sylvatica for each ensemble member and time slice (present day and future).In a further step, probabilities were transformed into favourabilities according to [9] in order to allow an easy comparison between the species.This is due to the fact that this transformation corrects the effect of the frequency of presences in the learning data (prevalence) on the probabilities [9].The output is a floating number between 0 (no) and 1 (maximum favourability) and describes the environmental favourability of a site for a species.
The last step was the merging of the separate results for each ensemble member to a final result.Two strategies were realized: an average over all appropriate ensemble results as well as a pseudodensity calculation.The average favourability (mean) of a grid box was only calculated if grid boxes contain valid values for all ensemble members.If a grid box contains a "Not a Number" (NaN) value, this 25 km rectangle was removed from the averaging procedure (e.g., at coast-lines).Additionally, the variance was calculated.The second way of integrating the different results was the calculation of a density map.Calculating density functions is a concept of merging ensembles with a large number of members.In this study we received 15 (far future), respectively, 21 (near future) members and therefore called the results due to the low number "pseudodensity function." Favourability output was converted into a binomial response (site is favourable or not) using a threshold of 0.5 (the prevalence transformed to the scale of favourability) according to [8].For each ensemble result, the favourability of ≥0.5 and ≤1 was taken as a one and otherwise as a zero.The pseudodensity is the proportion of ones at each grid cell.If all ensemble members have a one, then the ratio or pseudodensity is one; if half of the ensemble members have a one and the other half a zero, then the result is 0.5.The distribution shift was determined by subtracting the pseudodensity of the future predictions from that of the present time.

Results
Due to the chain from future scenario assumption to species distribution projection, several interim data were calculated.The following section presents results and figures that illustrate differences in the climate data ensembles and these differences-as a consequence-can be found again in the distribution models.Single response curves of the GAM approach try to illustrate this as well as averaged distribution maps of beech and spruce.The section ends with the two different averaging methods that are explained in Section 3 and shows the respective distribution shift for future time slices.

Climate Indices Ensemble.
The output of the climate model runs is afflicted with high uncertainties, but in this work the focus is on the handling of different climate model output and these uncertainties are not taken into account.Therefore, uncertainties of the final distribution maps within the ensemble approach first of all are due to differences in the climate input data.For example, precipitation data of the 21 ensemble members do have different maximum amounts of precipitation.They agree in high precipitation sums in Alpine regions but differ markedly in the amount of precipitation at the coastlines (not shown).Since the climate data differ, the derived indices differ as well.As an example, Figure 2 shows the length of the vegetation period calculated with a baseline of 10 ∘ C for a test region in Eastern Bavaria.While the ensemble mean shows a vegetation period of around 130-150 days in the present time slice, single ensemble members differ with a range of more than 80 days.This variation is decreasing in the future due to the reduction of datasets as only 15 datasets were available for the far future.Additionally, in the case of the 5 ∘ C-baseline the fixed limit of 244 days (vegetation period only between April and October) is reached in some years and reduces variation.As climate indices (precipitation sum and growing degree days) base on the length of the vegetation period, these indices inherit this variability.During the century, an extension of the vegetation period of 15 days is described by the climate data within the first half and of up to 35 days by the end of the century.Generally, growing season is longer for spruce as an evergreen conifer than for deciduous beech due to the different baselines (5 ∘ C, resp., 10 ∘ C) which results in an earlier start and later end day of vegetation period for spruce.Figure 3 shows a summary of the growing degree days calculated for spruce at the inventory plots with the data of the time slice 1971-2000.This sample illustrates the differences of the climate data ensemble members which in turn have consequences for species distribution modeling.The differences are even stronger for precipitation sum (not shown) because precipitation is harder to model than temperature.

GAM Ensemble.
The resulting single response curves of the GAMs of spruce are shown in Figure 4.Even though there is a general accordance of most models in the shape of the response curves, in some cases there are marked deviations due to the different climate data (e.g., Figure 3).The differences are pronounced where the data are scarce, for example, in the range of high precipitation sums which occur at high altitudes or at some coastlines.Additionally, response curves diverge due to different maxima, for example, of precipitation sum, and therefore there is a shift in the range of values.The marked differences of the input data lead to marked shifts of the response curves but the shape which describes the type of response seems to accord for many SDMs.Basically, GAMs agree in describing a distribution limit at mild winter temperature (3-5 ∘ C minDorm), high temperature sums in the vegetation period (>2000-2500 ∘ C sumVeg), and low precipitations sums (<200 mm sumVeg).The individual GAMs differ in describing the relationship of low temperature sums and high precipitation sums with the occurrences of spruce.In case of beech (not shown), response curves of minDorm are bell-shaped, sumVeg describes a distribution limit at high temperatures, and sumVeg shows a distribution limit due to drought. of favourability.Accordingly, the resulting maps represent averaged climate conditions for the three periods 1971-2000, 2021-2050, and 2071-2100.As we cannot present all results for each ensemble member, we show the final result as a union of the entire ensemble by averaging the single members of each time period (Figures 5 and 6).Additionally, we present the corresponding variability in the same figures.
The favourability map of beech shows a high accordance with the presences of the training data (white dots in the figure).This is underlined by high quality measures, for example, a mean AUC (area under the receiver operating curve; see [37]) of 0.90 with a minimum of 0.87 for the learning data.Using the prevalence of 0.32 as a threshold, the SDMs on average classify 90 percent of presences and 76 percent of absences correct which leads to 80 percent of the training data correctly classified.The Kappa value (see, e.g., [37]) is determined on average as 0.59 with a minimum of 0.53.Quality measures of the 10-fold cross-validation are in the same order of magnitude with the same mean and 2 to 9 per cent lower minimum values.For the validation data, the figures decrease: a mean AUC of 0.76, 85 percent of presences but only 53 percent of absences are on average correctly classified.The mean Kappa value is 0.35.Besides the present day distribution and quality of favourability maps, the figure characterizes the possible impact of climate change on the species distribution.Beech seems to be stable in middle Europe (e.g., Germany) until 2050 with only moderate shifts from west to north-east but the shift is anticipated to increase by the end of the century (calculation with 15 ensemble members).The changes mostly affect areas with low present day favourabilities of beech like western France.
Overall, the highest variability between the maps for beech is at the edges of distribution and it increases in the future especially in regions with a rising favourability in eastern and north-eastern Europe.Hotspots of variability can be found around the Alps, the Pyrenees as well as certain parts of the Carpathian arc, in the east Ukraine, and in a belt from St. Petersburg, Russia, around the gulf of Bothnia (northernmost arm of the Baltic Sea), via Sweden up to the Norwegian coastline.Generally, regions with a predicted favourability of 0.3 to 0.6 often also show a high variance of more than 0.06.The variance of the present day ensemble is low and thus points to a reliable result.
The situation for spruce is different (Figure 6).First of all, the distribution edge, for example, in Northern Germany or France appears to be not matched as good as in the case of beech (white dots with low favourability).These areas are highly affected by a strong northward shift of spruce's favourability by the end of the century.Second, the predicted distribution shift for spruce is larger than for beech and all resulting maps agree.That means that there is a low variance of the mean.Highest variances are found at the edge of the investigation area (a hotspot in Ukraine stretching into Belarus and Russia marking the southern limit of the distribution).Models have a mean AUC of 0.91 (minimum 0.89), a Kappa value of 0.67 (minimum 0.63), classify 83 percent correctly, a sensitivity of 0.89, and a specificity of 0.77 calculated with the training data and the prevalence as a threshold value (0.51).Again, 10-fold cross-validation results in the same order of magnitude with 1 to 6 per cent lower minimum values.Validation with the independent data provides the following quality measures: a mean AUC of 0.88, a mean Kappa value of 0.61, and 81 percent correctly classified.
4.4.Pseudodensity. Figure 7 presents the calculated pseudodensity functions and the differences between present day and future distribution of beech and spruce.There is a high accordance of the models describing the distribution expressed in a narrow transition belt from high density values to low ones and only a few differing GAM outcomes at the distribution edges (e.g., southeast Europe in case of beech).A few areas with low density in Italy or southwest France are backed as suitable by presences (white dots in Figure 7(a)).Interpreted as a measure of accordance or disaccordance, the pseudodensity helps to illustrate areas where the results are reliable due to high accordance (e.g., areas where 80 percent of the SDM agree).Additionally, the figures show the distribution shift due to climate change in terms of density shifts, allowing a clear picture of where favourable areas for growing beech and spruce decline (values below zero, red areas) and where new areas are gained (values above zero, green areas).Certainly, there is a high accordance between averages and pseudodensity maps.The results for beech as shown in Figures 7(a), 7(c), and 7(e) indicate a high accordance of different maps, in the middle of the century a slow north-eastern distribution shift that accelerates by the end of the century.In 2021-2050 beech stands at the western (western France) and southern distribution limit (Apennine Mountains in Italy, Hungary and some parts of southeastern Europe) will be affected first.These areas are widening by the end of the century also affecting the south of England, Ireland, and northern and northeastern Germany.On the other hand, new favourable regions for growing beech are gained in higher altitudes and northeastern Europe.In the case of spruce (Figures 7(b), 7(d), and 7(f)), the main spatial distribution in terms of pseudodensity function is situated around the Baltic Sea.The distribution includes the whole Scandinavian Peninsula except the higher altitudes of the Norwegian Mountains.The southern distribution border is characterized by a relatively large transition zone of lower model agreement along 50 ∘ N for the more continental regions and shifting northwards within an oceanic influence.Further south, the agreement of the entire ensemble is large again in the higher altitudes of the Carpathian Arc, the Swiss Alps, and in a few regions of the Pyrenees and the Balkan regions.Beside smaller parts of the Massif Central in France, the higher altitudes of the German low mountain ranges are also pointed out as suitable sites for spruce under present climatic conditions described with the SDMs.Nevertheless, the presences (white dots in Figure 7(b)) in northern Germany, the Benelux countries, and in France emphasizes the problems of many spruce stands in warmer climates today.They are in climatically unfavourable areas outside their natural distribution range and further warming as projected by scenarios indicates a high risk of decline.Sites with a raising favourability for spruce in the future can be found in the Scandinavian mountains and at the actual northern distribution edge.The distribution shift also shows that favourability for spruce declines in Germany (lower mountain ranges and alpine upland), Czech Republic, northeast Poland, Pyrenees and French Massif Central, and southern Sweden.

Discussion
This section provides a discussion of the major steps of the presented modeling chain including the climate data, the indices driving the SDMs, the SDM approach, and the resulting final maps.

Climate Data and Ensemble
Modeling.Using projection ensembles in SDMs is a common technique in order to reduce variance and uncertainty and to raise the accuracy of species distribution models [38][39][40].Whereas many studies focus on ensembles of SDM model classes, our focus is on the climate data, an approach becoming more important in the SDM literature [41].Perhaps due to limited access to climate data or problems in handling huge data volumes, ensembles often are limited to a couple of combinations of GCM and RCM (e.g., [1], but see [24]).Beside the assumptions made on the development of greenhouse gas emissions, the realizations of climate data via GCM and RCM differ depending on the particular GCM and RCM used.This introduces an additional range of uncertainty into the model chain.Regionalized downscaled climate data differ markedly, even when using only one emission scenario (Figures 2  and 3).We therefore concentrate in this study only on the different realizations of GCMs and RCMs and not on different emission scenarios.
A state-of-the-art ensemble model set is that of the ENSEMBLES Project.Lorenz and Jacob [42] compared the RCMs driven with quasiobserved lateral boundary conditions with gridded observational data [25] for near surface temperatures and concluded an underestimation of the past temperature trends as well as a wide range of the results of the RCMs.Jacob et al. [14] reported an uncertainty of several degrees as well as uncertainties in precipitation among the projected model data.The uncertainties are increasing during the century up to Δ 2.5 K for running averages and Δ 30% in case of precipitation.The changes vary during the seasons and different regions.The difficulty of judging the single members of the ensemble leads to an ensemble without weighting their members [22,43].
Not all ensembles members have a realization over the whole century; in seven cases the future projections were only resized to the year 2050.Therefore, a smaller variance at the end of the century does not mean a lower uncertainty.

Climate Indices:
Use of Daily Climate Data.We focused on three climate indices describing the environmental conditions driving species distribution.An alternative for using temperature and precipitation is the use of aridity indices (e.g., [44]).They have the advantage of combining two variables and therefore taking interactions into account.The disadvantage is that they are normally developed for climate data with a monthly or yearly resolution and that they are therefore weaker related to plant physiology.Two of the indices used in this work-temperature and precipitation in the growing season, sumVeg and sumVeg-depend on the length of this period, which can be estimated with daily climate data.Daily resolved data are suitable for the tackled questions but require methods of handling large data amounts.
As shown in Figure 2, the vegetation period is slightly enlarging with rising temperatures in the future projections which is in accordance with extrapolations made from phenological observations [16,45,46].The projected temperature rise for the 21st century has a twofold effect on the calculated indices.On the one hand, there are higher temperatures leading to a higher temperature sum.On the other hand, the reference base for the calculation of the temperature and precipitation sum-the length of the vegetation periodis enlarging too.This can intensify a trend in case of the rising temperature sums as well as compensate decreasing precipitation sum values due to more days which are taken into account to the end of the century.
We use the vegetation period for two of the selected indices because precipitation and temperature during the physiological active phase are most important for water supply and evapotranspiration and therefore tree growth [12] and thus for species distribution (e.g., see [47], a study on drought and mortality of spruce).Due to the large climatic differences within Europe, a fixed time period presumably has a weaker relation with tree growth.We use a simple but straightforward approach to define the growing season (mean temperature threshold and day length criterion).The distribution of beech and spruce are located in climate regions with thermal seasons.The detection defined by five consecutive days higher respective lower than a threshold based on a smoothed temperature line is a very simple but useful approach.There are more sophisticated ones (e.g., [45,46,48]), but due to the variances within the ensemble members there is no benefit expected with a more complex approach for vegetation period detection.
The chosen simple definition of the vegetation period may not be suitable for all European regions, that is, in southern Europe where the boundaries 1st of March and 31st of October are reached too often, due to their default exclusion.This happens especially in the case of spruce for which we assumed a temperature limit of 5 ∘ C for considerable transpiration and other physiological growth processes (Raspe pers.com., [45,[48][49][50]).We think our approach still does make sense because the distribution models are developed for species of the temperate to boreal zone and regions with high temperatures almost all over the year are far beyond the distribution range of spruce and beech.An important point is that we opened the door for using regionalized daily climate data.In future works more sophisticated indices concentrating on drought events in the growing season or drought indices can be implemented.

Use of Vegetation Map Data for the Treatment of Fallacious
Absence.Species distribution models are only as good as the calibration data.Hence, if environmental data do not explain the distribution well because of unaccounted effects like forest management practice (e.g., preference of one species over another due to cultivation traditions or economic reasons), postglacial migration borders due to competition, or land use effects, models may have high prediction errors.As forest inventory data are affected by human impact, we cannot assume an equilibrium between species traits and distribution leading to errors due to fallacious absences [28,29].That means that a species is absent although the site is favourable.The use of expert knowledge is one possibility to treat these fallacious absences (see [44]).So on the one hand, the used vegetation map helps to improve the distribution data.On the other hand, it might introduce some uncertainty due to its scale.In our opinion, the advantages outweigh the disadvantages because the distribution maps now give a general picture where spruce and beech potentially can grow.Whether the data describe the full potential of a species is not clear and it is obvious that our models are influenced by the human impact on species distributions.In case of spruce, we assume that plantations show the distribution limit in continental Europe.Ireland might be an example where our data is not complete and do not properly represent the potential distribution: models describe low to medium favourability for Ireland despite presences in the learning data.Absences dominates over the few presences in this area.The maps show the potential distribution area [51].This information is valuable for forest management decisions dealing with adaptation to climate change.Whether species will colonize this potential distribution depends on many factors like, for example, migration speed and barriers, dispersal and establishment limitations, and pests that, for example, can be described with mechanistic models [52].
The independent validation data were not corrected for fallacious absences and therefore reflect the realized distribution and show the differences between potential and realized distribution.Spruce is a tree species that is planted outside its natural distribution range with the help of forest management.Therefore, there are not so much fallacious absences and realized distribution is almost equal to potential distribution.Hence, validation results are in the same magnitude as cross-validation results (AUC, Kappa, percent correctly classified stands).The situation is different in case of beech, where quality measures of the validation with independent data are lower compared to cross-validation results.Beech is missing at sites where vegetation experts judge it to be within its natural distribution range (fallacious absences).It might be missing due to forest management and preference of spruce or oaks, for example, for economic reasons.Fallacious absences lead to higher error rate (percent correctly classified absences, Kappa value) when comparing these absences with a model that describes the potential distribution.That means that there could be a gap between actually realized distribution and potential distribution as depicted by vegetation experts because the cross-validation as the second validation method backs the model quality.Independent validation data could have been corrected for fallacious absences but then there would not be obtained more information of this method than by the cross-validation procedure.

Use of Climate Data as
Proxies for Water Supply.We used climate data to describe the distribution of spruce and beech.Soil influences species distribution, too (e.g., preference of sites with high base saturation in case of Acer campestre [6]).However, spruce and beech are highly competitive, widely distributed in Europe, and do not have a special preference to soil conditions (e.g., nutrients) or certain relief positions [10,11].Soil is also part of the hydrological cycle and therefore does matter when drought events are the reason for distribution limits.Generally, it is difficult to find adequate European soil data and the influence of soil physics on water availability on a European scale is small compared to the effects of precipitation and temperature.On a local scale it is much more important to incorporate soil data into SDMs.In case of the Level I data, the European demonstration project BioSoil might fill the gap of missing soil data in the future [53].Nevertheless, climate data as proxies for water supply in the growing season appear to be sufficient to reproduce the general distribution pattern of most tree species in Europe and to risk a glimpse into what climate change might mean for the distribution of these species.The picture will change and distribution models improve when taking into account frequency, intensity, and duration of drought [54].Daily climate model data prepare the technical ground to focus on extreme events.

Species Distribution Modeling Technique-GAM.
Following the argumentation of Wang et al. [41] and referring to the work of Marmion et al. [39], we concentrate on one SDM technique.Reference [39] recommends either a single best method or a consensus approach of different methods like realized within the Biomod framework [40].
Other than [41], we do not use a machine-learning technique like Random Forest or a Boosted Regression Tree technique despite their good performance [4,52,55], but use GAM, a well-established semiparametric regression technique [55].GAM can be handled to model smoothed continuous gradients meeting the continuum theory [56].Additionally, it can be run in a parsimonious way using only a few dose variables, whereas the strength of data driven machine learning techniques amongst others is to deal with a high number of variables and their interactions [4].It is especially important to avoid overfitting when balancing between model fit and predictive performance because overfitting may reduce transferability [4,57].We calibrated the GAMs in order to give a generalized picture of how the distribution of spruce and beech is determined by precipitation and temperature.Therefore, GAMs were run with low -values (basis dimension of the penalized regression smoother,  − 1 is the upper limit of the degrees of freedom and influences computational efficiency) which leads to a high degree of smoothness and therefore results in very general trends (see partial response curves in Figure 4).Thus, we avoid overfitting and the SDMs are as reliable as possible in terms of transferability or more precise forecasting [38,44].

Calibration of GAMs and Plausibility of Individual
Climate Data Experiments.One step in the presented species distribution modeling was the individual tuning of the single response curves in order to receive physiological plausible models according to some fundamental hypotheses about the relationship between climate and species distribution.In doing this, we follow [35].Elith et al. [58] also suggested to deliberately control the fit of models, for example, by using exploratory techniques like checking response curves for their contribution to the final outcome.They even went one step further and successfully integrated information from mechanistic models in order to enhance the reliability of correlative predictions.When calibrating the GAMs individually for each climate data set, we also tried to reach a high sensitivity (correctly modeled presences) because this is a requirement for describing the ecological potential of a species.Nevertheless, response curves differ which is a result of the different climate data and highlights the uncertainty of models in some areas, for example, at the edges of distribution where presences are scarce.As we cannot judge the quality of the different datasets, we use the ensemble technique in order to obtain a more complete picture of the distribution and a possible change with rising temperatures by the end of the century.A comparison of differences in input data (Figure 3) and SDMs (Figure 4) reveals the dependence of the SDMs on the climate data and can be used to cautiously judge the quality of input data for species distribution modeling: if there is no drought limit for the distribution of beech and spruce within a European data base, then this data set seems to be inappropriate, at least for a species distribution modeling approach.Overall, lower levels of temperature or precipitation leads to leftwards shift of single response curves in Figure 4. Thus, the SDMs do not differ markedly but the climate indices.

Transformation of Probabilities into Favourabilities.
We transformed the probability of occurrence as provided by a GAM into favourabilities according to [9] in order to compare the two species with different presence/absence ratios (prevalence) easily.According to [9], "favourability values can be regarded as a degree of membership of the fuzzy set of sites whose environmental conditions are favourable to the species." As we want to describe the potential distribution area of spruce and beech, favourabilities seem to be an adequate quantity to show where spruce and beech can grow today and in the future.

Effect of Climate Change on Species Distribution in
Europe.On a European scale it is obvious that spruce as a mountainous species is more affected by climate change than beech.This means that there are many spruce stands that will have a high risk of mortality due to warmer and drier stand conditions.There are different ways of managing risk in managed forests like shortening of rotation period, thinning [59], and in the long term transformation into stable mixed stands.We work on the species level not taking provenances into account due to the European scale and the fact that information on provenances are scarce.It can be expected that on a stand level, changes in vitality will be even more dramatic and adaptation to drier conditions will not keep pace with the projected warming.Therefore, our attempt treating the species as a whole and leaving genetic differences out might lead to an optimistic picture for single stands and a pessimistic one for long term adaptation processes.The latter is because we model the whole distribution range of the species and do not take a closer look on the distribution edges.With the help of forest management, today spruce is already grown beyond its natural distribution range at relatively warm sites.This means that stands that will fall outside the range in the future might survive for a certain time but the mortality risk is increasing strongly.Mortality might not be caused directly by drought but drought leads to lower vitality and subsequent pest attack.It is clearly visible that beech can be an alternative species on stands that will be too dry for spruce in the future because drought is a major discriminating factor at the warm and dry distribution limit between favourabilities for beech and spruce.Beech can withstand drier conditions because of deeper rooting and the chance to reduce leaf area in extreme situations in summer by a premature fall of leaves.Nevertheless, the distribution of beech has a drought limit too.With lower competition (e.g., growth) due to drought [60] beech fails to shade out other species.

Final Maps: Mean and Density.
Using ensemble modeling techniques entails handling of the multiple output generated.A simple but straightforward approach is to calculate a mean and the variance if there is no evidence of differences in quality of the resulting models or maps [22].Otherwise, a weighted mean (or median) should be favored [22].We present this approach (simple mean) as it leads to a map describing the whole range between suitable and unsuitable environments without losing information due to classification.Typically, species distribution models that are built on binomial distribution data lead to binomial output maps with the help of a threshold value.This leads to a clear picture of the final maps but has the disadvantage of losing information.A compromise is the use of multiple thresholds [44], linking the two approaches.We present another way of handling the ensemble output using one threshold, the prevalence [8], and thus combining binomial maps.we call it a "pseudodensity function" as the amount of ensemble members is limited to 20, respectively, 21, or 15, in case of the far future.The maps of the two approaches (mean and pseudodensity function) look similar to a certain degree but the density map has sharper distribution limits due to its binomial input maps.

Conclusion
Projections of the translocation of climatic favourability for beech and spruce can be done with SDMs forced by climate model and species presence/absence data because the growth of these tree species is mainly depending on climate conditions.Thus, the presented approach can be used for tree species favourability projections where edaphic factors can be neglected on a European scale.
The study shows consequences that can be expected for spruce and beech due to climate change with regard to their distribution in Europe.The general trend is a north-eastward and altitudinal shift of climatological favourability for both species, though stronger in case of spruce.These outcomes are in accordance with other studies [1,24] and can be used as a basis for risk assessment.The altitudinal shift is only relevant for beech because spruce already is part of the Alpine tree line and the area above this line with sufficient soil substrate for tree growth is limited.In general, the distribution shift will strongly affect the middle European spruce stands and underlines the need for stand conversion of pure spruce stands into mixed forests.Beech can be part of mixed stands and therefore used to stabilize spruce stands but certain beech stands at the warm and dry edge of distribution might need a conversion too.The presented work gives an example on how gridded daily resolved values of climate parameters, like temperature at 2 m and surface precipitation, can be used to describe the climatic site conditions for tree species.With these daily resolved climate model data very specific indices can be calculated.We present a rather simple approach to determine the vegetation period just to give an idea of the potential of temporal highly resolved data experiments.
The use of an ensemble with 21 members reduces the risk of a wrong projection.The study shows the importance and advantages of ensemble projections by pointing out regional differences of the uncertainties for each time slice.The outcome of an impact model like a SDM with climate data as input depends strongly on these differences.Beside a wide spread and common use of an arithmetic average combined with the appropriate variance, pseudodensity functions are another tool to visualize uncertainty originating from the climatic data ensemble.The individual check for physiological plausibility reduces bias and uncertainty deriving from the SDM framework, though this source of uncertainty cannot be eliminated completely although the used GAM technique, separately trained for each ensemble member, minimizes the individual bias of the different climate model members within the ensemble.

Outlook
The data handling with daily resolved climate model data and the connection to species distribution modeling is done with the presented work.Now that the frame is set, extension into different ways can be made: the focus can be set on variability and extreme events [54] and their influence on species distribution.For example, drought indices or deviations from long term mean might provide a better explanation of species distribution than simple mean values [54].Certainly, the models could be run for other species, or other SDM techniques.
The Coordinated Regional Downscaling Research Experiment (CORDEX http://wcrp-cordex.ipsl.jussieu.fr/) is the successor of ENSEMBLE-Project and is forced by the RCP scenarios [61] instead of those of the SRES scenario families used in the ENSEMBLE-Project.If these data are available, an ensemble with more members, higher spatial and temporal resolution, and new emission scenarios can be taken to force impact models like SDMs.Amongst other things, the CORDEX data will reduce the data preprocessing due to observing some standards.Future work will be done with these new climate data.
To reduce the workload or invite impact modelers with limited access to climate data, the shown GAM approach will be available in a virtual computing environment (C3Grid-INAD) and therefore can help to process large datasets under the control of the user (e.g., use other climate indices for sensitivity or extreme events studies or different definitions of the vegetation period).

( 1 )
shift of climatic favourability for European beech and Norway spruce in Europe; (2) incorporation of ensembles of gridded and daily resolved climate model variables into a SDM framework; (3) accentuation of uncertainties of the final results that derive from the use of an ensemble of climate models.

Figure 1 :
Figure 1: General framework of the separate methodical steps to model species distribution based on an ensemble of climate models with 21 respective 15 members.
with day = day in the year   = start of vegetation period (day)   = end of vegetation period (day) run = running mean with 11 timesteps of daily mean temperature.

Figure 2 :
Figure 2: Time series of length of vegetation period with the baseline of 10 ∘ C for the test region Bavarian Forest for all ensemble members (grey lines).Red line shows the mean.

4. 3 .Figure 3 :
Figure 3: Comparison of the climate index growing degree days ( ∘ C) at the inventory plots.The ensemble of climate data is shown for the period 1971-2000.Growing degree days are calculated depending on the length of the vegetation period defined as a period above 5 ∘ C. Base line of the temperature sum is 5 ∘ C (input data for modeling the distribution of Picea abies).

Figure 4 :
Figure 4: Single response curves of the GAMs for all climate data and species Picea abies.Single response curves are calculated for each climate index (minDorm, sumVeg, and sumVeg) while keeping the two others constant at the mean value.Single response curves give a hint on model behavior and physiological plausibility.Legend names are combinations of RCM and GCM.

Figure 5 :
Figure 5: Mean ((a), (c), and (e)) and variance ((b), (d), and (f)) of the spatial distribution of Fagus sylvatica (favourability according to [9]).White dots in (a) are beech occurrences in the learning data.Mean and variance of (a), (b), (c), and (d) underlie an ensemble of 20 climate data, whereas the ensemble of (e) and (f) consist only of 15 members.

Figure 6 :
Figure 6: Mean ((a), (c), and (e)) and variance ((b), (d), and (f)) of the spatial distribution of Picea abies (favourability according to [9]).White dots in (a) are spruce occurrences in the learning data.Mean and variance of (a), (b), (c), and (d) underlie an ensemble of 21 climate data, whereas the ensemble of (e) and (f) consist only of 15 members.

Figure 7 :
Figure 7: Shift of distribution of Fagus sylvatica ((a), (c), and (e)) and Picea abies ((b), (d), and (f)) expressed as difference of density distribution of the SDM ensemble in reference to time slice 1971-2000.White dots in (a) and (b) are beech respective spruce occurrences in the learning data.Density distribution in (c) and (d) underlies an ensemble of 20 respective 21 members, whereas the ensemble of (e) and (f) consists of only 15 members.
The experiments cover the time period from 1951 to 2100.Seven datasets only completed up to the year 2050; thus only 15 of the 22 realizations are available for the far future.The combination of RCM and appropriate forcing GCM is shown in Table1.For all ensemble members (single datasets), daily resolved values of near surface air temperature (at 2 m) and surface precipitation were taken.As indicated in the table, we excluded ensemble members, because of implausible values at least in some European regions (e.g., very high precipitation over Spain) which led to physiological implausible SDMs.