Snow depth sensitivity to mean temperature, precipitation, and elevation in the Austrian and Swiss Alps

Snow depth plays an important role in the seasonal climatic and hydrological cycles of alpine regions. Previous studies have shown predominantly decreasing trends in average seasonal snow depth across the European Alps. Additionally, prior work has shown bivariate statistical relationships between average seasonal snow depth and mean air temperature or precipitation. Building upon existing research, our study uses observational records of in situ station data across Austria and Switzerland to better quantify the sensitivity of historical changes in seasonal snow depth through a multivariate framework that depends on elevation, mean temperature, and precipitation. These historical sensitivities, which are obtained over the 1901–1902 to 1970–1971 period, are then used to estimate snow depths over the more recent period of 1971–1972 to 2020–2021. We find that the year-to-year estimates of snow depths, which are derived from an empirical–statistical model (SnowSens), that rely solely on the historical sensitivities are nearly as skillful as the operational SNOWGRID-CL model used by the weather service at GeoSphere Austria. Furthermore, observed long-term changes over the last 50 years are in better agreement with SnowSens than with SNOWGRID-CL. These results indicate that historical sensitivities between snow depth, temperature, and precipitation are quite robust over decadal-length scales of time, and they can be used effectively to translate expected long-term changes in temperature and precipitation into changes in seasonal snow depth.

1 Introduction

Snow on the ground is an important component of the hydrological cycle, the climate system, and mountain ecosystems throughout the world (Beaumet et al., 2021; Beniston et al., 2018; Gobiet et al., 2014; Notarnicola, 2022). The timing of snowfall, along with its accumulation, has profound implications for water resources (Viviroli et al., 2011; Colombo et al., 2023; Avanzi et al., 2024), mountain tourism (Elsasser and Bürki, 2002), and mountain hazards such as avalanches (Marty et al., 2017b). Understanding the impact that climate change has on this valuable resource is therefore essential to better assist in regional planning and preparedness.

Correctly quantifying changes in the climate system or hydrologic cycle generally requires robust measurements with sufficiently long time series of high data quality. Historically, there are two quantities that have been measured in situ by national hydrometeorological services that fulfill these criteria: (1) snow depth and (2) depth of snowfall. Depth of snowfall is defined as freshly fallen snow that accumulated on a snow board during a standard observing period of 24 h, while snow depth is the total accumulated snowpack from the ground surface to the snow top (Haberkorn, 2019). In the European Alps, there is a long history of snow depth and snowfall measurements that date back to the 19th century (Scherrer et al., 2013). Over the more recent past, satellite data have become increasingly important in providing information on specific spatial patterns (Hüsler et al., 2014; Lievens et al., 2019). However, the measurement record of satellite data is relatively short and unfortunately cannot provide the same quality of information as in situ observations when it comes to quantities involving snow depth (Lievens et al., 2022). Other variables describing snow characteristics, such as snow water equivalent (SWE), were introduced later in the measurement record and with a lower-density network (see, e.g., Haberkorn, 2019).

Matiu et al. (2021) conducted one of the first studies to provide an extensive and comprehensive analysis of changes in snow depth for the period of 1961–2020 that truly covers the entire region of the Alps. It is worth mentioning their great effort in merging many individual stations across different institutions and networks. In their study, Matiu et al. (2021) were able to show, predominantly, a decreasing trend in snow depth across the Alps. In addition to the regional differences that they found in the trends, they also showed a strong elevation dependence of snow depth trends. However, this elevation dependence of the snow depth trends is conflated, to some degree, with the fact that stations at higher elevations also typically receive more snow. As a result, there is also a benefit in investigating whether the relative changes in snow depth are increasing or decreasing with elevation (Laternser and Schneebeli, 2003; Marty and Blanchet, 2012; Marty et al., 2023).

Beyond quantifying historical trends in snow depth and snowfall themselves (Bertoldi et al., 2023), it is additionally useful to attribute these changes to a certain set of physical drivers. The accumulation of snow depth over a season is primarily driven by temperature and precipitation (Sippel et al., 2020; Pepin et al., 2022). There have been several prior studies that have linked changes in snow depth, at different elevations, to changes in air temperature and precipitation (Scherrer and Appenzeller, 2006; Morán-Tejeda et al., 2013; Sospedra-Alfonso et al., 2015; Scalzitti et al., 2016; Schöner et al., 2019; Monteiro and Morin, 2023). Overall, these studies have shown snow depth to be strongly related to air temperature at low elevations and to precipitation at high elevations (Morán-Tejeda et al., 2013; Schöner et al., 2019). However, these studies also experience, to some extent, a conflation between snow depth quantities and elevation. Furthermore, the statistical relationships shown are often correlations, which do not capture how much snow depth would change, for example, as a function of air temperature.

Our study aims to extend prior work in a number of ways. First, we start by computing anomalies of snow depth, temperature, and precipitation data for stations across Austria and Switzerland. This step removes regional and elevation-dependent climatological differences, thereby allowing us to better quantify anomalous or relative changes. Next, we observe how sensitive snow depth has been to anomalous changes in temperature and precipitation in the historical record for stations within specified elevation bands. Then, we use the historically derived sensitivities to construct an empirical–statistical model to estimate seasonal snow depth based on seasonal anomalies of temperature and precipitation. Lastly, we evaluate model performance. The model is calibrated over the period of 1901–1902 to 1970–1971, and evaluation is performed over the period of 1971–1972 to 2020–2021. The model evaluation is also compared to that of GeoSphere's SNOWGRID-CL model for the Austrian domain. Our primary objective is to provide an effective yet easy-to-interpret method to translate expected long-term changes in temperature and precipitation to changes in seasonal snow depth.

2 Data

Stations with daily measured snow depth were collated from the Austrian and Swiss meteorological services GeoSphere Austria, Hydrographischer Dienst (HD), MeteoSwiss, and the WSL Institute for Snow and Avalanche Research (SLF). Some of these stations began record keeping in the 1880s, while many others became active in the early part of the 20th century. At the time that the authors collected the data, these stations primarily had data coverage through the spring of 2021. While seasonal values of snow depth, mean temperature, and precipitation reflect the accumulations or averages spanning from one year to the next, for the duration of the paper, we simply use the year in which the season ends. So, for example, the year 2021 would refer to the season of November 2020–March 2021. There are a total of 291 snow stations, with 107 stations in Austria and 184 stations in Switzerland. The stations range in elevation between 121 and 2536 m, with a mean height of 1097 m (see Fig. 1). While we use the non-homogenized snow depth measurements in this study, it is not expected that the results would substantially change using homogenized snow data. Through personal communication, the authors of a recent homogenization study in the Alps (i.e., Resch et al., 2022) indicated that there are not any systematic changes in snow depth one way or the other as a result of the homogenization procedure (Marcolini et al., 2019; Buchmann et al., 2022). Furthermore, the primary focus of this study is to present a useful methodology or approach to quantify the influence that anomalous seasonal temperatures and precipitation have on snow depth.

