1 Introduction

The growing season of alpine plants above the climatic treeline is closely tied to snowmelt dates (Körner 2021; Hülber et al. 2010; Wipf and Rixen 2010). Therefore, the mostly long-lived alpine perennials only have about three months for their seasonal development, but there is a high local variability due to the complex topography and substantial snow redistribution by wind in winter. Most of an alpine plant’s growth occurs during the first few snow-free weeks (Körner 2021). In the core of the alpine vegetation belt of the Swiss Alps at 2500 m asl, peak aboveground biomass is commonly reached in early August; thereafter, photoperiod signals (daylength) trigger the resorption of nutrients from leaves (Prock and Körner 1996), followed by plant tissue senescence. The meristems of alpine plants (tissues for growth) are located a few centimeters below the ground surface and insulating snow during winter prevents soil freezing (Körner 2021). Freezing of alpine soils may alter nutrient relations by physical degradation of soil organic matter and retention of nutrients such as inorganic nitrogen (Edwards et al. 2007; Freppaz et al. 2008; Baptist et al. 2010a). Changes in the duration of the snow cover, particularly in snowmelt dates are thus relevant for plant development and nutrient cycling in alpine ecosystems. However, not all alpine species may profit from earlier snowmelt, as their development may be constrained by short photoperiods (Heide 2001; Keller and Körner 2003).

1.1 Decline in the snow cover

Recent climatic changes have affected the seasonal snow cover (Hock et al. 2019) and raised the concern of various stakeholders (Beniston et al. 2018), as the snowpack is crucially important for skiing and meltwater runoff affects hydropower production (Westaway 2000). In Switzerland, the number of days with snowpack reached a minimum in the late 1980s, based on 100-year snow series from five stations between 450 and 1860 m asl (Scherrer et al. 2013). These reductions were mainly explained by rising temperatures (Scherrer et al. 2004). By the end of the twenty-first century, a continuous snow cover of more than 30 days will be rare in the Swiss Alps below 1600 m asl, and the mean snow depth (Dec–Feb) is predicted to decline by 74%. Even the high-elevation station Weissfluhjoch at 2540 m asl may undergo reductions in snow depth of 23% (Schmucki et al. 2015a). This station’s snow depth has been considered to be representative for high-elevation sites (2000–2900 m asl) in Switzerland (Marty and Meister 2012).

Past trends in snowmelt dates, which are crucial for the alpine plant development, were assessed by Klein et al. (2016), but the analysis was restricted to 45 years and to one site in the alpine vegetation belt only. Furthermore, to predict the functioning of alpine ecosystems under climate change, it is necessary to project future snowmelt dates for the alpine vegetation zone.

1.2 Main objectives

We assess trends in snow onset and snowmelt between 1958 and 2019, focusing on 38 snow measurement stations between 1055 and 2555 m asl in the Swiss Central Alps. This is a core region for winter tourism, hydropower generation, but also longer-term research projects (Körner et al. 2019). Because most time series for snowmelt dates are not longer than 20 years at high-elevation sites, we elaborate the correlation of snowmelt dates at 15 automatic stations above 2100 m asl (daily median of 30 min intervals) with manually measured long-term series of daily snow depth data at 23 stations (mainly below 1500 m asl). A main goal of this study was to evaluate the role of air temperatures and of the mean winter snow depth on snowmelt dates in order to predict future trends in snowmelt dates. Owing to fast melting of the snowpack in spring, we hypothesise that snowmelt dates are largely determined by the time when mean air temperatures exceed a threshold, yet to be identified. We expect that the mean winter snow depth regulates snowmelt dates in years with exceptionally high or low snow depth. We then apply the identified drivers of past snowmelt dates to the Swiss climate change scenarios (CH2018) to predict the onset of the alpine growing season throughout the twenty-first century.

2 Methods

2.1 Station types

We analysed snow depth measurements of 38 stations between 1055 and 2555 m asl in the Swiss Central Alps to derive time series of dates of snow onset and melt (see Table 1 for list of stations). We used two types of stations: manually operated (MAN) and automatic stations. MAN stations are often below 1500 m asl, operated by the Federal Office of Meteorology and Climatology MeteoSwiss and by the WSL Institute for Snow and Avalanche Research (SLF). Snow depth readings are taken daily and many time series date back to the 1950s. Automatic stations were established to improve the Swiss avalanche warning system (Intercantonal Measurement and Information System IMIS). IMIS stations are typically located above 2000 m asl, often close to avalanche starting zones. For most of these stations, the first fully operational year was 2000. In addition to wind and temperature (air and ground), snow depth is measured by an ultrasonic sensor (SR50, Campbell Scientific, US; every 30 minutes). Both station types are typically located in relatively flat terrain. We utilised daily snow depth measurements from 23 MAN stations (morning readings) and 15 IMIS stations (daily median of 30 min snow depth values; Table 1). We included two stations outside of the Swiss Central Alps because of their long-term data series since 1936 and 1952, respectively (Fig. 1: 5WJ, 4MO).

