1 Introduction

The European Alps (hereafter referred to as “the Alps”) have warmed twice as much as the global average between the late 19th and the early twenty-first century (Auer et al. 2007; Ceppi et al. 2012), and a trend of faster warming at higher mountain elevations is predicted to continue until the mid-twenty-first century regardless of the choice of future climate scenario (IPCC 2018). Beyond the mid-twenty-first century, atmospheric warming in mountains will be stronger under a high greenhouse gas emission scenario (RCP8.5) and will stabilise at mid-twenty-first century levels under a low greenhouse gas emission scenario (RCP2.6) (IPCC 2018). While temperature trends in the Alps are clear, precipitation trends in the Alps show high spatial variability (Auer and Böhm 1994) with precipitation increasing in the NE and decreasing in the SW part of the Alps over the past decades (Brunetti et al. 2006; Auer et al. 2007). Taking into account the entire Alps, only slight shifts in average precipitation (Auer et al. 2007) and decrease in snowfall were recorded over the past decades (IPCC 2014). Summer precipitation reduction associated with increased convective rainfall is projected (Giorgi et al. 2016), whereas future trends in annual precipitation are uncertain (IPCC 2014).

The 4395 glaciers currently still present in the Alps with a total area of 1806 ± 60 km2 (Paul et al. 2020) are projected to lose a substantial part of their volume in the coming century and many of them are expected to disappear by 2100 regardless of which emission scenario is used (Haeberli and Hoelzle 1995; Zemp et al. 2006, 2007). Current regional and global estimates of the future evolution of glaciers in the Alps are based on models of variable complexity and rely on different climate projections. These models generally suggest a glacier volume loss of about 65–80% between the early twenty-first century and 2100 under a low greenhouse gas emission scenario (RCP2.6), 80–90% under a moderate greenhouse gas emission scenario (RCP4.5) and an almost complete disappearance of glaciers (volume loss 90−98%) under a high greenhouse gas emission scenario (RCP8.5) (Marzeion et al. 2012; Radić et al. 2014; Huss and Hock 2015; Zekollari et al. 2019).

Although the contribution of glaciers in the Alps to sea-level rise is low and is projected to account for less than 0.3 mm sea-level equivalent for the period 2010−2100 (Huss and Hock 2015), their relatively large volume loss could have a profound impact on river runoff, as already suggested in some earlier publications from the twentieth century (Collins 1987; Chen and Ohmura 1990a, b). Several studies have reported an increase in average winter runoff over the past decades in glacier-dominated Alpine basins and a decrease in summer runoff (Bosshard et al. 2013; Bocchiola 2014; Addor et al. 2014; Bard et al. 2015), with peak water discharge from the majority of Alpine glaciers already surpassed (Huss and Fischer 2016; Huss and Hock 2018). Similar trends in the amount and seasonality of river runoff are projected to continue over the twenty-first century and are expected to affect downstream water management, related hazards and ecosystems. There are also socio-economic risks related to rapid changes in mountain glaciers, such as the decline of tourism (Spandre et al. 2019; Steiger et al. 2019), deterioration of drinking water quality (Hodson 2014), impact on livelihoods (Haeberli and Whiteman 2015; Carrivick and Tweed 2016) and hydro-electric power generation (Farinotti et al. 2019).

One way to quantify the response of glaciers to changing climate is by monitoring their equilibrium-line altitude (ELA) (e.g. Lie et al. 2003; Zemp et al. 2007; McGrath et al. 2017). ELA is defined as the spatially-averaged altitude of the set of points on the surface of the glacier where the climatic mass balance is zero at a given point in time (Cogley et al. 2011). Although the ELA is generally understood as annual ELA, many different variations and approximations of the ELA exist in the literature (e.g. Lie et al. 2003; Zemp et al. 2007; Braithwaite and Raper 2009). Most of the previous studies on ELA report either the steady state equilibrium-line altitude (ELA0) of glaciers, calculated from direct glaciological mass balance measurements or the local topographic equilibrium-line altitude, also referred to as the effective ELA (hereafter referred to as “effELA”), which accounts for the effect of topography (Anderson et al. 2018). In contrast, the environmental ELA (hereafter referred to as “envELA”) is the regional altitude of zero glacier mass balance without the effects of shading, avalanching, snow drifting, glacier geometry or debris-cover (Anderson et al. 2018), averaged over a period that is long enough to be climatic, usually several decades. The envELA represents the theoretical altitude at which a glacier can form.

In this paper, we adopt the terminology of Anderson et al. (2018), who distinguished between the envELA and effELA. We investigate the past envELA over the entire Alps and project its evolution under different future climate forcing scenarios using a simple parameterization of the ELA in terms of summer temperature and annual precipitation initially proposed by Ahlmann (1924) and recently updated by Ohmura and Boettcher (2018) based on the theory of mass and energy conservation. Our approach has the advantage that it does not require detailed input data from various sources as in more sophisticated models, nor the parameterization of ice flow dynamics driven by the mass balance (e.g. Huss and Hock 2015; Zekollari et al. 2020).

We use 200 years of climate records and forecasts over the period 1901−2100 for the distributed modelling of the envELA at 5 arcmin grid cell resolution over the entire Alps. Historical envELAs are reconstructed based on the Long-term Alpine Precipitation reconstruction (LAPrec) dataset (Isotta et al. 2019) and temperature datasets from the Historical Instrumental climatological Surface Time series of the greater Alpine region (HISTALP) (Auer et al. 2007). Previous simulations of ELA changes over the Alps were mostly based on earlier greenhouse gas emission scenarios (e.g. Zemp et al. 2007) or relied on global climate series with low spatial resolution (e.g. Huss and Hock 2015). In contrast, our model is forced with high-resolution CORDEX regional climate model (RCM) simulations for the European domain (EURO-CORDEX, Jacob et al. 2014), which provide a better description of local climatic conditions. The envELA patterns are validated against existing end-of-mass-balance-year ELAs from the WGMS FoG database over the entire Alps for the period 1948−2017 (WGMS 2019). The projected envELA pattern under various climate scenarios is additionally validated for the period 2006−2019 using annual highest end-of-summer snow line altitudes (SLAs) derived from Landsat data for a subset area of the Alps using a previously developed semi-automated remote sensing method (Racoviteanu et al. 2019).