Figure 1 Study region with elevation. Stations with historical measurements of daily snow depth are plotted as the magenta circles. Monthly homogenized temperature and precipitation stations are plotted as the hollow orange boxes and the filled blue squares, respectively.

Monthly homogenized temperature and precipitation data for Austria are obtained from GeoSphere (previously ZAMG) as part of the HISTALP dataset (i.e., Historical Instrumentation Climatological Surface Time Series of the Greater Alpine Region, https://www.zamg.ac.at/histalp/dataset/station/csv.php, last access: 5 July 2023), while homogenized precipitation and temperature data for Switzerland are obtained from MeteoSwiss (https://www.meteoschweiz.admin.ch/service-und-publikationen/applikationen/ext/climate-tables-homogenized.html, last access: 5 July 2023). There are a total of 43 temperature stations and 48 precipitation stations in Austria (see Fig. 1), while, in Switzerland, there are a total of 29 temperature stations and 27 precipitation stations.

We evaluate the statistical model across Austria against a dynamical snow cover model, SNOWGRID-CL (Olefs et al., 2020), which is a simplified climate version of the operational snow model SNOWGRID (Olefs et al., 2013) of GeoSphere Austria. The model was developed for climatological simulations such as long historical runs and future scenarios. It relies on an extended degree-day scheme to approximate snow ablation from air temperature and the shortwave radiation balance (see Olefs et al., 2020, for further details). The model is forced with an observation-based gridded dataset of air temperature and precipitation (SPARTACUS v2.1, Hiebl and Frei, 2016, 2018) and simulates daily fields of snowpack properties (i.e., snow depth, SWE) at a spatial resolution of 1 km×1 km over the Austrian domain. The model output is updated daily and is publicly available on GeoSphere Austria's data hub (https://public.hub.geosphere.at/public/datahub.html?id=snowgrid_cl-v2-1d-1km/filelisting&anonymous=true#/snow_depth/, last access: 1 July 2023).

Evaluation metrics

The performances of different modeled time series of snow depth are compared using the root mean squared error (RMSE) statistic. This metric is used because it measures how well the estimates from the model covary with observations, while it additionally reflects whether or not there is any systematic mean bias between the two time series. This metric is computed using time series of both absolute (which are the raw modeled and observed values expressed in centimeters and should not be confused with the mathematical absolute value) and anomalous values. The modeled RMSE is calculated as follows:

$$\begin{array}{}\text{(1)}& {\mathbf{RMSE}}_{\text{MOD}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n(y_{\text{mod},i}-y_{\text{obs},i}{)}^{\mathrm{2}}},\end{array}$$

where y_mod,i and y_obs,i are, respectively, seasonal time series of modeled and observed snow depths at station i. Likewise, RMSE_CLIM is defined as

$$\begin{array}{}\text{(2)}& {\mathbf{RMSE}}_{\text{CLIM}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n(y_{\text{clim},i}-y_{\text{obs},i}{)}^{\mathrm{2}}},\end{array}$$

reflecting the error associated with estimates using climatology, where y_clim,i is the reference mean climatological snow depth at station i. One can think about y_clim,i as either a single value (i.e., the climatological mean) or a time series array, with all of the values being the same. The values of y_mod,i and y_clim,i change depending on whether we are computing the RMSE skill score using absolute values or anomalies. For the absolute values evaluated over the period of 1972–2021, y_obs,i, y_mod,i, and y_clim,i would all contain values expressed in centimeters, where y_clim,i is the mean seasonal snow depths, computed station by station, over the 1902–1971 calibration period. When computing the skill of the model estimated anomalies, y_obs,i, y_mod,i, and y_clim,i all contain values expressed as a percentage of normal (e.g., 120 % of normal, which is 20 % above normal), where y_clim,i is the mean seasonal snow depth at station i over the 1902–1971 calibration period (i.e., 100 %). Additionally, an RMSE skill score is also used to evaluate the performance of the models in terms of their ability to capture observed trends. In that case, a trend of 0 % per decade is treated as the climatological reference. Then, the RMSE skill score, SS_RMSE, is defined as follows:

$$\begin{array}{}\text{(3)}& {\mathbf{SS}}_{\text{RMSE}}=\mathrm{1}-{\displaystyle \frac{{\mathbf{RMSE}}_{\text{MOD}}}{{\mathbf{RMSE}}_{\text{CLIM}}}},\end{array}$$

where an SS_RMSE value of 1.0 would perfectly estimate the observations, while values between 0.0 and 1.0 reflect model estimates that perform better than the climatology, and values below 0.0 indicate that the model is less skillful than the climatology.

3 Methods

3.1 Seasonal snow depth

In this paper, we focus on snow depth averaged over the November–March season. Our first goal is to investigate historical empirical relationships between mean seasonal temperature, precipitation, and snow depth. Second, we use these historically derived empirical relationships to estimate changes in snow depth driven by changes in mean seasonal temperature and precipitation. During warmer months, and particularly at stations at lower elevations, an observable amount of precipitation will not always translate to a measured snow depth. This would result in trying to fit a predictor time series (i.e., precipitation) that varies with our predictand time series that does not (i.e., snow depth). We try to minimize the number of cases where there is zero measured snow depth by excluding the months of April and May from our seasonal average since those months often contain more stations with zero recorded snow depth. That way, we can have consistency in the lengths of our seasons for both the predictors, temperature (T) and precipitation (P), and our predictand, snow depth (SD). Figure 2a shows the percentage of data for which the 291 stations in this study measured an average monthly snow depth greater than 0.0 cm. For example, when considering all of the Januaries between 1901–2020, there were 20 263 station months with full data coverage. Of those, 500 recorded 0.0 cm for every day throughout January at a station. The percentage of zero measured snow depth for January is then 2.5 %, which is equal to 500$/$20 263 (this quantity is also equal to 100 % minus the percentage shown in the left-most bar of Fig. 2a). We used the season November–March because each of those 5 months contained less than 20 % zero measured snow depth. While the November–March season is somewhat shorter than what some other studies have used (e.g., Matiu et al., 2021; Morán-Tejeda et al., 2013), the average anomalous November–March snow depth (anomalies are represented as a percentage of normal; see Eq. 6) varies strongly with the average anomalous snow depth over the longer November–May season (see Fig. 2b). Thus, if one can skillfully estimate November–March snow depth then these will also be skillful for a longer season such as November–May. In Fig. 2c, average seasonal snow depth can be seen to vary with elevation, with higher elevations generally receiving more snow.

Figure 2 Panel (a) plots the percentage of data where the monthly average snow depth was greater than 0.0 cm. Panel (b) plots November–March average snow depth as a percentage of normal over the period of 1902–2021 against the longer November–May season. Panel (c) shows how the average seasonal snow depth varies with elevation, with higher elevations generally receiving more snow. Note the logarithmic scaling of the y axis in (c).

Observed changes in seasonal snow depth

Figure 3 shows the trends of snow depth anomalies (see Eq. 6) across the Austrian and Swiss Alps for the period of 1902–2021 using our four elevation bands. Here, we have computed linear trends. However, we should mention a couple of small caveats in doing this. First, the data are bounded by zero, and a negative trend line would eventually cross the origin to produce negative values of snow depth, which cannot physically happen. Second, the data do typically exhibit some level of positive skewness, and, as a result, the data cannot be considered to be truly Gaussian. In order to allow for easily interpretable trend information, linear trends are computed over the time series. That way, even though the values are bounded by zero and also cannot be considered to be truly Gaussian, it can be illustrated that values have decreased over time, i.e., show a negative trend. Additionally, we primarily want to illustrate that trends have decreased across each elevation band and that the trends have been greatest at lower elevations. For the stations below 500 m, the decrease in snow depth is −8.1 % per decade (Fig. 3a), while the stations at elevations above 1500 m exhibit less than half of the relative trend at −3.4 % per decade (Fig. 3d). Additionally, the regime shift at the end of the 1980s described by Marty (2008) for the Alps and by Reid et al. (2016) for the global level is nicely visible, especially for the two lower elevation bands. Please refer to Appendix A to see trends over a more recent period, as well as the spatial patterns of the historical changes.

Figure 3 The seasonal (November–March) snow depth anomalies for all stations and all seasons are plotted as the dots for each of the four different elevation bands. Trend lines and the relative change per decade are also plotted. Here, we use the entire period of record. All station–season values are used to compute the historical trend. For example, the trend in (a) is the least-squares regression fit to all of the points in (a). Bootstrapping is used to estimate 95 % confidence intervals, which are shown as the blue shading.

Different snow depth stations have different record lengths. As a result, the anomalies from the estimated mean seasonal snow depth can vary, to some degree, depending on which seasons are used to compute the average. In Table 1, we show a range of historical trends depending on how much data coverage we set as a threshold. The trends are not found to change very much when only the stations with more complete data coverage are used.

Table 1 Table values express trends and changes across different time periods for the four elevation bands. The trends are percentage changes per decade, given for the two different time periods from Figs. 3 and A1 in the Appendix, while the changes are the average percentage changes between two periods over the last 70 years. Stations have varying data lengths and coverages, which can influence mean estimation. The percentage of data coverage is evaluated over the 1951–1952 to 2020–2021 period, where 80 % coverage would mean that a particular station had at least 56 years of data. The first three sets of rows show how much the trends change when using stations with more or less data coverage. The bottom set of rows gives the number of stations used to compute the trends and changes. The units of the rows are provided within the square brackets.

3.2 Constructing homogenized temperature and precipitation time series at the snow stations

By itself, snow depth is not a measure of how much melted water is contained in the snow mass. Many different meteorological conditions can affect the density of the snowpack or the melting of snow. As a result, we cannot directly infer how much precipitation has fallen based on the snow depth measurements themselves. The precipitation could have fallen as rain, while precipitation falling as snow can accumulate at varying densities. Ultimately, we want to quantify historical changes in mean temperature and precipitation and how these have translated into changes in snow depth. Therefore, since many of the snow depth stations have neither temperature nor precipitation measurements, we construct these time series using monthly homogenized values of temperature and precipitation using nearby stations over the period of 1901–2021.

First, we obtain November–March sums of precipitation and averages of mean temperatures at all of the homogenized stations (see Fig. 1) over the years 1901–1902 to 2020–2021. Additionally, time series of standardized anomalies (i.e., z scores) are computed for each homogenized station (the mean and standard deviations are computed using the calibration period of 1902–1971). Then, for each snow depth station, we find the nearest five homogenized stations for temperature and, separately, for precipitation. Next, we compute two time series. The first contains the averaged standardized anomalies of temperature or precipitation as a function of the inverse distance of the nearest five stations. The second time series contains the averaged absolute values of temperature or precipitation as a function of the inverse distance of the nearest five stations. The inverse-distance-weighted standardized anomalies of the first time series are then adjusted to match the mean and standard deviation of the second time series. This is performed by simply reversing the steps of computing a z score, where the weighted standardized anomalies at each station are first multiplied by the standard deviation of the second time series at the same station, following which we add the mean of the second time series (also at the same station). This provides us with time series of seasonal (i.e., November–March) mean temperature and precipitation located at each of our 291 snow depth stations.

3.3 Sensitivity of snow depth to temperature and precipitation

Information concerning the bivariate correlations between mean temperature and snow depth and/or precipitation and snow depth can be useful. However, correlations by themselves do not provide information about the steepness of the slope between the two variables. For example, given a 1.0 °C increase in mean temperature or a 20 % increase in precipitation, what would the expected impact on snow depth be? Additionally, how might these expected changes be affected by elevation? Furthermore, what would the multivariate impact on snow depth be given some combination of mean temperature and precipitation changes? And finally, can we apply methods that do not have an underlying assumption of linearity?

In Fig. 4a and b, one can observe the spatial distribution of the bivariate correlations between snow depth and temperature and, similarly, the correlations between snow depth and precipitation. In Fig. 4c and d, these correlations are shown to vary as a function of elevation. Generally, we find the largest correlations (either positive or negative) at lower elevations for temperature and at higher elevations for precipitation. For example, a station below 500 m is more likely to see increases or decreases in temperature translate more strongly to decreases or increases in snow depth than stations at higher elevations. The opposite influence is observed for snow depth and precipitation. Therefore, stations (or regions) at lower elevations are primarily driven by changes in temperature, while stations at higher elevations are primarily driven by changes in precipitation. These relationships support prior findings, such as those by Morán-Tejeda et al. (2013) and Schöner et al. (2019).

Figure 4 The sizes and magnitudes of the Pearson correlation coefficients between NDJFM snow depth and temperature are plotted in panel (a). Panel (b) shows the same as (a) but using the seasonal precipitation time series instead of temperature. Downward- and upward-facing triangles reflect negative and positive correlations, respectively. The sizes of the triangles reflect the magnitude of the correlation. Panel (c) plots the relationship between station elevations and the correlations between snow depth and temperature (i.e., the same values as in a), while (d) plots the relationship between station elevations and the correlations between snow depth and precipitation (same values as in b).

The utility of the information from Fig. 4 can be improved in the following ways: (1) instead of only considering the bivariate statistical relationship between either temperature with snow depth or precipitation with snow depth, we can consider all three variables in a non-linear, multivariate framework, and (2) by computing anomalies in the data, we can leverage information across multiple stations to provide a more robust empirical–statistical relationship.

In Fig. 5, we illustrate why it can be important to use anomalous values with our methodology. Figure 5a–c shows bar plots for the distribution of values using absolute temperature, precipitation, and snow depth, respectively, for the Austrian and Swiss stations between 500–1000 m. The average station elevation of the Austrian stations used is 745 m, while it is 742 m for the Swiss stations. Hence, the average station elevations are not very different between the two regions. One can observe in Fig. 5a and b that the Swiss stations are generally warmer and wetter than their Austrian counterparts. At the same time, the Swiss stations have lower seasonal averages of snow depth. In order to stress our point regarding the usefulness of anomalies versus absolute values in our methodological context, we can take a further subset of these Austrian and Swiss data points over this 500–1000 m elevation band. The individual values of these subsets of data are shown as the scatterplots in Fig. 5d and e. A Student's t test shows that the means (for temperature, precipitation, and snow depth) of the subset of Austrian data points (Fig. 5d) are all statistically significantly different compared to the subset of Swiss data points (Fig. 5e). We find that, while this subset of data points in Austria has a greater absolute temperature and less absolute precipitation than the Swiss subset, the Austrian stations have significantly more absolute snow depth than the Swiss stations. As we decrease temperature and increase precipitation, we should expect snow depth to increase. However, this is exactly the opposite of what the absolute data are telling us. By simply using the absolute data alone, we can potentially get the wrong signal when comparing one region to another. This spatial climatological difference can be addressed through computing anomalies, station by station, for the different variables. After computing anomalies, we can then leverage the information across a larger region.

Figure 5 Panels (a–c) show bar plots for the distribution of absolute temperature, precipitation, and snow depth, respectively, for the Austrian and Swiss stations between 500–1000 m over the historical period of 1902–1971. The percentages of the blue and orange bars in each panel (a–c) sum up to 100 %. The bar plots are comprised of 1755 observed data points for Austria and 558 data points for Switzerland. A subset of these Austrian and Swiss data points are shown in the scatterplots in (d) and (e), respectively. These are historical, observed individual seasonal measurements that fall within the climatological regions bounded by the black rectangles. The size of the colored squares reflects the values of absolute snow depth, where larger squares correspond to larger snow depths. The average values of the Austrian and Swiss subsets of the data are provided in the inset text in (d) and (e).

To compute anomalies in the data, we begin with absolute (or raw) values of November–March seasonal temperature, precipitation, and snow depth data. Then, November–March average temperature anomalies, ${\mathbf{T}}_{x,t}^{\ast }$, for station $x\in (\mathrm{1},\mathrm{\dots },\mathrm{291})$ and year $t\in (\mathrm{1902},\mathrm{\dots },\mathrm{2021})$ are computed as

$$\begin{array}{}\text{(4)}& {\mathbf{T}}_{x,t}^{\ast }={\mathbf{T}}_{x,t}-{\stackrel{\mathrm{‾}}{\mathbf{T}}}_x,\end{array}$$

where T_x,t denotes the seasonal averages of absolute temperature, and ${\stackrel{\mathrm{‾}}{\mathbf{T}}}_{\mathbf{x}}$ is the time-averaged mean temperature over the calibration period of 1902–1971 at station x. Anomalies of November–March precipitation accumulations, P^∗, are computed as

$$\begin{array}{}\text{(5)}& {\mathbf{P}}_{x,t}^{\ast }={\displaystyle \frac{{\mathbf{P}}_{x,t}}{{\stackrel{\mathrm{‾}}{\mathbf{P}}}_x}},\end{array}$$

where P_x,t denotes the seasonal accumulations of absolute precipitation, and ${\stackrel{\mathrm{‾}}{\mathbf{P}}}_{\mathbf{x}}$ is the time-averaged precipitation over the period of 1902–1971 at station x. Similarly, anomalies of November–March average snow depths, HS^∗, are obtained using

$$\begin{array}{}\text{(6)}& {\mathbf{HS}}_{x,t}^{\ast }={\displaystyle \frac{{\mathbf{HS}}_{x,t}}{{\stackrel{\mathrm{‾}}{\mathbf{HS}}}_x}},\end{array}$$

where HS_x,t denotes the seasonal averages of absolute snow depth, and ${\stackrel{\mathrm{‾}}{\mathbf{HS}}}_{\mathbf{x}}$ is the time-averaged snow depth over the period of 1902–1971 at station x.

Once we have computed our data anomalies, we can plot in Fig. 6 the observed historical anomalous snow depths (i.e., HS^∗) along with temperature and precipitation anomalies (i.e., T^∗ and P^∗, respectively). The larger squares, bounded by the black lines, correspond to the anomalous measurements for one example station. That station, named “Feldkirch” with the number “11110”, has the coordinates of 47.27 ° latitude and 9.60 ° longitude and is situated at an elevation of 439 m. As one would expect, the average snow depth anomalies increase as the temperature anomaly decreases and as the precipitation anomaly increases. This figure can also be used to show that, for some increase or decrease in temperature, there can be a corresponding increase or decrease in precipitation that will yield approximately the same snow depth anomaly. Consider, for example, the snow depth anomalies in Fig. 6 between 2–3 °C below normal temperature and between 100 %–150 % of normal precipitation. Similar snow depth anomalies can also be observed between 0–1 ° below normal temperature and between 200 %–250 % of normal precipitation. This gives us a general idea of how sensitive snow depth anomalies are to temperature and precipitation anomalies. Even though Fig. 6 gives us a first look at the multivariate sensitivity of snow depth to temperature and precipitation anomalies, we can refine the approach by adding in a third variable. The sensitivities will change as elevation changes. Therefore, we can break up the multivariate sensitivities shown in Fig. 6 into different elevation bands. We have chosen to use four elevation bands: 0–500 m (containing 52 stations), 500–1000 m (75 stations), 1000–1500 m (91 stations), and >1500 m (73 stations).

Figure 6 Seasonal temperature anomalies for all stations and all seasons are plotted on the x axis against seasonal precipitation anomalies on the y axis. The color map corresponds to snow depth anomalies given the pairings of temperature and precipitation anomalies. Open black squares show the values for one example station, which is named “Feldkirch” with the number “11110” and which is located at the coordinates 47.27 ° latitude and 9.60 ° longitude. These anomalies are shown for the entire 1902–2021 period of record.

The points along the left column in Fig. 7 are like those from Fig. 6, except that the data are broken up by the four elevation bands, and the data are now only plotted for the calibration period of 1902–1971. These are the data that we will use to fit a model and to make model estimates for our 1972–2021 validation period.

Figure 7 Panels (a–d) are the same as in Fig. 6 except for the fact that the data are broken up by elevation band, and we only plot the data over the 1902–1971 calibration period. Panels (e–h) show the binned anomalies of snow depth using two-dimensional binning windows. The bins are centered using a 0.8 °C temperature binning window (with 0.2 °C increments of the binning window) and a 40 % precipitation window (with 10 % increments). Panels (i–l) show the fitted surface through localized linear regression using the nearest quartile of bins from (e–h) (see text). The dashed black lines show contours of 100 % of normal snow depth.

To begin, we calculate the averages of snow depth anomalies across two-dimensional bins of mean temperature and precipitation for different elevation bands. In our effort to construct the sensitivity diagrams, let us first consider the elevation band 0–500 m. In the lowest elevation band over the calibration period, the arrays T^∗, P^∗, and HS^∗ all have a maximum possible number of measurements of 3640 (which is 70 seasons over 1902–1971 multiplied by the 52 stations from that band). Next, we find all of the values in ${\mathbf{HS}}_{\text{0–500\hspace{0.17em}m}}^{\ast }$ that fall within a 0.8 °C window centered about a given temperature anomaly (with 0.2 °C increments of the binning window), ${\mathbf{T}}_{\text{0–500\hspace{0.17em}m}}^{\ast }$, and a 40 % window centered about a given precipitation anomaly (with 10 % increments), ${\mathbf{P}}_{\text{0–500\hspace{0.17em}m}}^{\ast }$. We did experiment with using different bin sizes, but we found that the choice of bin size does not strongly affect model performance (please refer to our comments in the paper discussion at https://doi.org/10.5194/egusphere-2024-1172-AC2). We calculate the average of all of the station snow depth anomalies that fall within this two-dimensional window of temperature and precipitation anomalies, given that there were at least 50 observed snow depth measurements that fall within that two-dimensional window. Then, we move the center of the window in order to perform the same set of operations across the range of temperature and precipitation anomalies. And lastly, we repeat the process for the other three elevation bands. The resulting averages using these two-dimensional bins are shown in Fig. 7e–h.

Next, the snow depth anomalies, resulting from the multivariate binning (Fig. 7e–h), are used to construct a smoothed and extrapolated sensitivity surface that will then allow us to estimate new values outside of what has been seen in the calibration record. For all binned values across the two-dimensional temperature and precipitation space, distances (in data space) are computed between each specific bin center and the center of all of the bins where we computed averages in the previous section (i.e., those are the grid cells which are colored in Fig. 7e–h). Multiple linear regression is then fit to the nearest quartile of values, where the bin centers of mean temperature and precipitation anomalies are the predictors and where the binned averages of snow depth anomalies are our predictands. To find the nearest quartile of values, we use a Euclidean distance measure which essentially equates the distances of a 10 % precipitation anomaly to a 0.2 °C temperature anomaly. Thus, a data point that is 0.4 °C warmer with 0 % of normal precipitation with respect to a point of interest and another data point that is 0.0 °C with 20 % of normal precipitation would be treated as the same distance. We did not find the model to be overly sensitive to providing more or less weight to the temperature or precipitation axes. Then, the regression coefficients are used along with the center point of the bin to obtain a snow depth anomaly. Our application of localized linear regression is only fit to the nearest quartile of data points for each bin center; therefore, it can accommodate a non-linear response surface across most of the data domain. At the same time, it smooths out the sensitivity surface (Fig. 7i–l) and also provides extrapolated values beyond what was observed in the calibration period.

3.4 A multivariate sensitivity model to estimate snow depth

The multivariate sensitivity plots shown in Fig. 7i–l, which are constructed using only data in the 1902–1971 calibration period, are now used to estimate snow depth anomalies at all stations and for all NDJFM seasons over the period of 1972–2021. For each station, we first determine the elevation band in which it falls (e.g., 0–500 m, 500–1000 m). Second, we locate the nearest bin center given each season's (e.g., November 1971–March 1972) anomalies of mean temperature and precipitation. The snow depth anomaly corresponding to that bin is then used as the estimated value for that season for that station. For example, consider a station in the elevation band 0–500 m which experienced a mean seasonal temperature at 2.0 °C above normal (where normal is defined with respect to 1902–1971) and a mean seasonal precipitation at 100 % of normal. The snow depth anomaly in Fig. 7i at 2.0 °C along the x axis and 100 % of normal along the y axis corresponds to 37.6 % (indicated by the color in that bin). So, given those meteorological anomalies at that station and in that season, the historical sensitivities would lead to an estimation of the snow depth being 37.6 % of normal. This procedure is repeated to produce the model estimates for all 291 stations and for all 50 seasons over the 1972–2021 validation period. The model estimates derived in this manner and driven by the historical sensitivities are referred to as the SnowSens model for the duration of the paper.

To be clear, we use known values of concurrent seasonal temperature and precipitation anomalies to estimate snow depth using the SnowSens model. As a result, our modeled values are not really forecasts. While we refer to the values produced by the SnowSens model as “estimates”, our broader aim is for the model to be capable of producing actual forecasts of future snow depth conditions given future projections of temperature and precipitation. True long-term forecasts of snow depth over some future period can be thought of as a modeling chain with two chain links. The first link in the modeling chain is made up of the projections of temperature and precipitation, while the second link is made up of the projections of snow depth driven by the temperature and precipitation projections from the first link. In this paper, our model validation is done by quantifying the skillfullness or the uncertainty associated with the second link of that modeling chain. With knowledge of how skillful the SnowSens model is, the model can then be forced or driven with a plausible range of future projections of temperature and precipitation. In this way, the SnowSens model can be run to produce actual forecasts of snow depth, which also contain the uncertainty of future projections from the first link of the modeling chain.

When evaluating the performance of the SnowSens and the SNOWGRID-CL models, we use both the absolute and anomalous values. The methodology outlined above provides SnowSens model estimates as anomalies for all stations and seasons. At the same time, the SNOWGRID-CL model provides estimates of absolute snow depths. Therefore, we must produce absolute estimated values for the SnowSens model and estimated anomalies for the SNOWGRID-CL model. Let us begin with the estimates of the SnowSens model, where the absolute snow depths, HS_MOD, are computed as

$$\begin{array}{}\text{(7)}& {\mathbf{HS}}_{\text{MOD},x,t}={\mathbf{HS}}_{\text{MOD},x,t}^{\ast }\times {\stackrel{\mathrm{‾}}{\mathbf{HS}}}_{\text{OBS},x},\end{array}$$

where ${\mathbf{HS}}_{\text{MOD},x,t}^{\ast }$ denotes the SnowSens- estimated anomalies at station x and time t, and ${\stackrel{\mathrm{‾}}{\mathbf{HS}}}_{\text{OBS},x}$ is the observed mean seasonal snow depth at station x over the calibration period of 1902–1971. Now, we have absolute estimates for both models. Next, we want to address any mean biases present in the models while simultaneously computing anomalies in the data. To do this, a common period of record is used. The SNOWGRID-CL model data begin in 1962, and, therefore, we use the common reference period of 1962–1971 to compute the anomalies and bias correction:

$$\begin{array}{}\text{(8)}& {\mathbf{HS}}_{{\text{MOD}}_{\text{BC}},x,t}^{\ast }={\displaystyle \frac{{\mathbf{HS}}_{\text{MOD},x,t}}{{\stackrel{\mathrm{‾}}{\mathbf{HS}}}_{{\text{MOD}}_{\mathrm{1962}-\mathrm{1971}},x}}}\times {\stackrel{\mathrm{‾}}{\mathbf{HS}}}_{{\text{OBS}}_{\mathrm{1962}-\mathrm{1971}},x}^{\ast },\end{array}$$

where ${\mathbf{HS}}_{{\text{MOD}}_{\text{BC}},x,t}^{\ast }$ denotes our estimated anomalies at station x and time t, which have been bias-corrected to remove any mean biases present over the common period of 1962–1971. Equation (8) is applied to both models to provide the estimated anomalies which have been mean-bias-corrected for all of the stations (i.e., all Austrian stations for SNOWGRID-CL and all Austrian and Swiss stations for SnowSens).

4 Results

4.1 Model performance

4.1.1 Comparing the SnowSens and SNOWGRID-CL models

The skill of the SnowSens model is evaluated with respect to the SNOWGRID-CL model. The SNOWGRID-CL model is run over the Austrian domain, and, as a result, the performance of the SnowSens model is evaluated in this section only using the 107 stations within Austria.

Prior to any bias correction, the SNOWGRID-CL estimates perform about as well as or worse than the climatology. The performance of the SNOWGRID-CL absolute model values over the period of 1972–2021 can be evaluated using different periods to calculate the climatological mean or what is treated as normal. When this period is 1902–1971, SNOWGRID-CL has an SS_RMSE value of 0.03, and when the 1972–2021 period is treated as the climatological normal, SNOWGRID-CL has an SS_RMSE value of −0.09. In contrast, the absolute estimates of the SnowSens model have greater skill, with SS_RMSE values of 0.25 and 0.16, respectively, for the two different climatological periods (the skills with respect to climatologies computed over the 1902–1971 period can be found in Table 2).

Table 2 A comparison of different skill scores for the two models over the Austrian domain. In the top two rows, the absolute year-to-year skill (i.e., using the raw model estimates) is shown for each elevation band and for all stations. The middle three rows give the year-to-year model skill using the estimated anomalies, which have been bias-corrected, for each of the two models and an ensemble average of the two. The last three rows provide skill scores of how well the models estimate relative changes between the period of 1997–2021 and the period of 1972–1996.

After applying mean bias correction to both models, the performance of the SNOWGRID-CL model is much improved. Now, using the estimated anomalies which have been bias-corrected, the skills (i.e., SS_RMSE) of the SNOWGRID-CL and SnowSens models over the 1972–2021 evaluation period are 0.34 and 0.26, respectively. Figure 8a–d show the modeled seasonal snow depth anomalies against the observed seasonal snow depth anomalies for the two models for each of the four elevation bands. While both contain statistically significant skill (p<0.01, where we used bootstrapping to assess the statistical significance), the bias-corrected SNOWGRID-CL model is found to be more skillful than the SnowSens model when it comes to modeling the year-to-year variability of the seasonally averaged snow depths.

Figure 8 SnowSens- and SNOWGRID-CL-estimated seasonal (November–March) snow depth anomalies (x axis) plotted against observations (y axis) for stations which fall in each of the four elevation bands, as shown in (a)–(d), respectively. The estimated changes in seasonal snow depth for the two models are plotted against observed changes in (e). Changes in (e) are the differences, at each station, in the anomalous seasonal snow depths between the more recent period of 1997–2021 and the prior period of 1972–1996. The black lines in (a)–(e) show the one-to-one lines.

Next, we want to know how well the observed trends over the evaluation period have been modeled by both SNOWGRID-CL and SnowSens. The relative changes are computed for all of the Austrian stations between the period of 1997–2021 and the period of 1972–1996. This is done for the two models and for the observations. These values are plotted in Fig. 8e. When it comes to correctly modeling the trend, the SnowSens model now outperforms SNOWGRID-CL. SNOWGRID-CL generally overestimates the changes over the last 50 years. The skill score, SS_RMSE, for the modeled versus observed changes over the evaluation period is 0.19 for the SnowSens model and 0.10 for SNOWGRID-CL (see Table 2 for a number of comparative skill scores).

We also computed an ensemble average of the two models using their estimated anomalies, and this was found to be more skillful than either model alone. This is true for both the skill scores of year-to-year variability and relative changes observed in the last 50 years. Given this result, using an ensemble such as that between these two models has the potential to further improve the projections of future seasonal snow depth.

4.1.2 Model performance for the entire domain

In the last section, we compared how well the estimates of the SnowSens and SNOWGRID-CL models performed at the locations of the 107 stations across Austria. We investigated the ability of the models to capture both the year-to-year variability and the historical trends of the observed records. Here, we present the results of the SnowSens model performance for all of the stations. Since we do not have SNOWGRID-CL estimates for the Swiss domain, we are now only showing the results for the SnowSens model.

In the previous section, we found that the SNOWGRID-CL model better captures the observed year-to-year variability at the station scale. However, we want to be clear that the SnowSens model still exhibits substantial skill for the year-to-year seasonal estimates at the station scale. Figure 9 shows all of the estimated seasonal snow depth anomalies against the observed anomalies for the 1972–2021 validation period. The Pearson correlation coefficient between estimated and observed anomalies for these 10 985 cases is 0.59, with an SS_RMSE of 0.25 (p<0.01). We used bootstrapping to assess the level of statistical significance, where different shuffled seasons from the calibration period are used as estimates for the validation period. That way, the spatial autocorrelative structure of the seasonal snow depths is preserved in our bootstrapping procedure. We evaluated the skill of 10 000 bootstrapped simulations and found the largest SS_RMSE value from these randomly generated simulations to be 0.022 (compared to the value of 0.25 from the SnowSens model). While we present our p value as less than 0.01, it is actually much smaller. Please also refer to Table 3 for some additional model skill scores using all of the seasons and all of the stations in the study domain.

Figure 9 All seasonal SnowSens model estimates versus observations over the validation period using all of the stations in the study domain. Again, the black line plots the one-to-one line.

Table 3 Skill scores of the SnowSens model for different elevation bands over the entire domain. The top and middle rows are the skill scores using the absolute and estimated anomalies, respectively. The bottom row is the skill in estimating the relative changes between the 1997–2021 and 1972–1996 periods. The three skill scores in the right-most column (which includes all stations over the entire domain) are statistically significant with p<0.01.

We would also like to know how well the model performs in its ability to extrapolate to new climatological terrain. Put another way, how well does the model perform in cases which it had not really seen before in the calibration period? In Fig. 10, we test the effectiveness of the SnowSens model in a climatological region that had rarely seen observations in the calibration period. Figure 10a shows the 95 cases where the average seasonal temperature in the validation period was greater than +1.0 °C and less than 50 % of normal precipitation. For these cases, the error of the SnowSens estimates is less than half that of the climatological estimates (indicated by SS_RMSE>0.50). The SS_RMSE for the points in Fig. 10a is 0.62, which is statistically significant with p<0.01. Note that nearly all of the modeled and observed values fall below 100 % of normal. Averaging across this set of cases (the larger green squares), we find that both the modeled average and the observed average are 33 % of normal. In Fig. 10b, we increase the sample size by using a threshold of less than 75 % of normal precipitation (instead of 50 %). This gives us 988 cases. The SS_RMSE in this case is 0.55, which is also found to be statistically significant with p<0.01. The modeled and observed averages over these cases are 42 % and 41 %, respectively. So, while we are extrapolating to “unknown” climatological terrain, we find that the model is quite capable of performing skillfully in that new terrain, especially when aggregating over a number of cases. We do advise, however, that one proceeds with caution when interpreting any of the individual model estimates.

Figure 10 Figure 10 shows the effectiveness of the SnowSens model in terms if its ability to estimate in new climatological terrain. Panel (a) plots the pairings of modeled versus observed anomalies for the 95 cases where the average seasonal temperature in the validation period was greater than +1.0 °C and less than 50 % of normal precipitation. Panel (b) plots the pairings for the 988 cases when using a threshold of less than 75 % of normal precipitation (instead of 50 %). The larger green squares show the modeled versus observed averages across the respective cases from (a) and (b). Again, the one-to-one lines are plotted in black.

A spatial plot showing the geographical distribution of the station-by-station skill scores is shown in Fig. 11a. With the help of Fig. 11b, it can be observed that the skill of the model estimates generally decreases with elevation (it is not shown, but the SNOWGRID-CL model also sees decreasing skill with increasing elevation). This makes sense, given that stations at lower elevations are more sensitive to temperature changes, and the ranges of these temperature changes at lower elevations are particularly large in the context of climate change. Bootstrapping is used again to test for the statistical significance of the model skill at each individual station (Fig. 11a and b). We find that there are only a total of 11 stations which either perform worse than the climatology or perform worse than randomly estimated time series (with p<0.01). Therefore, we find that the SnowSens estimates of snow depths from more than 95 % of the stations in our study exhibit statistically significant positive skill. Figure 11c plots the modeled and observed relative changes in seasonal snow depth between the periods 1972–2021 and 1902–1971, while Fig. 11d shows this between the periods 1997–2021 and 1972–1996. The yellow diamonds are relative changes averaged across the four elevation bands. The SS_RMSE of the band-averaged estimate changes for these two periods is 0.80 and 0.73, respectively. By implementing spatial averaging of the estimates across elevation bands, the skill improves dramatically. We should also note that we also applied the same methodology to snow water equivalent (SWE) values that we constructed via the approach outlined in Winkler et al. (2021). While not shown here, the skill of the SWE estimates follows the skill in estimating snow depth very closely.

Figure 11 (a) Skill over the 1972–2021 evaluation period at each of the stations. (b) Plotting how the skill varies with respect to elevation. (c) Modeled versus observed relative changes in snow depth between the periods 1972–2021 and 1902–1971. (d) Same as (c) but using the periods 1996–2020 and 1972–1996. The larger yellow diamonds in (c) and (d) show the changes averaged across the four elevation bands.

Figure 12 Average (November–March) estimated snow depth anomalies and observations are plotted as time series for each of the four elevation bands. The thinner and thicker lines show the average anomalies during the calibration and validation periods, respectively. The text in the upper right of each panel lists a select number of metrics (i.e., observed trend per decade, modeled trend per decade, SS_RMSE, and Pearson correlation coefficient) corresponding to the 1972–2021 validation period. We should also note what the correlation coefficients are between the modeled and observed time series when the trends of both time series have been removed over the validation period. In those cases, the correlations are 0.89, 0.83, 0.79, and 0.74, respectively, which are very similar to those that are reported in the figure with the trends present. The number of stations used to compute the band averages is 52, 75, 91, and 73, respectively (see Table 1).

In Fig. 12, estimated anomalies of seasonal snow depth are averaged across all of the stations which fall within each of the four elevation bands. The resulting band-averaged time series can be seen alongside the band-averaged observed time series in the panels of Fig. 12. Again, we find that SnowSens does not fully capture the observed year-to-year variability and is not able to reproduce the high and low extreme values. However, we do find that the model performance is further improved using these band averages.

Figure 13 This figure plots cross-sections of Fig. 7 at 80 %, 100 %, and 120 % of normal precipitation. Within a climate change context, this figure serves to provide greater ease in quantifying the changes in snow depth given a range of projected changes in temperature and precipitation. The vertical pink lines show how much warming has already taken place in each elevation band over the period of 1972–2021 with respect to the period of 1902–1971. The average (November–March) estimated and observed snow depth anomalies over the 1972–2021 validation period are plotted, respectively, as the open square and the “×”.

Figure 13 provides a useful and simple-to-interpret plot of expected future changes in snow depth as a function of temperature and precipitation anomalies. A user can take a range of expected future projections of temperature and precipitation and evaluate how these would translate into expected changes in snow depth. These expected changes in snow depth given different temperature and precipitation scenarios can be compared to prior studies, such as those of Schmucki et al. (2014) and Marty et al. (2017a). In Fig. 13, we have simplified the information provided in Fig. 7i–l. Figure 13 shows cross-sections of Fig. 7i–l for three different precipitation anomalies, and those are 100 % of normal, along with 80 % and 120 % of normal (i.e., 20 % below and above average). The vertical pink lines show how much warming has already taken place in each elevation band over the period of 1972–2021 with respect to the period of 1902–1971. The average estimated and observed snow depth anomalies are plotted, respectively, as the open square and the “×”. One can use this plot to gain a more detailed understanding of how something like an additional 2 °C would translate to snow depth anomalies at different elevations, given the assumption that precipitation stays about the same (100 % of the 1902–1971 normal). For elevations below 500 m, an additional 2 °C (which is 3.2 °C above the 1902–1971 normal) could lead to seasonally averaged accumulated snow depths being a thing of the past. Put another way, there would be nearly no snow depth accumulation projected at those temperature anomalies. Given 2 °C more warming in the other three elevation bands, we could expect snow depth anomalies of approximately 25 %, 50 %, and 80 % of normal, respectively.

5 Conclusions

Climate change has already had an observable impact on the average seasonal snow depths across the European Alps. Over the historical period of 1902–2021, stations across Austria and Switzerland showed a decrease in seasonally averaged November–March snow depth, ranging between −8.1 % per decade at elevations below 500 m and −3.4 % per decade for elevations above 1500 m. Over the more recent historical period of 1952–2021, these changes were even greater, with decreases ranging between −15.7 % and −3.6 % per decade for stations below 500 and above 1500 m, respectively (see Appendix A). Changes in seasonally averaged snow depth can primarily be attributed to changes in meteorological forcing variables such as mean temperature and precipitation. In some cases, blowing wind and sublimation of snow can greatly affect the snowpack, and, in these instances, temperature and precipitation alone may not suffice as predictors. In this paper, however, we focus exclusively on using anomalies of temperature and precipitation across different elevation bands as predictors of seasonally averaged snow depth.

Using historical observations of seasonally averaged temperature, precipitation, and snow depth at four different elevation bands over the period of 1902–1971, we constructed a multivariate empirical–statistical model, which is named SnowSens. Model validation, which was performed over the period of 1972–2021, shows that both the SnowSens and SNOWGRID-CL models can skillfully estimate year-to-year seasonally averaged snow depths across the Austrian domain. While SNOWGRID-CL is found to perform better in estimating the year-to-year variability of snow depth, the SnowSens model estimates historical trends better. The SnowSens model is not to be seen as a replacement for operational models such as the SNOWGRID-CL. Rather, this paper highlights how effectively historical sensitivities can be used in a multivariate framework to produce quite accurate estimates of long-term changes in snow depth. The model performs particularly well when the estimates are aggregated over a number of cases, such as across different elevation bands (as seen in Fig. 12). Furthermore, SnowSens relies on a comparatively simplified modeling framework, which lends itself well to easily translating projected changes in temperature and precipitation to changes in snow depth. Our results show that the historical sensitivities have been robust and persistent. If these sensitivities continue to remain persistent into the future and if future projections of temperature and precipitation are also skillful then this modeling approach can be expected to yield skillful forecasts for the next 50 years.

The impacts of a changing climate will vary from region to region. We developed multivariate sensitivities that are regionally specific to the Austrian and Swiss Alps. While outside of the scope of this study, the same approach can be applied to other mountain regions. How might the sensitivities of the Rocky Mountains or the Cascades of the United States compare to what we observe in the Alps? Another way our research can be extended relates to quantifying the changes in snow depth versus streamflow. Are specific reductions in snow depth at certain elevations noticeably affecting aggregated streamflow measurements? Or is it rather the timing of discharge that is impacted? And lastly, it would be valuable to investigate how capable a variety of different general circulation models (GCMs) and regional climate models (RCMs) are in capturing the observed sensitivities that we have produced.

Snow depth is a valuable resource that affects many communities adjacent to and downstream of mountain regions. Changes in snow depth can have broad impacts that range from water resources to snow tourism and avalanche preparedness. Climate change is expected to bring about further increases in temperature across the Alps, while it is less clear what the impact will be on precipitation. With improved tools, we can better quantify the impact that these meteorological changes will have on snow depth, thus allowing communities to better plan and prepare for the changes to come.

Appendix A

In Fig. A1, we plot the historical trends over the more recent period of 1952–2021. There are greater relative decreases in snow depth over a more recent historical record of 1952–2021 (in contrast to Fig. 3). In this more recent period, we observe decreases in snow depths ranging between −15.7 % per decade (Fig. A1a) for stations below 500 m and −3.6 % per decade for stations at elevations above 1500 m (Fig. A1d). Anthropogenic climate change is often more clearly recognizable in the recent past, and, hence, the trends derived for this period can help to improve our understanding of the expected future changes in snow depth.

Figure A2 plots the raw changes in snow depth over the last 30-year period of 1992–2021 compared to the period of 1962–1991 as a function of elevation. Greater absolute changes in snow depth are observed as elevation increases. However, this information needs to be placed in the context of differences in climatology. As is shown in Fig. 2c, the average seasonal snow depth scales with elevation. The relative changes in snow depth over the last 30 years are plotted against elevation in Fig. A2b. Using anomalies, we obtain a stronger relationship between the more recent changes in snow depth and elevation (compare the correlations between Fig. A2a and b). The explained variance between elevation and relative changes is approximately 34 %, while it is only about 10 % when using the absolute changes. In addition to giving us a better statistical relationship, we also get a clearer picture of where we can observe the greatest relative changes in snow depth. While stations at lower elevations saw smaller absolute changes in their snow depth over the last 30 years, these same stations saw a greater relative decrease over the same period of time. For example, stations below 500 m had an average absolute change of −1.58 cm, but this reasonably small absolute amount reflected a large relative change, which was 71 % of normal. In contrast, stations above 1500 m had an average absolute change of −8.38 cm, while the average relative change was 93 % of normal. Next, we observe where the changes in snow depth have been the greatest, or the least so, as a function of geographic location. To do this, we first isolate the influence of latitude and longitude (i.e., x and y space) by removing the dependence of these snow depth changes on elevation (i.e., z space). The dashed line in Fig. A2b shows the fitted line through least-squares linear regression of the data. The data are detrended with respect to the regression line while preserving the population mean (Fig. A2c). The downward-facing and upward-facing triangles in Fig. A2d show anomalous snow depth conditions over the period of 1992–2021, where the influence of elevation has been removed. Figure A2d, along with Fig. A2e and f, shows that, as one traverses the Austrian and Swiss Alps from east to west and from north to south, the stations experience slightly greater relative decreases in snow depth.

Figure A1 Same as Fig. 3 but using the 1952–2021 period of record.

Figure A2 Changes in snow depth (1992–2021 vs. 1962–1991). Panel (a) shows the absolute snow depth changes as a function of elevation. Panel (b) shows the relative or anomalous snow depth changes as a function of elevation. The dashed line is the fitted least-squares regression of the data. Panel (c) plots the snow depth anomalous changes where the elevation dependence has been removed by detrending the data with respect to the dashed line from (b). The solid lines in (b) and (c) show 100 % of normal with the elevation dependence removed, while the colors and direction of the triangles reflect the stations which experienced positive or negative changes after removing the elevation dependence. Panel (d) plots the detrended anomalous changes (from Fig. 8c) across the study region. The size of the triangles reflects the size of the anomalies. Panels (e) and (f) show the detrended anomalous changes against longitude and latitude, respectively.