Table 1 General description of all 38 stations used for identifying the snow onset and melt date, with elevation, coordinates and operational years in the Swiss Central Alps
Fig. 1
figure 1

Map of manual (MAN) and automatic (IMIS) stations and the elevational belts in central Switzerland (villages and cities are presented by an asterisk; digital elevation model of the CH2018 climate change scenarios, resolution 2x2 km)

2.2 Definition of snow onset and melt

We defined snow onset and melt dates as the first and the last day of the main snow period. The main snow period corresponded to the longest period of snow cover in a hydrological year (e.g. year 1991 runs from 1 Sept 1990 to 31 Aug 1991). Owing to the two measurement techniques for snow depth and the distinct elevations, we performed different approaches to identify the main snow period of MAN and IMIS stations, respectively.

At lower elevations with MAN stations, snow-free days during winter are not uncommon; therefore, the main snow period is less clearly identifiable. For these stations, we allowed the main snow period to contain days without snow cover, but not more than 5 consecutive days. We also tested maximum gaps of 1, 3 and 10 days.

IMIS stations at high elevations with low winter temperatures show a distinct main snow period. Thus, it was unnecessary to include snow-free days for the main snow period. However, the ultrasonic signal is scattered and often exhibits an offset, impeding a clear differentiation between snow cover and ground surface. Therefore, we implemented a procedure that recognises this offset and then deduces the main snow period, including snow onset and melt dates (R function in Online Resource 1). To suppress the scattering, we used 5-day running means of snow depth. The snow depth signal just after snowmelt was usually less variable than for snow cover; thus, we looked for periods with relatively constant ultrasonic signals to identify the ground surface (offset). Low standard deviations (SD) of the 5-day running windows, both in winter and summer, indicated a period of constant snow depth or a period without snow (ground surface). Of all winter periods with low SD that were longer than 20 days (empirically found), we selected the most frequently occurring snow depth to obtain the height of the ground surface (offset).

At Weissfluhjoch, a MAN (5WJ) and an IMIS (WFJ2) station run in parallel (Table 1), which enables us to compare the snow onset and snowmelt dates deduced from the definitions for the two station types.

2.3 Data quality and handling of data gaps

At MAN stations, the method for measuring snow depth remained unchanged during the observational period. Displacements of the measuring fields have rarely occurred and are not influential (Buchmann et al. 2020). In addition, there were no changes in the sensor type of the IMIS stations nor in measuring fields through dislocations. We visually checked the snow depth data of the IMIS stations and removed obvious outliers from the data series.

For both station types, the few dates with missing snow depth data were interpolated linearly. Overall, 5% of all observations between 1 October and 1 June were missing. However, some of the MAN stations report missing values over the snow-free season. For the period between the determined snow onset and melt date, only 0.6% of the data was missing.

Missing snow onset or melt dates of single years were either filled based on correlations with nearby stations at similar elevation and comparable snow depth, or if snow depth towards the end of the main snow period already went below 20 cm, snowmelt dates were extrapolated linearly (Online Resource 2 lists stations and years with interpolated data gaps). In total, our data set included 1198 observations for snowmelt dates, 2.3% of which were missing.

2.4 Statistical analyses

Any trends in time series (onset and melt) were analysed with Theil–Sen slopes and the Mann–Kendall test (R package trend v1.1.2, Pohlert 2020), both relatively robust against outliers. To calculate long-term trends in snow onset and melt dates, we chose the period 1958 to 2019, as it was the longest possible period with snow depth data including two stations above 2000 m asl. This period included 16 out of 23 MAN stations between 1055 and 2555 m asl with continuous measurements for snow depth.

As snowmelt dates are more decisive for alpine plant development than snow onset, we focussed on snowmelt dates for further analysis. However, both snow onset and snowmelt dates were needed to delineate the main snow cover duration. Short-term trends in snowmelt dates were computed as Theil–Sen slopes by running windows of 21 years for all 23 MAN stations, each covering at least the period 1958–2019 except for GOT at 2107 m asl (1936–1970). The trend in year y corresponded to the trend for the years y–10 to y+10. IMIS stations could not be included in the analysis of running window trends, as they did not cover the required 21 years. Thus, we compared high (IMIS stations) and low-elevation (MAN stations) sites by the Pearson correlation between their snowmelt dates. For each station we computed the yearly deviation of snowmelt dates from its mean between 2000–2019. Thereof, we calculated the yearly mean and SD of the deviation for each station type. A prerequisite for the correlation of time series is the so-called weak stationarity, where the expected mean, variance and covariance do not change over time. A commonly adopted approach is to subtract from each value the value from the previous year (Cryer and Chan 2008). For the so derived first-differences of the snowmelt series, we then calculated the Pearson correlation between snowmelt dates at MAN and IMIS stations.