Our goal is to use modelled envELAs as a proxy for climate change in the Alps. Specifically, our work is motivated by the following research questions:

  • How many glaciers in the Alps currently lie below the envELA?

  • What is the projected envELA rise in the Alps under different climate forcing scenarios?

  • How many glaciers in the Alps are likely to disappear (when envELA is above their maximum altitude) by the end of the twenty-first century under different climate forcing scenarios?

  • What is the spatial pattern of the envELA across the Alps?

2 Data and methods

2.1 Climate data

Our study area is defined between the latitude ~ 43° and 48° N and longitude ~ 4° to 17° E (Fig. 1) and covers the entire currently glacierised area of the Alps. We selected the long-term Alpine precipitation reconstruction (LAPrec) dataset (Auer et al. 2007; Isotta et al. 2014, 2019) and temperature historical instrumental climatological surface time series of the greater Alpine region (HISTALP) (Auer et al. 2007) (Table 1) to reconstruct the envELA back to 1901. We used EURO-CORDEX simulations (Jacob et al. 2014; https://www.euro-cordex.net/) to produce projections of the future envELA evolution by 2100.

Fig. 1
figure 1

Study area defined as the Alpine Convention perimeter (https://www.alpconv.org/en/home/). Glacierised areas in ~ 2003 are from Randolph Glacier Inventory (RGI v6.0) (RGI Consortium 2017). Grid cells at 5 arcmin spatial resolution cover locations of every glacier polygon centroid and represent a basis for the Alpine environmental ELA distributed modelling. Rectangle shows the subset extent used for remote sensing ELAs

Table 1 The historical gridded climatology datasets used in this study

LAPrec is a new gridded monthly precipitation dataset of the Alps on a regular 5 km grid that is spatially referenced in ETRS89/ETRS-LAEA (European Terrestrial Reference System 1989-Lambert Azimuthal Equal Area). LAPrec integrates information from two data sources: the HISTALP station data archive (Auer et al. 2007) which extends over the full period of interest, and the high-resolution station-based spatial analysis EURO4M-APGD (Isotta et al. 2014), which spans the period 1971−2008 and incorporates more than 5500 rain-gauge measurements on average per day. This is one of the densest in-situ observation networks over a high-alpine topography worldwide. Two variants of LAPrec exist. The longer version starts in 1871 and uses 85 almost continuous series, whereas the shorter version starts in 1901 and uses 164 almost continuous series. Both versions end in 2017. We decided to use the LAPrec1901 version due to its higher reliability (Isotta et al. 2019).

The original measurements used for LAPrec1901 have not been corrected for precipitation gauge undercatch and as a result are expected to generally underestimate precipitation, especially in winter months and at high elevations (e.g. Sevruk 1985; Adam and Lettenmaier 2003). In order to account for this effect, we added a gauge correction to the LAPrec1901 dataset by first estimating a mean monthly undercatch correction from the global precipitation data set of Legates and Willmott (1990) (UDel Version 3.01) (Table 1), spanning the period 1900 to 2010, and then using it to correct the monthly precipitation values in the LAPrec1901.

The HISTALP gridded dataset comprises monthly homogenized records of mean air temperature for topography height on a regular grid with a spatial resolution of 5 arcmin. The dataset was created by merging high resolution gridded climate means (1961−1990) for each month (Hiebl et al. 2009) and long-term monthly data from 140 stations (Chimani et al. 2013). The time period covered is 1780−2014.

For the climate model simulations, twelve sets of driving GCM and corresponding RCM experiments carried out in the frame of EURO-CORDEX at 0.11° resolution and three different emission scenarios for the variables near surface air temperature (tas) and precipitation (pr) are considered (Table 2) (Jacob et al. 2014; https://www.euro-cordex.net/). EURO-CORDEX RCM simulations are forced with different scenarios, denoted by Representative Concentration Pathways (RCPs) (van Vuuren et al. 2011; IPCC 2014). RCPs are identified by their approximate total radiative forcing in year 2100 relative to 1750. The three considered RCPs in this study include the RCP2.6 low end scenario leading to a very low forcing level (2.6 W m−2), the RCP4.5 stabilization scenario (4.5 W m−2), and the RCP8.5 high end scenario with maximum greenhouse gas emissions (8.5 W m−2). For RCP2.6, the radiative forcing peaks at 3 W m−2 before 2100 and then declines to 2.6 W m−2 by 2100; for RCP4.5 it stabilizes by 2100, and for RCP8.5, it does not peak by year 2100. Historical runs for the selected models were used for systematic bias correction. We compared the monthly precipitation totals and mean monthly temperature of the observed and simulated datasets over the period 1976−2005 (Fig. S1) and estimated a bias correction following the procedure described in Huss and Hock (2015). The bias is then assumed to remain constant in time also during the projection period and is used to bias-correct temperature and precipitation on a monthly basis from 2006 to 2100. For precipitation, we use a linear scaling bias correction method in which, for each month m of a year y of the projection period, precipitation is corrected as:

$$ P_{m,y, corrected} = \frac{{\overline{{P_{{m, LAPrec1901 \left( {1976 - 2005} \right)}} }} }}{{\overline{{P_{{m, EURO - CORDEX \left( {1976 - 2005} \right)}} }} }} \quad P_{{m,y EURO - CORDEX\quad \left( {2006 - 2100} \right)}} $$
(1)
Table 2 EURO-CORDEX RCM simulations at a 0.11° resolution examined in this study

where \(\overline{{P}_{m, LAPrec1901 (1976-2005)}}\) is the average monthly precipitation total from the LAPrec1901 dataset over the period 1976−2005, \({P}_{m,EURO-CORDEX (1976-2005)}\) is the average monthly precipitation total for the EURO-CORDEX historical simulation dataset over the period 1976−2005, and \({P}_{m,y EURO-CORDEX (2006-2100)}\) is the monthly precipitation total for every year of the EURO-CORDEX RCM simulation for the period 2006−2100.

Variance scaling after Huss and Hock (2015) was used to bias-correct temperature by which, for each month m in a year y, we apply the following equation:

$${T}_{m,y, corrected}= \overline{{T}_{m, 25}}+\left({T}_{m,y EURO-CORDEX (2006-2100)}- \overline{{T}_{m, 25}}\right) \frac{{\sigma HISTALP}_{,m (1976-2005)}}{{\sigma }_{EURO-CORDEX,m (1976-2005)}}+ \left(\overline{{T}_{m HISTALP \left(1976-2005\right)}}- {\overline{{T}_{m EURO-CORDEX (1976-2005)}}}\right)$$
(2)

where \(\overline{{T}_{m, 25}}\) is the centred running mean of monthly temperature in a 25-year period for y, \({T}_{m,y EURO-CORDEX (2006-2100)}\) is the mean monthly temperature for every year of EURO-CORDEX RCM simulation for the period 2006−2100, \({\sigma }_{HISTALP,m (1976-2005)}\) is the standard deviation of monthly temperatures of the HISTALP dataset over the period 1976−2005, \({\sigma }_{EURO-CORDEX,m (1976-2005)}\) is the standard deviation of monthly temperatures of EURO-CORDEX historical simulations over the period 1976−2005, \(\overline{{T}_{m HISTALP \left(1976-2005\right)}}\) is the monthly average of HISTALP over the period 1976−2005 and \({T}_{m EURO-CORDEX (1976-2005)}\) is the monthly average of EURO-CORDEX historical simulation over the period 1976−2005. In this way, the EURO-CORDEX temperature series are adjusted to account for differences in year-to-year variability between the historical (HISTALP) and projected (EURO-CORDEX) time series.

Because the historical and simulated datasets do not use the same horizontal grid, we interpolated all EURO-CORDEX RCM simulations and LAPrec historical data onto the HISTALP equidistant grid of 5 arcmin (~ 6 × 9 km) resolution. Both downscaling (for EURO-CORDEX datasets) and upscaling (for LAPrec dataset) interpolation was performed using the distance-weighted average remapping method in the Climate Data Operators software (CDO, https://code.zmaw.de/projects/cdo), which uses the four nearest neighbour values for each grid cell from the original field. According to Torma et al. (2015), the distance-weighted average method produces the most consistent spatial patterns among several interpolation methods (e.g. bilinear, bicubic, field conserving) applied to the Alps area. The interpolation was done as a first-step operation on the original field before any further computation (undercatch correction, bias correction, calculation of yearly precipitation sum and mean summer temperature).

The final 200-year time series of temperature and precipitation fields, comprising 105 years of historical data (1901−2005) and 95 years of bias-corrected projections (2006−2100), were used in the distributed modelling of envELA.

2.2 Environmental ELA distributed modelling

For the sake of comparison with other glacier-related modelling studies (e.g. Huss and Hock 2015; Huss et al. 2017; Zekollari et al. 2019), glacierised areas were obtained from the Randolph Glacier Inventory (RGI) v6.0 (RGI Consortium 2017), which provides a snapshot of glacier extents at the end of the summer 2003 (more than 96% of those glaciers in the Alps that are included in the RGI) for the Alpine area. During that year, there were 3892 glaciers covering an area of 2089 km2 and altitudinal range of 1269−4776 m above sea level (asl). Each individual glacier polygon centroid was rasterized to the resolution of HISTALP dataset (5 arcmin) and then used as a mask area for the calculation of envELAs.

The envELA for every “glacierised” grid cell was computed over the period 1901−2100 using the new precipitation/temperature diagram of Ohmura and Boettcher (2018). This precipitation/temperature curve is based on data from 104 glaciers worldwide and takes the form of a fourth-order polynomial. However, for wide applications such as ours, the approximation of the P/T relationship via a quadratic function, that was derived by Ohmura and Boettcher (2018), is more practical:

$${P}_{a}=5.87{{\mathrm{ T}}_{sum}}^{2}+230 {\mathrm{ T}}_{sum}+966 \, \mathrm{ S}.\mathrm{E}. =648\mathrm{ mm}$$
(3)

where \({P}_{a}\) is annual precipitation (mm) and \({\mathrm{T}}_{sum}\) is summer (JJA) air temperature (°C). This equation was applied to every 5 arcmin grid cell comprising the glacier centroids. First, the value of \({\mathrm{T}}_{sum}\) was obtained using the gridded LAPrec/RCM precipitation dataset (\({P}_{a}\)). Then, the height difference (ΔHTsum) between the simulated \({\mathrm{T}}_{sum}\) and HISTALP/RCM \({\mathrm{T}}_{sum}\) was computed by applying a summer environmental lapse rate of 0.65 °C 100–1 m representative for the greater Alpine region (Rolland 2003; Rubel et al. 2017) to the digital elevation model (DEM) that underlies the HISTALP climatology. Finally, the envELA was calculated by subtracting ΔHTsum from DEM (Fig. 2). In contrast to a generally clear relationship between temperature and altitude, the variations in precipitation with elevation are very complex and less robust (Frei and Schär 1998). The spatial variation in precipitation gradients with elevation has been found to vary from 2 mm m−1 at the northern and southern foothills of the Alps and in the outer regions, to 0−0.6 mm m−1 in the inner regions and southern slopes, while small negative gradients have also been detected (Schwarb et al. 2000). Because of these uncertainties, no additional adjustments to the precipitation of each grid cell were made except the one related to the equation’s standard error of 648 mm, which served to determine the envELA error.

Fig. 2
figure 2

The workflow of envELA calculation. See text for full explanation

In our computations of the envELA, we assume that the glacier centroid masks remain constant throughout the model simulation. In other words, we are essentially computing the lowest boundary of climatic glacierisation, which is dependent solely on climatic factors, i.e. annual precipitation and mean summer temperature.

2.3 Model validation

We validated the envELA time series pattern over two time periods using two ELA datasets representing the effELA as per the terminology used in this study: a) a dataset of glacier-wide end-of-mass-balance-year equilibrium-line altitudes from the WGMS FoG database (WGMS 2019) and b) annual end-of-summer SLAs derived from Landsat time series.

In the first validation exercise, we compared the envELA gridded data (5 arcmin spatial resolution) with point values from the WGMS dataset (a), both averaged over the entire study area for the historical climatological time series (1948−2017). The WGMS dataset used here (hereinafter referred to as “wgmsELA”) consists of 1255 ELA values spread over 70 years (1948−2017) for 62 glaciers, well distributed geographically and climatologically over the Alps (Fig. 1). Since time spans for individual glaciers in the WGMS FoG database were variable, we excluded those ELAs that were outside the glacier altitude range reported in the database, as well ELAs before the year 1948 and after the year 2017 with a too small annual sample size. However, to match the end of the envELA historical time series (1901−2005), only wgmsELAs for the period 1948−2005 were taken into account in the model performance analysis.

In the second validation exercise, we compared the envELAs derived from climate projections for the period 2006−2019 with the end-of-summer SLAs (b), both averaged over the subset area covering the western Alps (Fig. 1), where adequate images Landsat cloud-free were available. SLAs were extracted using a previously developed automated method (Racoviteanu et al. 2019) and Landsat data 5, Enhanced Thematic Mapper Plus (ETM +) and Operational Land Imager (OLI) (30 m spatial resolution), freely obtained from the USGS. We used level 2 (surface reflectance) data, which is orthorectified, analysis-ready imagery. We selected all the cloud-free images during the ablation season for each year (July to September) (a total of 50 images for the western Alps). Out of these, only 36 images were suitable for the calculation of the highest SLAs. Assuming that for mid-latitude mountain glaciers, superimposed ice can be considered negligible, the end-of-summer remote sensing SLA is an accurate approximation of the glacier ELA (hereafter referred to as “rsELA”) (Braithwaite 1984; Rabatel et al. 2005, 2013). The performance of the rsELAs over the subset region was checked visually using Landsat colour composites, as well as against manually-derived ELAs for 43 glaciers in the French Alps by Rabatel et al. (2013) for the period 2003−2010. The rsELAs were checked for errors due to shadows in the accumulation area of glaciers, and manually corrected when necessary. Uncertainties in the rsELAs were calculated as root mean square error (RMSEz) on the basis of DEM (EU-DEM v1.1 2017) accuracy, buffer size used to derive snowlines and differences with manually edited ELAs.

3 Results

3.1 Past climate and future projections

Here we present the analysis of mean summer temperature and annual precipitation, which are used as inputs for the envELA computations. The climate results are presented as average values over the area inside the Alpine Convention perimeter (Fig. 1) unless otherwise stated.

The increase in mean summer temperature is evident in the 50-year segment from 1901 to 1950, then a stable to slight cooling in the subsequent 25-year section until the mid-1970s, followed by a trend towards a strong warming in the last 25 years of the twentieth century. Mean summer temperature in the Alps in the period 1901−1930 was 13 °C, whereas it was 0.1 °C higher in the period 1911−1940, 0.5 °C higher in the periods 1921−1950, 1931−1960 and 1941−1970, only 0.2 °C higher in the period 1951−1980 and 0.4 °C higher in the period 1961–1990. In the period 1971−2000, mean summer temperature was 13.9 °C, which is the highest in the twentieth century. The coolest summer during the historical period (1901−2005) was recorded in 1913 with an average temperature of 11.7 °C and the warmest in 2003, when the average temperature over the Alps was 17.7 °C. Mean summer temperature anomalies (2071−2100 vs. 1971−2000) range from + 1 °C to + 2 °C under RCP2.6, from + 2 °C to + 4 °C under RCP4.5 and from + 4 °C to + 7 °C under RCP8.5. Projections also suggest that the warming will be stronger in the western Alps with respect to the eastern and this anomaly is particularly evident under RCP4.5 and RCP8.5 (Fig. 3).

Fig. 3
figure 3

a Summer air temperature in the period 1971–2000 for HISTALP historical datasets. Absolute change in summer temperature in 2071–2100 relative to reference period (1971–2000) for b RCP2.6, c RCP4.5 and d RCP8.5 scenarios. e 11-year centred running mean region-averaged summer air temperature of HISTALP temperature dataset (1901–2005) and bias-corrected EURO-CORDEX multi model mean (2006–2100) for RCP2.6 RCP4.5 and RCP8.5 simulations (Table 1)

No clear temporal and spatial trend in precipitation exists for the historical period (1901−2005). The only evident dry period was recorded near 1950 in correspondence with an increase of summer temperature. However, this period was also characteristic of increased transport of convective moisture at high elevations (Auer et al. 2007). The driest year during the historical period (1901−2005) was recorded in 1921 with 965 mm of precipitation and the wettest in 1960, when the average precipitation over the Alps was 1906 mm. Annual precipitation is projected to increase (2071−2100 vs. 1971−2000) by 5−10% in the eastern part and only by 0−5% in the western part of the Alps under RCP2.6. Projections suggest that under RCP4.5 and RCP8.5, the annual precipitation will decrease by 14% in the S-SW and increase by 22% in N-NE (Fig. 4).

Fig. 4
figure 4

a Annual precipitation in the period 1971–2000 for LAPrec historical datasets. Percent change in annual precipitation in 2071–2100 relative to reference period (1971–2000) for b RCP2.6, c RCP4.5 and d RCP8.5 scenarios. e 11-year centred running mean region-averaged annual precipitation totals of LAPrec1901 precipitation dataset (1901–2005) (corrected for gauge undercatch) and bias-corrected EURO-CORDEX multi model mean (2006–2100) for RCP2.6, RCP4.5 and RCP8.5 simulations (Table 1)

As for the climate across glaciated areas (see ELA calculation grid cells in Fig. 1), the summer mean temperature in the period 1901−1930 was 6.8 °C and 7.6 °C in the period 1971−2000. Projections suggest that the temperature will rise by 1.6 °C under RCP2.6, 2.8 °C under RCP4.5 and 5.4 °C under RCP8.5 by the end of the twenty-first century (1971−2000 mean vs. 2071−2100 mean). The mean annual precipitation in the period 1931−1960 was 1541 mm and 1592 mm in the period 1971−2000. About 4% increase in precipitation for all three emission scenarios is projected by 2071−2100 compared with 1971−2000.

3.2 Past and future environmental ELA evolution

Our reconstructions of past envELAs suggest that in the first 30-years of the twentieth century (1901−1930) the mean ELA over the Alps was 2980 m asl and stayed almost the same (2977 m asl) in the period 1911−1940. The subsequent period (1921−1950) showed a dramatic rise in envELA by 101 m with respect to 1901−1930, whereas the envELA stabilized during the periods 1931−1960 and 1941−1970 at 3063 m asl and 3066 m asl, respectively. The envELA in the period 1951−1980 rose only by 19 m with respect to the period 1901−1930 and by 60 m in the subsequent 30-year period (1961−1990). In the period 1971−2000, the mean envELA was 3094 m asl and rose by 1.6 m yr−1 between 1901−1930 and 1971−2000. In the last 30-years (1991−2020), the mean envELA was 3234 m asl (Fig. 5). This was calculated as average of the mean ELA for the period 1991−2005 (calculated from HISTALP and LAPrec datasets) and mean ELA for all three RCPs for the period 2006−2020 (calculated from RCMs dataset).

Fig. 5
figure 5

The envELA evolution (11-year centred running mean) for historical period (1901–2005) and various EURO-CORDEX RCM simulations (2006–2100) for a RCP2.6, b RCP4.5 and c RCP8.5. The red thick line corresponds to envELA computed with HISTALP and LAPrec historical datasets. Other thick lines are the RCM means for individual scenario and thin lines are individual RCM simulations (12 per RCP) (Fig. S2). The transparent bands correspond to the standard error of Eq. 3 (i.e. 648 mm). The numbers marked are 30-year mean envELAs (in m asl) for the periods 1901–1930, 1991–2020 and 1971–2100

Our projections, which are based on the multi model mean of 12 EURO-CORDEX RCM simulations (Table 2), suggest that the envELA will rise at a very similar rate across scenarios by 2050, whereas in the second half of the twenty-first century, the envELA is largely determined by the RCP. The mean envELA for the first half of the twenty-first century (2001−2050) is projected to reach 3283 m asl under RCP2.6, 11 m higher under RCP4.5 and 36 m higher under RCP8.5. For the second half of the century (2051−2100) the mean envELA is projected to rise to 3293 m under RCP2.6, 172 m higher under RCP4.5 and 457 m higher under RCP8.5. In contrast to envELA evolution under RCP2.6 and RCP4.5, a substantial part of the rise will take place in the second part of the twenty-first century under RCP8.5, and while the envELA seems to stabilize by 2100 under RCP2.6 and RCP4.5, it will continue to rise under the RCP8.5. By the end of the century (2071−2100), the mean envELA in the Alps rises to 3288 m asl under RCP2.6, 3484 m asl under RCP4.5 and 3880 m asl under RCP8.5 (Fig. 6). According to our calculations, the envELA in the Alps will rise by 194 m, 390 m and 786 m by the end of the century (2071−2100) under RCP2.6, RCP4.5 and RCP8.5, respectively, compared with the reference period 1971−2000, which amounts to a rise of 1.9, 3.9 and 7.9 m yr−1.

Fig. 6
figure 6

a Average envELA in the period 1971–2000. b envELA rise between 2071–2100 for RCP 2.6 relative to reference period (1971–2000). c envELA rise between 2071–2100 for RCP 4.5 relative to reference period (1971–2000). d envELA rise between 2071–2100 for RCP 8.5 relative to reference period (1971–2000). Black dots represent the remaining glaciers with envELAs below their maximum elevation. e Number and percentage of glaciers that will remain (lower part)/disappear (upper part) by 2071–2100 according RCP 2.6, RCP 4.5 and RCP 8.5

3.3 Environmental ELA variation across the Alps

Here we present past and future environmental ELA patterns across the Alps (west to east, and south to north). Figure 7 shows that the envELA rises from west to east with one anomaly between 7° and 9°E, where the envELA is on average lower due to higher precipitation amounts with respect to nearby sectors. The mean 1971−2000 envELA in the easternmost sector (Fig. 7e) was 2567 m asl, about 742 m lower than in the westernmost sector (Fig. 7a). The relative difference between the two periods when the envELA dropped during the twentieth century (around 1915 and 1960), is 50−130 m greater in the eastern Alps (between 9° and 15°E) compared to the western Alps (between 5° and 9° E). This corresponds to a larger change in precipitation in the east (3−6%) compared to the west (0−2%) (1910−1920 mean vs. 1950−1960 mean). While sectors in the east show a gradual rise in envELA since around 1960, the western sectors exhibit a steady envELA between the 1960s and 1980s, with only a slight drop towards 1980. The rise in envELA towards the end of the twenty-first century is more pronounced in the west compared to the east under all scenarios, but especially under RCP8.5 (+ 9 m year−1 in the westernmost sector vs. + 7 m year−1 in the easternmost sector). This is related to a negative precipitation anomaly and a higher positive temperature anomaly in the western Alps, the latter being especially evident under RCP8.5 (Figs. 3 and 4). We projected a rise in envELA for the periods 1971−2000 vs. 2071−2100 in the westernmost sector of + 899 m under RCP8.5, compared to + 695 m in the easternmost sector (Table S1).

Fig. 7
figure 7

(ae) The envELA evolution (11-year centred running mean) by longitude (every 2° between 5° and 15°E) for various EURO-CORDEX RCM simulations. Thin lines are individual RCM simulations (36 in total). Thick lines are the RCP means and the transparent bands correspond to multi model standard deviation. The red line corresponds to the envELA computed with HISTALP and LAPrec historical datasets. f Map of division of sections by longitude

Regarding the envELA variation by latitude, our model suggests that the envELA rises from north to south (Fig. 8, Table S2), with the mean 1971−2000 envELA of 2861 m asl in the northernmost sector (Fig. 8d), which is 520 m lower compared to the southernmost sector (Fig. 8a). We found different past ELA evolution patterns between the northern and southern sectors. Glaciers in the northernmost sector (Fig. 8d) had two envELA peaks of similar magnitude at around 1930 and just before 1950, while other sectors (Fig. 8a–c) showed only one prominent peak around 1950. The latter corresponds to a decrease of 10−12% in annual precipitation and an increase in summer temperature of + 0.4−0.5 °C at that time (1941−1950 mean vs. 1901−2005 mean), whereas the drop in annual precipitation and an increase in summer temperature in the northernmost sector was 5% and + 0.2 °C, respectively. While northern sectors (Fig. 8c–d) show gradual rise in the envELA since around 1965, the southern sectors (Fig. 8a-b) exhibit a slight drop in the envELA between 1955 and 1980. We project a faster rise in the envELA in the south compared to the north by the end of the twenty-first century (2071−2100 mean), under all scenarios. While the rise is very similar for all sectors under RCP2.6 (1.8−2.3 m year−1), there are larger differences among sectors under RCP4.5 (3.3−4.8 m year−1) and especially under RCP8.5 (6.9−9.5 m year−1). This is primarily related to the precipitation anomaly, namely an increase in precipitation by more than 5% under RCP4.5 and more than 10% under RCP8.5 in the north and an increase by less than 5% or even decrease by more than 5% under RCP4.5 and RCP8.5 in the south (Fig. 4).

Fig. 8
figure 8

ad The envELA evolution (11-year centred running mean) by latitude (every 1° between 44° and 48°N) for various EURO-CORDEX RCM simulations. Thin lines are individual RCM simulations (36 in total). Thick lines are the RCP means and the transparent bands correspond to multi model standard deviation. The red line corresponds to envELA computed with HISTALP and LAPrec historical datasets. e Map of division of sections by latitude

3.4 Model performance

The mean envELA for the period 1948−2005 was 3089 m asl, which is 110 m higher than the wgmsELA over the same time period (2979 ± 183 m asl). For the reference period 1971−2000, the mean envELA was 3094 m asl, 74 m higher with respect to the wgmsELA over the same time period (3020 ± 193 m asl). Despite a relatively low coefficient of determination (r2 of 0.47), the envELA reproduces fairly well the time series pattern of the wgmsELA for the period 1948−2005 (Fig. 9a) with a root-mean-square error (RMSE) of 190 m and a median absolute deviation (MAD) of 140 m.

Fig. 9
figure 9

a Comparison between the envELA averaged over the grid cells as defined in Fig. 1 and average regional wgmsELA for 62 glaciers for the period 1948−2017. The last 12 years (2006−2017) of the envELA (dashed line) represent average of all three RCPs envELAs. Thick lines correspond to 11-year centred running mean, thin lines represent yearly variations, and transparent bands correspond to the standard error of Eq. 3 (i.e. 648 mm) for the envELA and sample standard deviation for wgmsELA. (b) Comparison between the rsELA derived from Landsat images and envELA derived from climate projections under three different scenarios (RCP2.6, RCP4.5 and RCP8.5) for the period 2006−2019 over the rsELA extent as defined in Fig. 1

For the time period over which we validated with the rsELA (2006−2019), the mean envELA averaged over the subset area of the western Alps was 3314 ± 252 m asl, 3288 ± 244 m asl and 3315 ± 242 m asl (uncertainties are reported as multi-model mean standard deviations) under RCP2.6, RCP4.5 and RCP8.5, respectively. The rsELA derived from the Landsat time series, averaged over this subset area over the same time period was 3165 ± 44 m asl. The year to year variability shows that the overall temporal pattern of the rsELA for the period 2006−2019 is best reproduced by the envELA derived under RCP2.6 (RMSE of 168 m, MAD of 127 m and r2 of 0.24) (Fig. 9b). Under this scenario, the mean envELA is on average 149 m higher than the rsELA over the 2006−2019 period. A comparison of the rsELA with manually derived ELAs (based on end-of-ablation season SLAs) from 43 glaciers in the French Alps over the period 2003–2010 (Rabatel et al. 2005) yielded an RMSE of 41 m based on year-to-year comparison, which is within the uncertainty calculated for the rsELA; the temporal pattern of the rsELA is reproducing well the manually derived ELAs.

4 Discussion

4.1 Environmental versus effective ELA

The results presented here concern the envELA over an entire mountain range, based purely on climate data, which represents the first study of this kind. In contrast, the effELA that is here represented by both wgmsELA and rsELA, and was used to validate the envELA, takes into account the topoclimatic effects as well. Consequently, the envELA is ~ 75−150 m higher than the regional effELA when averaged over a longer climate period, e.g. 15−30 years. The difference between envELA and wgmsELA might also be related to the fact that the glaciers selected for measuring annual glacier mass balance by the WGMS are those with easier access and thus tend to be located on lower altitudes, consequently having on average lower ELAs. While the envELA is consistently higher than effELA, the time series pattern is in general well reproduced (Fig. 9). However, when comparing year-by-year values, we note some “extreme” values in the envELA with respect to the regional effELA. In the year 2003, the envELA is the highest recorded over the historical period 1901−2005 (Fig. 9a) and overestimates the regional wgmsELA by 668 m. The 2003 envELA is 4.1 standard deviations above the 30-year mean (1988−2017), whereas the regional wgmsELA deviates only by 2.0 standard deviations. This anomaly is related to the 2003 summer heatwave in Europe (Fink et al. 2004), when average summer temperatures were far above the long-term mean and were accompanied by annual precipitation deficits of up to 300 mm (Parry et al. 2007). Since the envELA is computed as a direct function of these two climate variables, it may immediately react to these climate anomalies. This is not the case for effELA, which is a result of the snow accumulation and metamorphism from previous years as well. Winter accumulation in steep alpine terrain can be heavily dominated by avalanches and wind-blown snow supply (Winstral et al. 2002; Dadic et al. 2010; Purdie et al. 2015). As a consequence of both of these drivers the additional accumulation on a glacier can be more than two times the winter precipitation (Kuhn 1995; Lie et al. 2003), which will enhance snow thickness locally and lower the effELA. Small glaciers are particularly susceptible to these two phenomena and their ELA can become progressively decoupled from climate (Dahl and Nesje 1992; Colucci 2016; Scotti and Brardinoni 2018). However, a reliable estimate and small deviation from the effELA can be achieved once the envELA is averaged over a longer climate period.

4.2 Model uncertainties

Model uncertainties are mainly related to the description of cloud and precipitation processes and to possible dependencies on elevation (e.g. Frei and Schär 1998), which are not included here. In addition, the input for the envELA model is annual and not winter precipitation, since the former represents annual accumulation better because precipitation in summer is often solid at high altitudes (Ohmura and Boettcher 2018). However, even at high altitudes a small portion (on average 9% at ELA according to Ohmura and Boettcher (2018)) of the total annual precipitation will fall in liquid form. Even though part of the liquid precipitation might remain in the snow and firn layers as superimposed ice, this amount represents a potential loss that would be unaccounted for as a contribution to the accumulation. Model uncertainties related to precipitation are relevant for both historical and future envELAs, whereas an even bigger source of uncertainty for the latter is probably the ensemble spread (for both temperature and precipitation) within each RCP.

4.3 Link between environmental ELA and historical observations of glacier fluctuations

The envELA is a direct, undelayed signal of climate change, whereas glacier length, and to some extent the effELA, lag their equilibrium response to a climate trend. The lag in geometric glacier response has been estimated at 50 ± 28 years for glaciers in the Alps longer than 1 km (Zekollari et al. 2020). Glaciers shorter than 1 km and ice bodies of the smallest size (i.e. ice patches and glacierets with an area of 0.1−0.5 km2), which account for almost 90% of all glaciers in the Alps (RGI Consortium 2017), are likely to respond much faster under changing climate (Kuhn 1995; Huss and Fischer 2016). These ice bodies are often so small that can repeatedly lie entirely in the ablation or accumulation zone on an annual basis. Despite the lag in glacier response, the envELA variation can be correlated, to some extent, with the front variation (e.g. Colucci 2016) and effELA of individual glaciers (Fig. 10), but quantitative comparison of the first two is not presented here. Although the majority of large valley glaciers in the Alps have retreated continuously since around 1865 (Painter et al. 2013), there are three documented advances of smaller mountain glaciers in the 1890s, 1920s and 1970−80 s (Zemp et al. 2008). The two advances in the twentieth century can be recognized from the drop in the envELA around that time. The first drop in the envELA is recorded from the beginning of our studied period until around 1915, when the envELA was the lowest on record in the twentieth century. This is followed by a substantial and constant rise until the 1950s, when the envELA starts to decline again. A similar rising pattern in the effELA until the 1950s can be recognized for the Claridenfirn and Silvretta glaciers (Fig. 10). The second drop takes place for about 10 years and is followed by ~ 25-year period of stable values. While the effELAs of most of the glaciers presented in Fig. 10 follow the envELA curve and display in general stable values for the period between 1960 and 1985, when many Alpine glaciers were relatively close to equilibrium (Zekollari et al. 2020), few others (e.g. Allalin, Hohlaub, Schwarzberg) exhibit a substantial drop in the effELA in the 1970−80 s, when many smaller mountain glaciers slightly readvanced in the 1980s (Zemp et al. 2006, 2008). After that period, a rise in the envELA occurs until the end and beyond the historical period. Over this period, glaciers showed an area reduction of about 22% mainly occurring after 1985 (Zemp et al. 2006). Around 2010, glacier retreat in the Alps started to catch up and outpace the climate warming signal (Zekollari et al. 2020). The main processes leading to glacier decline are disintegration and down-wasting (Paul et al. 2004). Only glaciers of the smallest size show any sort of equilibrium, mostly in areas with high mean annual precipitation (e.g. Colucci and Žebre 2016; Scotti and Brardinoni 2018).

Fig. 10
figure 10

The wgmsELAs of 14 glaciers located in the sectors a between 7° and 9°E longitude and b 9° and 11°E longitude. Lighter lines equate to northern located glaciers and darker lines to southern located glaciers. In red is the envELA averaged over that same sector. All ELAs are presented as 11-year centred running means

4.4 Comparison with modern and future ELA projections from other studies

Results from previous global studies (Huss and Hock 2015; Huss et al. 2017) rely on a method in which geometry changes are parameterised; they suggest that the regionally averaged ELA was at 3147 m asl in 2010 in Central Europe and at ~ 2900−3000 m asl in the period 1980−2000 in the Alps. The same studies projected a rise in ELA by ~ 100−500 m depending on the emission scenario by the end of the twenty-first century. Our model suggests that the envELA in the Alps in 2010 was on average at 3190 m asl and between 3023 and 3291 m asl for the different RCPs, the mean 1980−2000 envELA was 3131 m asl, while the envELA is projected to rise by 91−683 m (1981−2010 mean vs. 2071−2100 mean). According to the regional study by Zemp et al. (2007) the mean ELA of all Alpine river basins over the period 1971−1990 was 2951 m asl, 107 m lower than our envELA for the same period. This difference can be attributed to a different empirical relationship between temperature and precipitation and higher (3 arc sec) spatial resolution used in the ELA distributed modelling in the work of Zemp et al. (2007). Unlike the study of Huss and Hock (2015), 35 glaciers in Pyrenees, Apennines and Balkans are not taken into account in our calculations. These glaciers are located at very low altitudes, where the topoclimatic effect plays a dominant role in their existence (e.g. Hughes 2010; Serrano et al. 2011; Colucci and Žebre 2016; Gachev et al. 2016), contributing to the lower regional ELA calculated by Huss and Hock (2015). In addition, in Huss and Hock (2015), the ELA is considered as a regional effELA, and this generally underestimates the envELA. Dissimilarities might also be attributed to different climatological time series used for the ELA calculations. Huss and Hock (2015) projections are based on 14 individual GCM time series from the Coupled Model Intercomparison project Phase 5 (CMIP5) (Taylor et al. 2012), while ours relay on 12 EURO-CORDEX RCM simulations of higher resolution. However, the trend of our future envELA (Figs. 5 and S3) agrees well with the results of Huss and Hock (2015), who projected a slight decrease in the rate of ELA rise beyond 2050 under RCP2.6 and stabilization under RCP4.5, and a nearly linear rise throughout the projection period under RCP8.5 (Fig. S9 in Huss and Hock (2015)).

4.5 Implications for future glacier state

Most glaciers in the Alps are in disequilibrium with current climate (e.g. Christian et al. 2018) and part of their future mass loss is committed even without further warming (e.g. Bahr et al. 2009; Mernild et al. 2013; Marzeion et al. 2017, 2018; Santin et al. 2019). In 2001, the committed volume loss was 34% of the total ice volume at that time, increasing to 42% in 2010 (Zekollari et al. 2020). A glacier volume loss of about 65–80% under RCP2.6, 80–90% under RCP4.5 and 90–98% under RCP8.5 is projected for the period between the early twenty-first century and 2100 by different models (Marzeion et al. 2012; Radić et al. 2014; Huss and Hock 2015; Zekollari et al. 2019). Here, we assume that a glacier becomes unviable and likely to disappear when the envELA lies above its maximum altitude, meaning on average no accumulation area on a glacier for at least a 30-year period. At the date of the Randolph Glacier Inventory v6.0 (RGI Consortium 2017) (~ 2003), 51% of all glaciers in the Alps were already placed completely below the mean envELA over the period 1974−2003 (3113 m asl). Our projections suggest that under RCP2.6, RCP4.5 and RCP8.5, the envELA will exceed the maximum elevation of 69%, 81% and 92% of the glaciers in the Alps by the end of the twenty-first century (2071−2100 mean), respectively (Fig. 6e). Therefore, glacier retreat and disappearance in the Alps will potentially have a considerable impact on river runoff and hydrology.

5 Conclusions

Using a simple model, we explored 200-years of the environmental ELA evolution over the entire Alps. Our results are in line with available observations and qualitatively in line with previous work employing more sophisticated models that include glacier dynamics. We find that the environmental ELA between the beginning and end of the twentieth century (1901−1930 mean vs. 1971−2000 mean) rose by 114 m mainly due to an increase of 0.8 °C in summer temperature, whereas the annual precipitation stayed largely unchanged. Our results suggest that projected summer warming of 1.6 °C, 2.8 °C and 5.4 °C under the RCP2.6, RCP4.5 and RCP8.5 scenarios, respectively, will likely drive significant environmental ELA rises of about 194 m, 390 m and 786 m by the end of the twenty-first century (2071−2100 mean) for glaciers in the Alps compared with the period 1971−2000. A slight increase (~ 4%) in annual precipitation for all three emission scenarios, derived from a winter increase and a summer decrease, will only have a minor effect on the average environmental ELA over the entire Alps. However, local precipitation anomalies will drive different changes in environmental ELA over different regions of the Alps, suggesting in general a faster environmental ELA rise in the west and south with respect to the east and north. These future projected levels of the environmental ELA will result in the unviability and possible subsequent disappearance of between 69 and 92% of all glaciers in the Alps by the end of the twenty-first century. This will have substantial implications for water resources in the region, particularly during the summer season.

This is the first study to investigate the environmental ELA of the entire Alps over such a long period and provides a strong basis to better understand regional differences in glacier response to a changing climate. Our modelling approach needs a minimum set of input data and can thus be easily applied to larger areas and longer time periods. This approach can be particularly useful in Quaternary science for quick and robust palaeo-ELA reconstructions. However, there is still potential to improve the simulations, for example by incorporating solar radiation as the third variable in the model and increasing the spatial resolution of the input data.