To project future snowmelt dates, we unravelled the contributions of temperature and snow depth to snowmelt dates by fitting a model including both parameters and their relationships with snowmelt dates. For snow depth, we used the mean snowpack height HSmean between snow onset and snowmelt date, available from the daily snow depth measurements. Automatic weather stations (MeteoSwiss) measure the daily mean air temperature 2 m above ground. For nine of the MAN stations, air temperature measurements from the same location (village), usually within less than 100 m elevational difference were available: GOS (27; means 27 years of overlapping snow depth and temperature data at a given site), GTT (37), 4UL (39), GRH (50), GUE (58), DIS (59), 5WJ (60), 4MO (68) and 2AN (79). The temperature data have already undergone a semi-automatic quality control, and we used homogeneous data series for GTT, 4UL, GRH, GUE, DIS, 5WJ and 2AN. The 14 IMIS stations recorded air temperatures every 30 minutes during 20 years (WFJ2 was excluded due to the parallel station 5WJ, see Table 1 for abbreviations).

Snowmelt dates differ between stations (elevations) and fluctuate over the years. Thus, the time when temperatures are high enough for snowmelt to occur also varies between stations and years. As we observed that the snow depth rapidly decreases during the last weeks before the snowmelt date, we assumed that the snowmelt date would only occur when air temperatures were high enough for a certain time period. We thus adopted a temperature measure for the first day (day of year DOY) when the daily mean air temperature exceeded 5 °C in a running window of 21 days (T21d ≥ 5 °C). We empirically found a high correlation of T21d ≥ 5 °C with the snowmelt dates (r = 0.87, p < 0.001). As snowmelt occurs at temperatures above 0 °C we also tested a threshold of 0 °C and 3 °C, and for an upper limit we used a threshold of 10 °C. Moreover, we tested running windows of 7, 14 and 30 days, and also degree-days with thresholds >0 °C and > 5°C, all yielding weaker correlations (shown in Online Resource 3). Moreover, our prior exploration of the data implied a high linear correlation of the snowmelt dates with the square root of HSmean (r = 0.85, p < 0.001). We thus fitted the linear mixed model:

$$ {SM}_{station}=a\ast {T}_{21d}\ge 5{}^{\circ}C+b\ast \sqrt{HS_{mean}}+{intercept}_{station} $$
(1)

where a indicates how much the snowmelt date was delayed by an additional day of T21d ≥ 5 °C (unit day day–1); b (day cm–1/2) describes the non-linear relation of snowmelt dates and the square root of HSmean, and interceptstation is the intercept of stations used as random factor. We visually checked residuals of the linear mixed model and excluded 33 data points (outliers) with large residuals and high leverage from the analysis. Because T21d ≥ 5 °C and HSmean may be correlated, we also calculated the variance inflation factor. This factor indicates how much variance of a variable is influenced (inflated) due to collinearity with another variable. To assess their individual effect on snowmelt dates, we also fitted linear mixed models for either of the two fixed effects, HSmean or T21d ≥ 5 °C. For each model, we derived the marginal R2m that showed the percentage of variance explained by the individual fixed effect, along with the conditional R2c that explained the variance by the fixed and random effects combined (R package MuMIn v1.43.15; Bartoń, 2020). R2m enabled us to quantify the contributions of both, HSmean and T21d ≥ 5 °C to snowmelt dates.

2.5 Snowmelt dates in the twenty-first century

To estimate trends in snowmelt dates for the twenty-first century, we applied the above linear mixed model to the CH2018 Climate Change Scenarios for Switzerland (CH2018 Project Team 2018). These include three Representative Concentration Pathway scenarios, referred to as: RCP2.6 (strong emission reductions, warming limited to 2 °C, 12 simulations), RCP4.5 (declining emissions after 2050, >2 °C warming compared to pre-industrial times, 25 simulations) and RCP8.5 (unabated emissions, 31 simulations). For these simulations, a total of 31 regional climate models of the EURO-CORDEX (Jacob et al. 2014) were used, all exhibiting resolutions of 12 or 50 km. Downscaled scenarios (quantile mapping method) are available for seven of the MAN stations considered here (2AN, 4MO, 4UL, 5WJ, DIS, GRH, GUE), and also as a 2×2 km grid over Switzerland (CH 2018). For further analysis, we used the scenarios for these seven MAN stations (between 1197 and 2540 m asl) and for the IMIS stations, we extracted the simulation parameters from this 2×2 km grid. All simulations provided projections for air temperature, thus T21d ≥ 5 °C was estimated accordingly. To account for elevational differences between IMIS locations and grid cells, we applied an air temperature lapse rate of –0.5 K per 100 m increase in altitude (results for lapse rates between –0.2 and –0.8 K per 100 m are provided in the Online Resource 4). As the scenarios include projections for precipitation, we deduced HSmean from accumulated winter precipitation. For the seven MAN stations, we first determined the sum of the measured daily precipitation starting at 1 October until T21d ≥ 5 °C (DOY). We only considered days with daily mean temperatures ≤ 0 °C during the period between 1982 and 2019. We calculated Pearson correlations of the so derived cumulative precipitation and of the HSmean for the seven stations. We then computed a station-specific proportionality factor of the cumulative precipitation in relation to the HSmean (mean accumulated precipitation per cm of snowpack height). For IMIS stations, which commonly do not measure liquid precipitation, we used the mean proportionality factor of the seven MAN stations (1197–2540 m asl).

2.6 Model validation

To assess the model’s performance, particularly with measured HSmean versus parameterised HSmean (Eq. 1), we refitted the linear mixed model with empirical data of the seven MAN stations from uneven years only. We subsequentially used the measured and the parametrised HSmean to predict snowmelt dates of uneven years and compared the model results to the observed snowmelt dates. In a second step, we applied the original model to the reference periods of the climate change scenarios, 1990 to 2019 (MAN), and 2000 to 2019 (IMIS), respectively. For these 30 and 20 years, we compared estimated and observed snowmelt dates as well as trends (Theil–Sen slopes) of the snowmelt date time series (estimated vs. observed).

We then projected future changes in snowmelt dates for the three 30-year periods 2020–2049 (beginning of century), 2045–2074 (mid-century), and 2070–2099 (end of century) as well as for the twenty-first century (2020–2099) to achieve future longer-term trends of snowmelt dates, determining the onset of the alpine growing season. All calculations and analyses were performed with R 3.6 (R Core Team 2019).

3 Results

3.1 Comparison of definitions for snow onset and snow melt

We carefully assessed the suitability of our definitions for snow onset and melt dates. The main snow period at MAN stations was only slightly affected by the maximum number of consecutive snow-free days we allowed for. In 53% of all cases, snowmelt dates were the same, irrespective of the number of snow-free days. The 5-day-gap we used for further analysis yielded the same snowmelt dates as gaps of 3 days in 86% of all cases.

Snow onset at both station types at the Weissfluhjoch (Table 1) differed by no more than one day in 11 out of 19 overlapping years. However, late autumn snowfall in two years was classified as part of the main snow period at the MAN station, resulting in a difference of 10 and 17 days compared to the IMIS station. Snowmelt was slightly earlier at the MAN than at the IMIS station, with a median difference of one day, and a maximum difference of four days in 2000 and 2010. This variation in snowmelt dates was most likely caused by micro-topographical differences (slight topographical depression at the IMIS station).

3.2 Trends in snow onset and melt dates: 1958 to 2019

Out of the 38 stations, 16 MAN stations had continuous snow depth recordings for the past 62 years or more. To assess relevant shifts across time and elevations, we calculated the trends in snow onset and melt dates for the period 1958 to 2019. There was a tendency towards later snow onset that was indicated by persistently positive Theil–Sen slopes. These trends were statistically significant for two stations only: 2ME, 6RI (Table 2). In contrast, there was a clear trend towards earlier snowmelt and for ten stations, this shift was statistically significant (p < 0.05). Across all stations, snowmelt on average occurred 2.8 ± 1.3 days (mean ± sd) earlier per decade. The magnitude of the trends in snowmelt did not depend on elevation (p = 0.46; R2 = 0.04) and at the two stations in the alpine vegetation zone (GUE, 5WJ), snowmelt occurred 3.1 and 4.0 days earlier per decade (Table 2). With 6.3 days per decade, the station 2ME (1320 m asl) showed the most distinct shift in snowmelt dates.

Table 2 Trends (Theil–Sen slopes) for the snow onset and snowmelt dates between 1958 and 2019 with the corresponding p-value (Mann–Kendall test)

3.3 Trends in snowmelt dates: 21-years running window trends

Trends of 21-year running windows revealed the temporal variability of snowmelt date trends between 1946–2009. We required continuous snowmelt date time series of at least 21 years, excluding the IMIS stations. Our analysis unveiled a period with later snowmelt dates during the 1940s to 1980s, but trends were rarely statistically significant. Strikingly, the late 1980s and early 1990s were denoted by a sharp, significant trend towards earlier snowmelt, with a mean of 15.1 ± 3.8 days per decade in 1988 (1978–1998, nstations = 21; Fig. 2). This pattern was evident for nearly all stations and along all elevations, even at the highest location (5WJ, 2540 m asl). Only two stations (1MR, 6CB) did not show this shift. After the 1990s, trends in snowmelt dates did not exhibit a consistent direction and they were mainly statistically non-significant. The clear trend towards earlier snowmelt for the whole period between 1958 and 2019 was caused by the abrupt shift during the late 1980s and early 1990s.

Fig. 2
figure 2

Snowmelt trends in 21-year running windows at MAN stations, ordered by elevation. The colours indicate the direction of the trend. The symbol size corresponds to the slope of the trend (expressed in days per decade) and opaque fillings represent statistically significant trends at p < 0.05 (see Table 1 for station abbreviation and full station name, respectively). For instance, a trend for the year 2009 (last dot) corresponds to the Theil–Sen slope for the period 1999 to 2019

3.4 Comparison of snowmelt dates at MAN and IMIS stations

As high-elevation stations were scarce in the past, we compared snowmelt dates of IMIS and MAN station during 2000 and 2019, with measurements for both station types. Between 2000 and 2019, average snowmelt at the IMIS stations (2101–2555 m asl) was in mid-June around the longest day of the year, whereas at the MAN stations (1055–1980 m asl without GUE and 5WJ) snowmelt occurred in mid-April. The Pearson correlation for the yearly mean snowmelt dates of the two station types was 0.83 (nyears = 20, p < 0.01; Fig. 3). Even without the two MAN stations above 2000 m asl (5WJ and GUE), this correlation was strong (0.82, nyears = 20, p < 0.01). This means the changes in snowmelt dates at stations above the climatic treeline (IMIS stations) are closely related to these at lower-elevation MAN stations. There was one exception in 2002 with the second warmest February since the start of measurements at that time, followed by a warm spring, causing early snowmelt at low elevation (12 days earlier than the mean at MAN stations). Whereas at high elevation, large precipitation amounts in February, March and May 2002 fell as snow, thus delaying snowmelt (six days later than the mean at IMIS stations).

Fig. 3
figure 3

Yearly deviation of snowmelt dates from the 10-year average (2000–2009) for manual (MAN; (red, n = 22) and automatic (IMIS) stations (blue, n = 15; mean [line] ± sd [bands]). The deviation at single stations is shown in Online Resource 5

3.5 Influence of HSmean and temperature on snowmelt dates

Air temperature and HSmean, as a proxy for accumulated winter precipitation, directly influenced snowmelt dates along the whole elevational range of the stations. We quantified both effects by analysing snowmelt dates in relation to the observed T21d ≥ 5 °C and HSmean (Fig. 4; stations at 1055-2540 m asl; variance inflation factor of 1.04). For each week that this 5 °C threshold was attained earlier, mean snowmelt dates occurred 3.9 days earlier (a = 0.6 days day–1 in Eq. 1). The model with T21d ≥ 5 °C explained 44% of the variance in snowmelt dates (Rc2 0.81). A higher HSmean delayed snowmelt, proportionally to the square root of HSmean (Fig. 4b; b = 4.6 days cm −1/2 in Eq. 1). For instance, an increase of HSmean from 40 to 90 cm caused a delay in snowmelt dates by 14 days, but an increase from 140 to 190 cm postponed snowmelt by nine days only. In the model with HSmean only, still 30% of the variance in snowmelt dates was explained (Rc2 0.80). Together, T21d ≥ 5 °C and HSmean explained 80% of the variance in the snowmelt date (Rc2 = 0.89).

Fig. 4
figure 4

Snowmelt a) in relation to the first day (day of the year; date) when the running mean air temperature of a 21-day window reaches a threshold of 5 °C (T21d ≥ 5 °C), and b) in relation to the mean snow depth between snow onset and melt date (HSmean). The 755 data points are from 10 MAN (1941–2019) and 14 IMIS stations (2000–2019), but not every station covers the entire time span)

3.6 Snowmelt throughout the twenty-first century

We projected future snowmelt dates by applying the linear mixed model to the CH2018 Climate Change Scenarios (CH2018 Project Team 2018). For the seven MAN stations with scenarios for temperature and precipitation, a HSmean of 100 cm on average corresponded to an accumulated winter precipitation of 560 ± 70 mm (minimum of 460 mm in 4MO, maximum of 650 mm at GRH) and disclosing no effect of elevation. The correlation of HSmean and the total winter precipitation for the period 1982 to 2019 was high, with Pearson correlation coefficients between 0.74 (DIS) and 0.86 (4UL). Only at the station Guetsch (GUE) it was with 0.65 lower.

Compared to the observed snowmelt dates, the predicted snowmelt dates for even years (model refitted with data of uneven years only) had a standard error of less than one day (0.82 days), highlighting a good model fit. With HSmean estimated by the station-specific proportionality factors, the standard error increased by 15% to 0.94 days.

For our model validation, we used the reference period of the climate change scenarios: 1990 to 2019 (MAN) and 2000 to 2019 (IMIS). At MAN stations, medians of the observed and simulated snowmelt dates deviated by zero to seven days only. The median snowmelt date of IMIS stations was also well reproduced by the model outputs. Deviations of one (TUJ2) to 14 days (OBW2) were overall larger than at MAN stations, because the model parameters extracted from a 2×2 km grid did not always mirror the exact station conditions (estimated vs. observed snowmelt in Online Resource 6). The time series analysis of the observed snowmelt dates revealed no significant trend during the reference period, except for one station (10 days decade–1 at TUJ3, p = 0.03). Most CH2018 simulations for the reference period did not display significant trends in the snowmelt date either. On average 11% of the simulations at MAN stations featured significant trends towards earlier snowmelt, and 7% at IMIS stations, but without uniform direction.

After the validation, we simulated the future snowmelt dates based on CH2018 projections for the beginning of the twenty-first century (2020–2049), for the mid-century (2045–2074), and for the end of the century (2070–2099) as well as the overall trend for the twenty-first century (2020–2099; Fig. 5). The trends of the simulated time series of snowmelt dates within each RCP were relatively concise regarding their direction. A 1000 m elevational gradient caused a decrease in median snowmelt trends by 0.5 days decade–1 only; thus, the magnitude of the effect of elevation was negligible (p = 0.17, R2 = 0.10). We here present the rates of change (trends) in snowmelt dates throughout the twenty-first century as the median trend across stations (7 MAN and 14 IMIS stations) and of a total of 68 simulations (12 for RCP2.6, 25 for RCP4.5, and 31 for RCP8.5). Extreme trends in snowmelt (outliers in Fig. 5) were clearly identifiable as results of single simulations and not because we aggregated stations of different elevations. As there were no large differences in trends between MAN and IMIS stations, we combined the results (Fig. 5), and we rounded our estimations to whole numbers (Table 3).

Fig. 5
figure 5

Projected trends in snowmelt dates of 21 (7 MAN and 14 IMIS) stations throughout the 21st century for the three emission scenarios RCP2.6, RCP4.5 and RCP8.5. In total, 68 simulations based on downscaled regional climate models: 12 for RCP2.6, 25 for RCP4.5, and 31 for RCP8.5, whereby the uneven sample sizes are responsible for the variability of trends between RCP scenarios in the early century. The coloured box contains 50% of all simulations per RCP scenario and time period. The region shaded in grey indicates the 25th and 75th quantile trend during the 30-year period 1990 to 2019

Table 3 Projected deviation of snowmelt dates compared to the reference period (MAN 1990–2019, IMIS 2000–2019; median, 25th and 75th quantiles, expressed as number of days: negative values mean earlier snowmelt, positive later snowmelt) for the future 30-year periods and three RCP scenarios. All values are rounded to whole days

In RCP2.6, predicted snowmelt dates did not show any change throughout the twenty-first century. Under RCP4.5, median snowmelt dates throughout the century will occur 2 days earlier each decade and 81% of the trends were statistically significant. This overall trend was mainly caused by a pronounced advancement during the beginning of the century (2020–2049, Fig. 5), with 3 days per decade. This shift dropped to 1 day per decade by the end of the century. For the whole century, the RCP4.5-models projected the median snowmelt date to advance by 10 days by the mid-century, and by 13 days by the end of the century (Table 3).

With unabated emissions (RCP8.5), snow will melt 4 days earlier each decade between 2020 and 2099 (trends of all simulations at all stations statistically significant). In contrast to the other two RCP scenarios, we found that snowmelt dates will accelerate over the century. While at the beginning, snowmelt will occur 2 days earlier per decade, towards the end of the century it will be 6 days earlier per decade. Median snowmelt dates will be 14 days earlier in the mid-century and 27 days earlier by the end of the century (Table 3).

To further illustrate the effect of these trends, we took a closer look at 2AN at 1440 m asl, a touristic skiing resort in central Switzerland, and at the two highest IMIS stations BED2 (2450 m asl) and LUK2 (2555 m asl). During the reference period (1990–2019), the median snowmelt date in 2AN occurred in late April, and only once before 1 April (2011). Under RCP4.5, snowmelt before April could become normal in one out of 12 winters, and under RCP8.5, this would be the case in one out of four years. At the IMIS stations BED2 and LUK2, the earliest snowmelt was at the beginning of June during 2000 to 2019, and the median snowmelt date was in the first third of July. In one out of four winters, snowmelt will occur in the first half of June (RCP4.5) or even before 1 June (RCP8.5).

4 Discussion

4.1 Past trends of snowmelt dates

Between 1958 and 2019, snowmelt occurred 2.8 ± 1.3 days earlier per decade. This was substantially more conservative than the 5.8 days per decade stated by Klein et al. (2016) for the period 1970 to 2015. We could associate these different rates to the specific time periods. Trends in snowmelt dates varied considerably throughout the past decades, as revealed by trends in 21-year running windows. The late 1980s and early 1990s were marked by substantially earlier snowmelt dates. This so-called regime shift manifested itself by earlier snowmelt and caused a 20% reduction of snow days above 1300 m asl in Switzerland (Marty 2008). Moreover, rapid climatic shifts during the late 1980s were evident in other datasets worldwide, for example, ocean and air temperatures, sea ice extent and grape ripening dates (Reid et al. 2016). In Switzerland, March to May temperatures increased by 0.39 K decade–1 between 1959 and 2008 (Ceppi et al. 2012), and by 0.84 K decade–1 in the period 1975 to 2004 (Rebetez and Reinhard 2008), suggesting that spring temperatures caused the regime shift in snowmelt dates described above. For Europe, this regime shift has recently been explained by natural, coincidental anomilies in the atmospheric circulation that were possibly related to the Pacific Decadal Oscillation and the Arctic Oscillation. This short-term warming superimposed a long-term trend of rising temperatures (Sippel et al. 2020). Locally, as for example in Switzerland, the decline of these atmospheric anomalies even led to a (short-term) trend towards lower temperatures (Ceppi et al. 2012; Bader and Fukutome 2015; Saffioti et al. 2016). Indeed, MeteoSwiss highlighted the dominant role of cyclonic low-pressure weather situations over Europe during winters (Dec-Feb) since the 1990s. Accordingly, we do not see any clear trends in snowmelt dates after the regime shift. Colder winters due to large-scale weather phenomena over Europe may thus have masked the recent warming and slowed the retreating snowmelt dates.

4.2 The role of elevation

Most studies on snow cover changes (e.g. Laternser and Schneebeli 2003; Scherrer et al. 2004; Marty 2008; Klein et al. 2016) focussed on stations at lower elevation, including only a few in the alpine vegetation belt (e.g. Gr. St. Bernard, Weissfluhjoch). Here, we included the relatively new snow depth series of IMIS stations and clearly demonstrated that the snowmelt dates correlated well with MAN stations, despite the two-months earlier snowmelt dates at lower elevations. We did also not find a pronounced elevational gradient in the trends of snowmelt dates. Scherrer et al. (2004) observed an only weak elevational gradient for the snow days (Dec–Feb) along an elevational gradient between 400 and 2500 m asl. Laternser and Schneebeli (2003) highlighted that the snow cover duration shortened fastest between 1000 and 1600 m asl, because temperature increases close to the zero-degree isotherm had the largest impact on the snow cover (Ceppi et al. 2012). The advancement of snowmelt dates at lower elevations may be slowed down by the reduced amount of solar radiation that is available when snowmelt dates recede towards earlier in the year (Musselman et al. 2017). At higher elevations, where the snowpack may persist beyond the mid-summer, receding snowmelt dates may move to a time with highest solar radiation and thermal energy and trends may therefore be accelerated (Essery et al. 2020).

Counter to our expectations, our proportionality factor describing the relation of winter precipitation to HSmean was unaffected by elevation. We would have assumed a temperature dependence, due to a higher fraction of liquid precipitation at lower elevations. However, under-catch and wind-induced snow redistribution at high elevations may potentially balance out the effect of an elevational gradient.

The RCP scenarios displayed slightly more moderate trends at higher elevations, but with less than 0.5 days decade–1 for 1000 m elevational difference, this was negligible.

4.3 Effect of temperature and precipitation on snowmelt dates

Temperature and precipitation, the two key factors for determining snowmelt dates, are highly affected by climate change (IPCC 2013). Morán-Tejeda et al. (2013) emphasised temperature as the driving factor for the snowpack duration below 1200 m asl. Above this elevation, accumulated precipitation and the winter-history of temperature and precipitation were more influential. For snowmelt, we found T21d ≥ 5 °C to be 1.5 times more influential than the precipitation proxy HSmean. Moreover, snowmelt dates related to the square-root of HSmean, revealing a decreasing influence of snow depth at higher elevations. Our measures for temperature and precipitation thus differed from monthly mean air temperatures and accumulated precipitation considered in Morán-Tejeda et al. (2013). T21d ≥ 5 °C represents a threshold when temperatures are high enough for the whole snowpack to melt. Although HSmean correlates with accumulated winter precipitation, it does not always directly relate to monthly precipitation. Both our measures account for the temperature and precipitation history during winter.

4.4 Snowmelt dates throughout the twenty-first century

Accelerating snowmelt dates above the climatic treeline can only be halted by abating emissions enormously, as our simulations predict a roughly one-month earlier snowmelt for RCP8.5 and two-week earlier snowmelt for RCP4.5 by the end of the twenty-first century. Compared to previous projections of snow cover changes, our estimations are rather conservative. Bavay et al. (2009) projected snowmelt dates in eastern Switzerland to occur 40 days earlier by the end of the century, and in the large Aare catchment (3190 km2) the snow duration at 2000 m asl may shorten by 2.5 months compared to the reference period (1999–2012; Marty et al. 2017). Such early melt could be interpreted as a direct consequence of lower snow depth and snow water equivalents, as simulated by Schmucki et al. (2015b) and Steger et al. (2013). In contrast to these studies, our model is fully based on empirical observations and not on physical processes. Our data used for model fitting covered a period with massive changes in snowmelt dates; nevertheless, the relation of T21d ≥ 5 °C and HSmean to the snowmelt date was robust.

4.5 Consequences for alpine plants

Our simulation results also enable us to put alpine snow manipulation experiments into the context of climate change scenarios. Some experiments realised 18-day earlier snowmelt dates by reducing the snow depth in spring, summarised by Wipf and Rixen (2010). Under RCP8.5, these experiments would thus reproduce conditions we may observe in the mid-century (2045–2074).

A one-month earlier snowmelt may expand the alpine growing season by more than one third by the end of the century (RCP8.5), most likely with drastic consequences for alpine plants. These may include lower flower numbers, reduced leaf growth and a lower survival, as observed for the evergreen forb Gentiana nipponica by a one-month earlier snowmelt (Kawai and Kudo 2018), and snowbed species produced less seeds after an earlier snowmelt of 18 days (Tonin et al. 2019). According to an experiment with 23 nival and alpine species, it is assumed that one third to one half the alpine flora may be sensitive to photoperiod and hence may not profit from a longer growing season because of too short daylengths (Keller and Körner 2003). The 4h-difference in daylength they applied in the greenhouse study is comparable to snowmelt in mid-March instead of late June, a shift of roughly three months. Such a shift exposes plants to higher frost risks. Recurrent frosts killed flower buds in alpine herbs (Inouye 2008), impaired growth snowbed species (Baptist et al. 2010b) as well as in three dwarf shrub species (Wipf et al. 2009). Francon et al. (2020) further suggested a strong elevational impact, disclosed by growth ring analysis of the shrub Rhododendron ferrugineum. Longer growing seasons at 2400 m asl enhanced its growth, while at 1800 and 2000 m asl early snowmelt dates led to frosts and ultimately limited growth.

Because a one-month earlier snowmelt will not reduce the daylength sufficiently for photoperiod restrictions, we assume that in the long-term, frost events at the beginning of the growing season will play a bigger role than photoperiod limitations and may diminish the effect of longer growing seasons.

5 Conclusions

We show a high correlation between snowmelt dates at high-elevation IMIS stations and lower-elevation MAN stations, implying that snowmelt dates are advancing with similar rates along an elevational gradient of 1000 m asl to 2500 m asl in the Swiss Alps. Our empirical approach highlighted the dominant effect of temperature over a slightly weaker influence of mean snow depth (HSmean) on snowmelt dates. With T21d ≥ 5 °C, we found an appropriate temperature measure that correlated well with snowmelt dates at different elevations. By applying both, T21d ≥ 5 °C and HSmean to CH2018 Climate Change Scenarios, we could anticipate snowmelt dates throughout the twenty-first century.

By the end of the century, we expect an earlier snowmelt by up to one month. Such an early melt will affect hydropower production, winter tourism in the Swiss Alps and will expose the alpine flora and soils to a higher frost risk with additional consequences on biogeochemical cycles.