Local- and Regional-Scale Forcing of Glacier Mass Balance Changes in the Swiss Alps

Abstract

Glacier mass variations are climate indicators. Therefore, it is essential to examine both winter and summer mass balance variability over a long period of time to address climate-related ice mass fluctuations. In this study, we analyze glacier mass balance components and hypsometric characteristics with respect to their interactions with local meteorological variables and remote large-scale atmospheric and oceanic patterns. The results show that all selected glaciers have lost their equilibrium condition in recent decades, with persistent negative annual mass balance trends and decreasing accumulation area ratios (AARs), accompanied by increasing air temperatures of ≥ +0.45 °C decade⁻¹. The controlling factor of annual mass balance is mainly attributed to summer mass losses, which are correlated with (warming) June to September air temperatures. In addition, the interannual variability of summer and winter mass balances is primarily associated to the Atlantic Multidecadal Oscillation (AMO), Greenland Blocking Index (GBI), and East Atlantic (EA) teleconnections. Although climate parameters are playing a significant role in determining the glacier mass balance in the region, the observed correlations and mass balance trends are in agreement with the hypsometric distribution and morphology of the glaciers. The analysis of decadal frontal retreat using Landsat images from 1984 to 2014 also supports the findings of this research, highlighting the impact of lake formation at terminus areas on rapid glacier retreat and mass loss in the Swiss Alps.

1. Introduction

Glacier mass balance, the difference between accumulation and ablation, is an important prerequisite for monitoring direct and immediate impacts of climate forcing on the cryosphere [1]. Glaciers react to climatic conditions by losing or gaining mass [2,3,4]. This reaction is mainly controlled by two processes: snowfall precipitation within the accumulation zone in winter and snow–ice melt, runoff, and sublimation during the summer. In connection with ongoing global atmospheric warming and rising air temperatures, the ablation season has expanded and resulted in the acceleration of glacier volume loss [2,5]. Therefore, changes in seasonal mass balance can be considered a clear indicator of climatic variations in high-mountain regions [6]. Since glacier mass balance is linked directly to local meteorological conditions and large-scale atmospheric circulation, as well as glacier dynamics and hydrology, a thorough understanding of glacier mass balance and its interactions with climatic variables is critical to local communities and regions because of hydropower supply [7,8,9,10,11,12], irrigation systems [13,14], hazard monitoring [15,16,17,18], and opportunities like geotourism [19,20,21,22].

There is a long-standing tradition of measurements of European glacier’s surface mass balance conditions. The first mass balance concepts involved a series of accumulation and ablation measurements on Nordic glaciers in the 1920s and 1930s [23]. However, glacier monitoring started in 1894 with a more simplified focus on glacier length changes [24]. After the late 1940s, interests were shifted to collecting data on glacier-wide mass balance [25]. The first standardized data publication on change in glacier length, area, volume, and mass started with the ‘Fluctuations of Glaciers’ series in the late 1960s [2,26]. Since the 1980s, data on glacial changes has been documented in a relational database, Worldwide Glacier Monitoring, for access by scientists, policy makers, and the public [2,27]. The World Glacier Monitoring Service (WGMS) plays a critical role by compiling available data through a scientific and international network based on contributions from thousands of observers from all over the world. These archives are based on common measurement protocols and are currently available through the WGMS and Global Terrestrial Network for Glaciers (GTNG) as part of the Global Climate Observing System. Mass balance data are obtained either from direct glaciological methods (field measurements) or from satellite/aerial monitoring change on glacier surface topography, which converts the volume to mass change using a density assumption [28,29].

Compared to other mountain ranges around the world, the European Alps have the longest record of mass balance data [23]. According to available reports and previous studies, these glaciers have demonstrated great mass loss since the Little Ice Age. Zemp et al. [30] modeled glacier mass balance fluctuations for the European Alps dating back to 1850 using the inventory from the WGMS and the Swiss Glacier Inventory of 2000. They reported an overall estimated loss of 35% in glaciated areas from 1850 to the 1970s, followed by a 50% acceleration in mass loss from 1970s to 2000, visible mainly as terminus retreat. It has been suggested that such fluctuations in glaciers are the indirect and lagged response to past climatic conditions [2,31,32]. Huss et al. [31] re-analyzed seasonal glacier-wide mass balance in the southeastern Swiss Alps from 1900 to 2008 based on field measurements. Using a distributed accumulation and temperature-index model, the authors concluded that the rate of mass loss in this time frame varied significantly between adjacent glaciers because of the dynamic adjustment to the post-Little Ice Age climate change, namely warm air temperatures. Their study also projected a 63% decrease in glacier extent through 2050 and increased annual runoff over the next three decades for the deglaciating catchments. Fischer et al. [1] extended a re-evaluated inventory of geodetic mass balance, surface elevation change, and volume change for the entire glaciers in the Swiss Alps from 1980 to 2010. The estimated average mass balance was −0.62 ± 0.07 m w.e. yr⁻¹, with an overall volume loss of −22.51 ± 1.76 km³. The findings of this research also revealed the mass change for the glaciers in Valais Alps (Rhone) to be lesser (−0.59 m w.e. yr⁻¹), and that for the southern glaciers (e.g., Maggia, Ticino, Maira, Poschiavino) to be greater, than that for the entire Swiss Alps. Following that study, Huss et al. [33] re-evaluated 19 series of the Swiss glacier-wide seasonal mass balance from 1920 based on point observations to establish a homogenous mass balance series. This dataset, which includes winter (1 October–30 April) and summer (1 May–30 September) observations back to the 1980s, offers a strong basis for analyzing the response of glacier mass balance components to climate variability and change.

Previously, numerous studies have addressed the sensitivity of mass balance to various climatic factors including changes in surface air temperature and precipitation [3,34,35,36,37,38,39]. Huss and Fischer [40] developed a model, the Glacier Evolution Runoff, to evaluate the response of very small glaciers (<0.5 km²) in Switzerland to current and future atmospheric warming. According to their study, Swiss glaciers are projected to disappear in the next 25 years. Glacier mass balance, in general, has also been linked to oceanic and atmospheric variability through various climate indices (e.g., the North Atlantic Oscillation (NAO), the East Atlantic (EA) pattern, and the Atlantic Multidecadal Oscillation (AMO)) [3,41,42,43,44]. The impact of such climate forcing on Swiss glaciers have been addressed in a few studies [6,32,45]. Huss et al. [32] associated North Atlantic climate variability with the 100-year surface mass balance record for thirty Swiss glaciers. The results found the AMO index inversely correlated with glacier mass balance, suggesting that positive AMO values (warm sea surface temperatures) impacted significant mass loss from 1908 to 2008.

Even with a rich literature dedicated to monitoring long-term glacier trends, there is still a further need to perform time series analysis to quantify Swiss glacier mass balance changes and drivers as they relate to local and remote climatic conditions, specifically for glaciers included in newly developed datasets [33]. Moreover, there is a special need to study the development and evolution of catchments through ice melt and associated features such as new lakes, since they are closely tied to glacier ice loss [46]. As such, our study employed glaciological, remote sensing, and climate datasets and assessed glacier mass balance components and hypsometric characteristics with respect to interactions with local meteorological variables and remote large-scale atmospheric and oceanic patterns to more comprehensively understand physical processes driving glacier mass balance variability.

Precisely, we analyzed the behaviors of the southwestern Swiss glaciers in the context of post-Little Ice Age climate change. We also conducted a long-term (1970–2017) time series analysis of annual glacier mass balance (1 October–30 September), as well as winter (1 October–30 April) and summer (1 May–30 September) mass balances, derived from glaciological reports (1881–2014) contained in the WGMS [26]. Additionally, we analyzed and compared the time series data to identify trends, fluctuations, and similarities in glacier annual and seasonal mass balances. We further provided hypsometric characteristics of these glaciers using swissALTI^3D Digital Elevation Model (DEM) and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global DEM (GDEM) Version 3 (V3) and examined the possible impacts of geometry on the recorded responses. Using Landsat satellite images covering a span of 30 years, we investigated the changes in the rate of terminus retreat for these glaciers. We further explored impacts of local climatic variations, including the air temperature fluctuations during ablation season (May–September) and precipitation during accumulation season (October–April). We then evaluated multidecadal variations of glacier mass balance and local climate to forcing from regional atmosphere–ocean patterns, including the AMO, NAO, and EA, as well as the Scandinavian pattern (SCA) and Greenland Blocking (described by the Greenland Blocking Index (GBI)), since the atmospheric patterns described by these indices play a key role in modulating European weather and climate [47].

2. Study Area

This study is focused on seven glaciers (Grosser Aletschgletscher, Allalin, Hohlaub, Schwarzberg, Gries, Gietro, and Rhone) of varying size, geometry, hypsometry, aspect, and slope, in the Upper Rhone catchment (Figure 1 and Table 1). Grosser Aletschgletscher, which is located within the Bernese Alps, is the longest (22.6 km) and largest (83.0 km²) ice body in the European Alps, with an extensive archive of climate data and glaciological records dating back to 1841 [22,48]. The accumulation area of this glacier faces south and contains two summits, Jungfrau and Mönch. The main body is divided by three branches (Grosser Aletschfirn, Jungfraufirn, and Ewig Schneefeld) at Konkordiaplatz with reported ice thickness of 890 m at its summit. According to the reports, it has lost significant volume since the Little Ice Age with a fast acceleration since the late 1980s [26].

Rhone is a south-facing valley glacier with a climatic condition similar to Aletsch. Gries is different in this regard than the other glaciers as it is placed south of the main Alpine crest at the border of Italy within a small valley exposed in a northeastern direction with relatively high precipitation. According to Huss et al. [6], Gries and Aletsch have lost 39% and 16% of their ice volume in the 20th century, respectively. Mass balance monitoring programs were established on Gries in the 1960s with hydropower projects, whereas the data collections on Aletsch and Rhone were taken as a part of glaciological research projects.

The Mattmark region includes three glaciers (Allalin, Hohlaub, Schwarzberg), with the longest mass balance record dating back to 1955, established during the construction of Mattmark dam in the Saas Valley. The Laboratory of Hydraulics, Hydrology, and Glaciology at ETH Zürich (VAW/ETHZ) has been observing the glaciers in corresponding catchments, including Gietro and two other glaciers, Allalin and Schwarzberg, in the Mattmark area (Saastal). Allalin is a large valley glacier flowing in a northeast direction and covers 9.6 km² of the Saas Valley. This glacier has the highest equilibrium line altitude (ELA) of 3575 m a.s.l. owing to the little solid precipitation in the valley [26]. According to the Swiss glaciers report in 2014/2015, the least amount of mass loss, relative to the inventory, was recorded for this glacier, with average thickness changes of approximately −0.6 m w.e. [26].

Two additional glaciers, Silvretta and Clariden, located outside of the catchment in the Central and Glarus Alps, respectively, have the longest records in the Swiss glacier database since the 1900s, and are included to understand the regional variations between the glaciers and provide information on the temporal and spatial patterns of mass balance fluctuations within the Swiss Alps. The glaciers vary in size from 2 to 83 km² as documented by VAW/ETHZ glaciological reports (1881–2014), which includes stake-derived mass balance measurements as well as geometries, aspects, and local climatic conditions of individual glaciers.

3. Data

A wealth of data is available for the Swiss glaciers, especially in the Upper Rhone catchment. The combination of mass balance measurements (Table 1) and climate data allows us to draw an integrative picture of the local and regional forcing of mass balance variations in this region.

3.1. Mass Balance Data

For the Swiss Alps, two direct mass balance series are available that span five continuous decades (Table 1, Figure 2) [26,33].

These mass balance time series are calculated by glaciological field methods including dense snow probings and pits dug in April–May that measured the winter snow accumulation, and ablation stake readings in September [40]. The longest mass balance time series (since 1900) for 50 glaciers in the Swiss Alps was reconstructed by Huss et al. [33] based on the combination of field measurements and ice volume changes by comparison of DEMs. The extrapolation of individual measurements to the entire glacier surface was performed using a mass balance model which considers all important processes controlling glacier mass balance distribution. There are seven glaciers in the Upper Rhone catchment with a long sequence of mass balance measurements recorded since 1884. In this dataset, winter mass balance (Bw), corresponds to the accumulated snow/ice mass measured from October to the end of April. Summer mass balance (Bs) is the result of difference between annual mass balance (Ba) determined at the end of the ablation period (end of September) and Bw. The length of annual mass balance and time series of Ba, Bw, Bs, and AAR for each glacier are presented in Table 1 and Figure 2a,b.

3.2. Climate Data

Monthly mean air temperature and precipitation data are used to determine the relationship between the mass balance components and local climate. Long-term homogenized series of average daily air temperatures (two meters above surface) and total precipitation from the Swiss National Basic Climatological Network (NBCN) stations considering the underlying topography, with coverage since 1961 over the entire Swiss Alps, were provided by MeteoSwiss in a 2 × 2 km grid structure [49,50]. The daily mean meteorological time series is extracted from the gridded area covering the glaciers.

In order to assess the relationship between the components of mass balance and the larger-scale atmosphere and oceanic conditions, seasonal (October–April and May–September) and annual times series of the NAO, SCA, and EA were obtained from the Climatic Prediction Center (CPC), and those of the GBI [44] and AMO were downloaded from the National Oceanic and Atmospheric Administration (NOAA). The GBI represents the normalized 500 hPa mean GPH anomalies across Greenland, 60–80° N, 20–80° W [51,52,53,54]. The NAO has been defined in several ways, and we use two versions here. The Hurrell NAO index used here is based on an empirical orthogonal function (EOF) analysis of North Atlantic sea-level pressure, which tends to capture spatial variations in the action centers generally found over Iceland and the Azores [55]. To evaluate the upper-air circulation features over a similar geographic area, the CPC NAO index is used supplementally and is identified as the first EOF derived from a rotated EOF analysis of 500 hPa anomalies from 20–90° N. The corresponding second and third EOFs, the EA and SCA, are also analyzed to understand a broader array of North Atlantic synoptic settings potentially linked to the Swiss glaciers variability and changes. The EA pattern tends to have an action center near 80° N and 20–35° W, while the SCA coincides with a high-pressure center over the Scandinavian countries during its positive phase [56,57]. The AMO index used here is based on detrended North Atlantic sea surface temperature anomalies from the Kaplan SST dataset [58,59]. This teleconnection is linked to continental climate variations in both North America and Europe [58]. Positive AMO values are associated with warm near-surface air temperature anomalies in Europe, and prevailing low-pressure anomalies over the Atlantic and Europe during the winter (DJF) [60].

3.3. Remotely Sensed Data

For hypsometric analysis, we use swissALTI^3D DEM obtained from the Swiss Federal Office of Topography (Swisstopo) for the seven glaciers in the catchment and the newly released ASTER GDEM Version 3 (V3) data for the two additional glaciers located in the Central and Glarus Alps. This swissALTI^3D DEM (2019 edition) represents the surface glaciers at a spatial resolution of 10 m with ±1 m to 3 m (1σ) vertical accuracy for the stereo correlation and ±0.1 m to 1 m for manual stereo measurements. Areas below 2000 m a.s.l. were created using airborne laser scanning data with an accuracy of ±0.5 m (1σ) [61]. The ASTER GDEM V3 was released in 2019 and displays similar level of accuracy as ASTER GDEM V2. However, V3 has lesser void area due to the increase in ASTER stereo image data since the previous release and improved processing [62].

To observe the glacier termini changes over a period of three decades (1984–2014), we use 12 Landsat 5 Thematic Mapper (TM) and Landsat 8 Operational Land Imager (OLI) images downloaded from EarthExplorer data portal [63]. The images are selected if they are completely cloud-free over the termini of the glaciers and show least seasonal snow cover for the selected timeframe. The collected images are then spread across a 10-year time interval, except for Gietro (Table 2), due to the unavailability of cloud-free coverage of its terminus within the same years as for other glaciers.

4. Methods

To achieve the research objectives, we adopted a coherent approach combining statistical methods for glaciological and climate time series with geospatial methods for remote sensing data (Landsat and DEM) analyses to make meaningful inferences. Below, we have compiled the methods within various sub-sections to detail on the key steps for any such analyses in future works.

4.1. Statistical Analysis

4.1.1. Exploratory Data Analysis

Using the WGMS database, time series of Ba, Bw, Bs, and AAR (the ratio of accumulation area at the end of melt season to total glacier area) are calculated and assessed for inhomogeneities (i.e., missing values (Figure 2), and any change in unit measurements), testing assumptions to select the analysis strategy (e.g., presence of autocorrelation within the time series), trends and outliers, fluctuations, and the general pattern of glacier behaviors (i.e., mean and standard deviation of Ba, Bs, Bw, and AAR of glaciers) within the common period of observation (1970–2017). Exploratory analyses are initially conducted, including time series plots, boxplots, autocorrelation functions, correlation matrices, and smoothing techniques (i.e., moving averages). To identify similarities between glaciers and address the seasonal aspect of mass balance, the Ba series of all the glaciers are intercompared and linked with Bw and Bs characteristics using Spearman correlation analysis. Prior to conducting correlation analyses, data series are detrended using a triangular moving average approach to detect and highlight intradecadal changes [44,64,65]. A seven year level of triangular average is applied as it yields optimal correlations between the annual and seasonal mass balance of all glaciers versus seasonal and annual climate indices.

4.1.2. Structural Change Model

Various approaches have been developed in the literature for testing structural changes and identifying breakpoints in time series. Most tests examine coefficient changes of a linear regression model where the assumption is a single change in time and the type of change is known [66]. Another method referred to as multiple structural changes has been widely used [66,67,68,69,70] to test the presence of multiple shifts over time. We consider a linear regression model as:

(1)$y_i=x_i^T{\beta }_i+{\mathrm{\mu }}_i\left(i=1,\dots , n\right)$

where at time $i$, $y_i$ is the observation of the dependent variable, $x_i$ is a k × 1 vector of regressors, and ${\beta }_i$ is the k × 1 vector of regression coefficients [71]. Here, the null hypothesis is that the regression coefficients are constant, while the alternative hypothesis is that at least one coefficient changes over time:

(2)$H_0 :{\beta }_{i }={\beta }_0\left(i=1,\dots , n\right)$

This assumes that there are m breakpoints, where the coefficients deviate from one stable regression relationship to a different one [66]. Thus, there are m + 1 segments in which the regression coefficients are constant, and model (1) can be rewritten as:

(3)$y_i=x_i^T{\beta }_i+{\mathrm{\mu }}_i\left(i=i_{j-1}+1,\dots , i_{j },j=1,\dots , m+1\right)$

where j is the segment index, and I_m,n = {i₁, …, i_m} refers to the set of breakpoints or the number of m partitions [66].

Overall, two approaches for testing structural change are used: (i) F statistics [72], and (ii) generalized fluctuation tests [73]; the former approach tests against a single-shift alternative of unknown timing, whereas the latter technique, also known as empirical fluctuation process (efp), is concerned with constancy in a graphic method reflecting the fluctuations in residuals and coefficient estimates [74]. The null hypothesis of the generalized fluctuation test is based on the central limit theorem, and if the efp boundaries are crossed by fixed probability α, the fluctuation is large, and the null hypothesis of no structural break is rejected. The efp function from the fluctuation test framework was programmed in R and the cumulated sum of standard OLS residuals (CUSUM) function was specified as the type of fluctuation process to be fitted in the argument. CUSUM detects change points in time series that often results from structural changes based on the limiting process of standard Brownian motion [75]. Since the number of breakpoints m is not known in advance, it is necessary to compute the optimal breakpoints for m = 0, 1, …n breaks and choose the model that minimizes some information criterion such as the Bayesian information criterion (BIC) [76]. Following this exercise, a Chow test [77] is also performed to confirm the statistical significance of resulting change points (10-year sliding window) at the 95% confidence level when p ≤ 0.05.

4.1.3. Trend Analysis

Mann–Kendall and Sen’s Slope Estimator

To evaluate whether mass balance values increase or decrease over time, the Mann–Kendall test is applied. This test, which is a nonparametric form of monotonic trend regression analysis, is robust to missing values and outliers in time series. The test assumes that the data are independent, and the distribution of data remains consistent [78]. Therefore, autocorrelation (ACF) and partial autocorrelation (PACF) tests are conducted, and serial correlation is not found in the series which further motivated application of this approach.

The Mann–Kendall test calculates the difference between the later value and all previous values, (x_j − x_i), where j > i, and assigns the value of 1, 0, or −1 to positive, no differences, and negative differences, respectively [79,80]. The test statistic, S, is then calculated as:

(4)$S={\displaystyle \sum }_{i=1}^{n-1}{\displaystyle \sum }_{J=i+1}^{n}sgn\left(x_j-x_i\right)$

where sgn (x_j − x_i), is equal to +1, 0, or −1. When S is greater than 0, x_j > x_i, an upward trend is indicated, while negative S values, xj < x_i, indicate a downward trend. If x_j − x_i is 0, trend is not present within the time series. The test statistic τ can be computed as:

(5)$\tau =\frac{S}{\mathrm{n}\left(\mathrm{n}-1\right)/2}$

which has a range of −1 to +1 and is analogous to the correlation coefficient in regression analysis. S and τ values that are significantly different from zero indicate the presence of trend and reject the null hypothesis of no trend. Following that, the rate of change or the power of the trend can be obtained using Sen’s slope estimator [78]:

(6)${\beta }_1=\mathrm{median} \left(\frac{x_{j }-x_i}{y_{j }-y_i}\right)$

in which the slope of all the pairs is computed and the median of these slopes is denoted.

Multiple Linear Regression

We used multiple linear regression to account for the association of relevant variables including mean summer temperature (May–September) and winter precipitation (October–April), with the annual mass balance response variable. The regression model that is estimated specified annual mass balance (B_a) as a function of time (t), summer temperature (T_s), and winter precipitation (P_w) is as follows:

B_a = β₀ + β₁ t + β₂ T_s + β₃ P_w +ε(7)

where β₀ is the intercept, ε is the unexplained model error, and β_1..3 refers to the slope coefficients of time, summer temperature, and winter precipitation, respectively. The null hypothesis related to time is that the parameter β₁ = 0. Rejecting this null hypothesis in favor of a nonzero β₁ would indicate a trend. Similarly, the null hypotheses related to parameter β₂ and β₃ are that both are equal to zero (i.e., no trend). Evidence against these hypotheses (i., β₁ > 0) suggests that they are helpful in understanding temporal variation in annual mass balance. If the null hypotheses are not rejected, however, a simple regression model should be applied. The key assumptions in this model include normal distribution of residuals, independent variables, and homoscedasticity, which are tested and satisfied, justifying use of this technique.

4.2. Remotely Sensed Data Analysis

In the following Section 4.2.1 and Section 4.2.2 we use remotely sensed data to link the mass balance analysis and the local and regional climate forces with the hypsometry and decadal frontal retreat of glaciers to confirm the statistical results from the previous sections.

4.2.1. Hypsometric Index (HI) and Hypsometric Integral (HIn)

We use Equation (8) as suggested by Jiskoot et al. [81] to perform the HI calculations.

HI = (Hmax − Hmed)/(Hmed − Hmin), If 0 < HI < 1 then HI = −1 / HI(8)

Here, H_max and H_min are the maximum and minimum elevations, respectively, and H_med is the median elevation. Based on the thresholds suggested by Jiskoot et al. [81], glaciers can be divided into five categories: (1) very top-heavy (HI < −1.5), (2) top-heavy (−1.5 < HI < −1.2), (3) equidimensional (−1.2 < HI < 1.2), (4) bottom-heavy (1.2 < HI < 1.5), and (5) very bottom-heavy (HI > 1.5).

HIn, also known as the ‘elevation–relief ratio’, represents the area under the normalized elevation curve with values between 0 and 1. HIn is widely used in hydrological and geomorphological studies and is calculated using Equation (9) [82,83].

HIn = (H_mean − H_min)/(H_max − H_min)(9)

Here, H_max and H_min are the maximum and minimum elevations, respectively; and H_mean is the mean elevation. Low HIn values represent a concave hypsometric curve, whereas larger values depict a convex curve. Both HI and HIn are used to give a better understanding of the erosion and fluvial activity along with the shape of the glacier basin [81,83,84].

4.2.2. Glacier Terminus Delineation

The terminus outlines are delineated manually based on visual observations of temporal Landsat images. The approach defined by Bjørk et al. [85] is adopted for reporting the recession rates, where the ice front is first divided by generating equally spaced points. We then calculate the mean distance between points with respect to a fixed point placed up glacier. To maintain uniformity in this approach related to data spatial resolution, we make all the observations on 30 m Landsat images. We follow the approach suggested by Hall et al. [86] and use Equation (10) to estimate the uncertainty associated with terminus recession rates:

U = ((I1)²+(I2)²)^1/2 + E_reg(10)

where U is the total uncertainty in terminus position, I1 is the pixel size of the image from earlier year, I2 is the pixel size of the image from the later year, and E_reg the co-registration error between images. The geo-registration of Landsat Tier 1 scenes is consistent and within prescribed image-to-image tolerances of ≤12 m radial root mean square error (RMSE). Thus, considering 12 m as the E_reg value provides us with a total uncertainty of U~ 54.43 m.

5. Results and Discussion

5.1. Analysis of Annual Mass Balance

Trend calculations on the annual mass balance dataset for the selected glaciers in the Swiss Alps, from the early 1900s to 2017, reveal a negative downward trend and two significant breaks in the 1980s and 2000s (Figure 2b and Table 3). However, the Ba of Gries and Schwarzberg glaciers show three break points. Accordingly, the series of Schwarzberg includes three phases: (i) the years from 1956 to 1970 with a near-equilibrium condition (μ_Ba = −69 mm w.e., σ = 344 mm w.e.), (ii) the period from 1971 to 1980 of mass gain (μ_Ba = 621 mm w.e., σ = 933 mm w.e.), and (iii) the period after 1980 with a strong imbalance due to mass loss (μ_Ba = −640 mm w.e., σ = 650 mm w.e.). Similarly, Gries shows breaks in: (i) the period from 1962 to 1980 with a nearly balanced condition, (ii) from 1981 to 2000 (μ_Ba = −885 mm w.e., σ = 541 mm w.e.), and (iii) a rapid mass loss period from 2001 to 2017 (μ_Ba = −1527 mm w.e., σ = 590 mm w.e.). Among the other glaciers, the rates of change after break years (1980s) become negative for Allalin and Hohlaub. In Gietro and Silvretta, the rates have changed from −200 to −1000 mm w.e., indicating a trajectory toward the disappearance of these glaciers in the near future. Based on these results, we can conclude that the 1970s was a decade of mass gain, and the mid-1980s onward, especially the 2000s and 2010s, have been the decades of rapid mass loss for many of the glaciers in the Swiss Alps. According to [87], the melt extremes of the 2000s occurred under high shortwave irradiance and large downward latent and sensible heat fluxes in association with increased atmospheric moisture content. In contrast, melt events at the end of the 1980s were associated with lower downward turbulent heat fluxes explained by high downward shortwave radiation under relatively dry atmospheric conditions.

The correlation analysis of glaciers reveals a significant relationship among the Ba of the glaciers (Figure 3). Accordingly, Silvretta, Clariden, Gries, and Gietro can be categorized into one group, and Allalin, Schwarzberg, and Hohlaub in another group. The similarity of Allalin, Schwarzberg, and Hohlaub can be attributed to the proximity of these glaciers and corresponding geomorphic (e.g., elevation range, degree of debris coverage, surface slope, and surface facing) and local climatic factors (e.g., precipitation and temperature).

Although the glaciers in the first group are located apart from each other, they have exhibited similar rapid rates of change in mass balance in the available series. In this group, the size of these glaciers, except Gietro, is less than 5 km² (see Table 1), which makes them more sensitive to warming air temperatures.

5.2. Analysis of Seasonal Mass Balance and Accumulation Area

The analysis of Bs and Bw suggests that the changes in mass balance are primarily influenced by increases in ablation rates during summer (Table 4), which agrees with previous studies [3,33].

The comparison of summer mass balance and AAR series indicates that the glaciers’ accumulation area has decreased at the expense of increased ablation rates and areas. Among the glaciers, the accumulation area of Gries approached zero from 2000 to 2007 and again from 2014 to 2017. Gries is a medium-sized steep glacier that is faced to the northeast. Unlike Gries, the highest increase in AAR was observed in Clariden, which is likely related to the greater amount of solid precipitation received during the accumulation months. Overall, the accumulation areas of larger glaciers have been more consistent than smaller glaciers. For example, the AAR of Aletsch glacier, the largest of those surveyed, has fluctuated around 50% through the entire period and has been less sensitive to climate variations. The larger glaciers’ state can be explained by their ability to dynamically adjust to climate warming by retreating to higher elevations [33]. Further, the decrease in surface albedo and bare terrain exposure at the terminus coincides with an increase in upward longwave radiation from and downward shortwave absorption to accelerate glacier frontal retreat and melt in the ablation area [33,87].

The time series of glaciers is continuous within the common period of 1970 to 2017 (Figure 2b). With the comparison of glaciers in this period, lower Bs mean and higher Bs variation are observed for Gries with −2232 mm w.e. and 838 mm w.e., respectively. During the accumulation season, glaciers received an average of approximately 1000 mm w.e., with the highest value for Clariden (1385 mm w.e.). In the same period, the lowest average annual mass balances of −919 mm w.e. for Gries and −509 mm w.e. for Silvretta were detected. Similarly, the mean AARs were significantly lower for Gries and Silvretta (26% and 36%) and higher for Clariden (56%). The comparison of seasonal and AAR values indicates that the glaciers like Hohlaub, Allalin, and Schwarzberg are roughly in equilibrium relative to other glaciers due to higher altitudinal ranges. Switzerland’s air temperatures are highly dependent on elevation, and much of the observed warming has been found at lower elevations due to increased incoming shortwave radiation under more frequent anticyclonic conditions [88].

5.3. Climatic Controls

5.3.1. Local Meteorological Variables

The overall analysis of days with the maximum temperature > 0 °C per year for all selected glaciers indicates that the total number of days has considerably increased from 130 days in 1970 to 175 in 2015 with a change rate of ~9 days per decade (Figure 4).

The strongest summer warming trends were during the period of 1970 to 2017 (Figure 5b,d,f,h). The largest trends were found for glaciers located in the southern Rhone including Gietro, Hohlaub, Allalin, and Schwarzberg (0.58 °C decade⁻¹) (Figure 5b), and the smallest trend was found for for Silvretta (0.46 °C decade⁻¹) (Figure 5f). Following the same pattern of the Italian Alps [3], three discernible phases can be detected from a 11-year running mean applied to the time series. The first phase starts from 1970 with the mean air temperature anomalies around −1 °C. Then, after 1985, the anomalies reached 0 °C and remained stable until 2000s. Since the 2000s, air temperature anomalies over glaciers have been accelerating beyond 0 °C and have approached +1 °C in recent years with the unprecedented peak of roughly +3 °C for all glaciers in 2003.

The warming trends during the winter season vary from 0.32 to 0.43 °C decade⁻¹, with the lowest trend associated with Silvretta and the highest to the southern Rhone glaciers. The fluctuations in temperature values over the selected period were different for the glaciers located in southern Rhone during the accumulation season. The corresponding moving averages show that the positive trend lasted up to 1985, then gradually declined for five years until it reached 0–1 °C above average in 1990s. Over the 2010s, the temperature anomalies of this group of glaciers oscillated between +2 and +3 °C during the winter months. The warmest peaks of these glaciers happened in 1989 and 2014/2015, whereas a dramatic drop (below 1 °C for the southern glaciers and below −1 °C for other groups) occurred in 2010. Despite the dissimilarity of Silvretta and Clariden during summer (Figure 5f,h), the temperature trends are quite similar during the accumulation season.

The total precipitation anomalies do not reveal significant trends in accumulation season (Figure 6a,c,e,g). The moving averages of precipitation have fluctuated close to the long-term mean. However, for the northern Rhone glaciers (Figure 6c), Silvretta (Figure 6e), and Clariden (Figure 6g), two phases are apparent. In the first phase, periods with above-average precipitation were frequent prior to the second half of the 2000s, with the maximum of ≈40% and the minimum of ≈−40% in 1985 and 2005 and the maximum of 50% at the beginning of the 1980s. In the second phase, running average precipitation remained steady prior to 2010 and was then followed by a ≈10% drop below the long-term mean after 2010.

Precipitation during the ablation season showed a slightly different pattern (Figure 6b,d,f,h), especially for the glaciers located in the southern Rhone (Figure 6b). This group of glaciers has experienced an upward trend over the selected period, including three phases as highlighted by the moving average. The period of the 1970s and the first half of the 1980s with precipitation 10% below the average was followed by a drastic drop to 20% below average that lasted for 20–25 years. From 2000 to 2015, only weak fluctuations in precipitation anomalies occurred as values remained around the long-term mean. The variability of precipitation in summer months is considerably higher in the northern Rhone glaciers (Figure 6d) compared to Silvretta (Figure 6f) and Clariden (Figure 6h). Maximum observed summertime precipitation, 50% and 90% above the mean, occurred in 2002 in the northern and southern Rhone glaciers, respectively, while minimum precipitation, 40% below mean, is evident among all glaciers in 2003. The investigation of local meteorological time series supports the results in Section 5.2 that indicated that increased ablation duration and the associated feedbacks, e.g., surface elevation and albedo reduction, are the main culprits of the annual mass imbalance for the Swiss glaciers.

Reduction in snow/ice-albedo creates a positive feedback mechanism that amplifies snow melt as the temperature increases [87]. Exposed dark mountain surfaces absorb incoming shortwave radiation, which in turn augments near-surface air temperatures. This finding is supported by the higher significant correlation of Ba with May–September and annual air temperatures than the seasonal and annual precipitation (Table 5 and Table 6). In contrast to this overall pattern, the Ba of Clariden and Silvretta show a higher correlation with May–September precipitation and October–April air temperature. Correlations may vary between glaciers as orographic precipitation and accumulation in higher elevations may largely differ from one glacier to another. Therefore, the altitude differences between glaciers and the corresponding orographic effect, along with wind drift and avalanches from near slopes, appear responsible for the accumulation differences [44,89].

A multiple linear regression model is applied to winter precipitation and summer temperature to assess the combined impact of these two components on annual mass balance. The results indicate that the winter precipitation effect is weaker than the May–September temperature effect (Table 7), meaning that at a given precipitation level, the temperature has a stronger, more negative association with the annual balance. Among the glaciers, the Ba of Clariden is more highly influenced by winter precipitation (r = 0.27, p < 0.05) than by summer air temperature (r = −0.06), whereas the Ba of Gries and Schwarzberg glaciers have a significantly inverse relationship with May–September temperature, r = −0.40 (p < 0.01) and r = −0.33 (p < 0.05), respectively.

5.3.2. Large-Scale Atmospheric Circulation

The selected glaciers and the two other glaciers in the central and eastern Alps, Clariden and Silvretta, showed regional mass loss consistency throughout the 20th century. Therefore, it is important to further examine the relationships between seasonal components of glaciers mass balance and large-scale North Atlantic atmospheric circulation and ocean surface temperature conditions that impact the temperature and precipitation regimes of the glaciers.

The correlation analysis indicates that there is predominantly a positive relationship between winter mass balance and seasonal and annual teleconnection patterns (Table 8). Although the corresponding values are not highly significant, there are negative correlations between winter mass balance and winter GBI and the summer NAO_Hurrell index. The correlations between summer balance and the selected indices are mostly significant (p < 0.05). A strong relationship is found between the annual AMO and EA indices and summer mass balance with correlation coefficients of −0.86 and −0.88 (p < 0.05), respectively. The ten lowest Bs had summer AMO and EA mean values of 0.22 and 0.46, respectively. Further, the summer GBI index is significantly anticorrelated with summer mass balance patterns in Switzerland. The summer GBI tends to direct the North Atlantic storm track and associated moisture over Northern Europe and concurrently influence a warming and drying effect over much of the central/southern continent [53]. In contrast to other indices, both summer NAO_CPC and SCA show a positive correlation with summer mass balances, r = 0.69 and 0.66 (p < 0.05), respectively. However, the lagged associations between the teleconnections and summer and annual mass balance should be taken into account, as cold season precipitation (presumably snow) influences the melt rates of the following summer.

In addition to the correlation between summer and winter balances, we also examined the relationship between annual balance and teleconnection patterns. In general, seasonal and annual AMO, GBI, and EA show statistically significant, negative correlations with the annual balance of the Swiss glaciers. In comparison with the other indices, we conclude that the summer and annual balance of the Swiss glaciers are largely driven by AMO and EA behaviors (Table 8). Figure 7 highlights this strong negative association between summer mass balance and summer AMO (r =−0.86) and annual EA (r =−0.88), as well as the consistency of EA and AMO staying above 0 °C since the late 1990s. This is well corroborated by Huss’s (2010) findings, showing that the decadal AMO and mass balance of Swiss glaciers exhibit a strong and negative correlation of 0.78 [32]. Our findings similarly suggest that the positive (warmer-than-average) North Atlantic sea surface temperature (SST) anomalies amplify glacier mass losses (Figure 7), perhaps through warmer regional air temperatures and increased rain versus snow at high altitudes. In contrast, negative (cool) SST anomalies are associated with periods of cooler air temperatures and more snowfall, resulting in positive mass balance.

The oceanic surface temperature conditions are linked to the mass balance of the Swiss glaciers through atmospheric circulation. For instance, under a positive NAO regime, low pressure near Iceland and high pressure over the Azores work in tandem to promote a zonal jet stream and westerly winds that transport warm, moist air from the North Atlantic Ocean to continental Europe.

Comas-Bru and McDermott (2014) noted the association of positive NAO, EA, and SCA indices with below-normal winter precipitation and colder winters in southern Europe [90]. The authors also examined the relationship between the indices and found that winter precipitation patterns are affected by the state and combination of the NAO with either the EA or the SCA; NAO and EA patterns of the same sign coincide with poleward migration of the Icelandic Low and Azores High over North Atlantic waters. When positive phasing occurs between the patterns, the northward shifted jet stream leads to decreased precipitation over central and southern Europe. On the other hand, in-phase NAO and SCA patterns tend to coincide with low pressure atop Greenland with the subtropical high situated over Western Europe. This meridional jet stream configuration produces southwesterly flow and increased precipitation across Scotland and northern Scandinavia while drier conditions persist across much of central Europe [90].

These spatial precipitation patterns, in part modulated by the ocean–atmosphere climate modes, are supported by the strong anticorrelation found herein between the EA pattern and Swiss glaciers’ mass balance characteristics (e.g., positive EA tends to yield less precipitation to the region). Future work will further explore observed seasonal mass balance responses to teleconnection phasing and related ocean–atmosphere interactions, including patterns of temperature advection off the mid-to-high latitude North Atlantic Ocean.

5.4. Hypsometric Characteristic and Terminus Retreat

Hypsometric analyses of the studied glaciers provide some new insights (Figure 8). In general, the shape of the hypsometric curve explains the temporal changes in the slope of the original basin. Strahler classified the drainage basins as youth (convex upward curves), mature (S-shaped curves that are concave upward at high elevations and convex downward at low elevations) and peneplain or distorted (concave upward curves) [91].

In addition to describing the stages of the landscape evolution, these hypsometric curves also provide an indication of erosion status for a watershed. HIn, in general, gives an estimate of the distribution of landmass volume above a reference plane [92], and for a glacier, it may represent the geometry and corresponding volumetric distribution of ice and snow. The inset plot in Figure 8 shows a significant congruence between HI and HIn for the selected glaciers, further highlighting the agreeable inferred geometry of the glaciers represented in Table 9.

HI and corresponding shapes are compiled in Table 9. We note that the two DEMs used for this analysis are based on availability and differ in origin and source data, spatial resolutions, accuracy, and processing steps (Figure 9). As described in Section 5.1, the rate of mass balance change after the 1980s became negative for Allalin and Hohlaub. According to Table 9 and Figure 8, both glaciers are bottom-heavy with N-NE facing aspects. In Gietro and Silvretta, the rates also significantly changed from −200 to −1000 mm w.e., indicating heavy mass loss for these glaciers in near future. It is interesting to note that these two are the only equidimensional glaciers in the study area and both have N-NW aspect (Table 9 and Figure 8).

Similarly, the correlation analysis of Ba presented in Section 5.1 reveals two distinct classes. Silvretta, Gries, and Gietro can be categorized into one group as per correlations (Figure 3) and equidimensional morphologies. Allalin, Schwarzberg, and Hohlaub show a stronger correlation with each other, and interestingly, all three glaciers display bottom-heavy or very bottom-heavy morphologies. Aletsch is the only glacier with ELA below the mean elevation [26,33], and the fluctuation of this glacier near mean is well-supported by its very top-heavy morphology (Table 9). The dashed curves in Figure 8 represent the two glaciers outside the catchment, and the remarkably distinct morphologies of these two glaciers are apparent in the shape of the curves.

Silvretta, in particular, is located in a very narrow elevation range of 672 m, and its hypsometric curve is nearest to a straight line, highlighting its equidimensional shape. The three glaciers with the highest AAR (Table), i.e., Aletsch, Rhone, and Schwarzberg, have similar S-shaped curves that are distinct from the rest of the glaciers (Figure 8). A comparative analysis between Clariden and Gries reveals another interesting facet. These two glaciers are within ~2400–3300 m, have similar area, and have NE facing aspect (Table 1). Clariden displays a convex curve, making it very bottom-heavy, while Gries has a distinct concave profile for the majority of its middle elevations, thus making it an equidimensional glacier (Figure 8). During the accumulation season, Clariden on average showed the highest accumulation (1385 mm w.e.), while in the same period, the lowest average annual mass balance of −919 mm w.e. was observed for Gries. Similarly, the mean AAR was the lowest for Gries (26%) and the highest for Clariden (56%) from 1970 to 2017. The very bottom-heavy morphology of Clariden, its gentler slope than Gries (Figure 9), and significant winter accumulation are well coupled with Ba and Bw (Table 4). Thus, these observations highlight the observed correlations and mass balance trends in agreement with hypsometric distribution and morphology of glaciers. Future morphological changes due to continued mass loss in these glaciers have the potential to further accelerate the rate of negative mass balance.

We further observed the effect of local climate controls and hypsometry on the glacier termini changes over the past three decades (1984–2014) (Figure 10 and Table 10). The frontal retreat patterns reveal a significant spatiotemporal disparity (Table 10). Apart from Allalin and Silvretta, all the glaciers show higher retreat rates in the recent decade (2004–2014). The retreat rate is even stronger for the two bottom-heavy glaciers in the Mattmark region, i.e., Hohlaub and Schwarzberg, at ~50% and ~647%, respectively, higher than the previous decade (Figure 10, Table 10). The dramatic retreat of Schwarzberg can be attributed to a substantial mass imbalance since the 1990s (μBa = −640 mm w.e., σ = 650 mm w.e., rate of change = −641 mm w.e.) that reciprocated in its extreme terminus retreat rate in the next decade. As mentioned above, the three glaciers with the highest AAR (Table 1), i.e., Aletsch, Rhone, and Schwarzberg, have smoother S-shaped curve that are distinct from the rest of the glaciers (Figure 8).

Table 10 reveals that Rhone has also followed the same substantial retreat path of Schwarzberg in the recent decade (2004–2014) with the retreat rate increasing by ~424% compared to the previous decade. It is notable that the decadal retreat of Aletsch was considerably higher (>1100 m) from 1984 to 2014 (Table 10). WGMS report also confirms the accelerated thinning of this glacier since the late 1980s [26]. Being the largest glacier in the Alps, it indicates varying climatic controls over its accumulation and ablation zones. Although the local climate of Aletsch’s front is characterized as dry, the accumulation area receives a high amount of precipitation due to regional advection effects [49].

6. Conclusions

This study effectively used the available glaciological, terrain, and remote sensing data to provide updated analyses of long-term annual and seasonal mass balance time series of glaciers across the Swiss Alps. The purpose of this analysis was to identify trends/changes in the mass balance series and study the relationship between the components of mass balance, and local climatic variables and large-scale oceanic and atmospheric patterns. The most prominent result of this research is the significant glacier mass loss and acceleration in ablation due to the warming signals influenced by North Atlantic climate patterns. According to the results emerging from this analysis, the following conclusions can be reached:

• -. All selected glaciers have lost their equilibrium condition in the recent century, and persistent negative annual mass balance trends, accompanied by decreasing AARs, have been observed for the glaciers located in the south of Switzerland.

• -. Such imbalanced behavior is the product of ablation induced by warmer temperatures, increased melting, and acceleration of the ice-albedo feedback process. Although the precipitation trends are not significant, the increase in the number of days with the air temperature above 0 °C supports the reduction in solid precipitation over the glaciers. Additionally, rain-on-snow events may also be another contributor to the ice mass loss in these regions.

• -. The analyzed glaciers reveal that the annual mass balances are mainly controlled by the summer mass balance, which is statistically attributed to the significant anticorrelation with summer air temperature.

• -. Atmosphere–ocean teleconnection patterns, including AMO and EA, are strongly linked through the time to mass balance characteristics of the southern Swiss Alps, and these large-scale climatic forcings will likely continue to influence surface mass balance regime.

• -. Most of the observed glaciers are losing accumulation area given that AAR percent has dropped from approximately 75% in 1970 to 25% in recent years, which highlights the acceleration of the melt rate.

• -. Although a decadal frontal retreat analysis during 1984–2014 reveals significant spatiotemporal variations, nearly all glaciers have considerably receded during 2004–2014. This highlights the possibility of continuing high melt rates in the next decades and corroborates the mass balance–climate coupling in the region.

• -. Based on the current research and the corresponding results, the negative mass balance trends cannot be simply attributed to the direct response of global warming. Hence, we emphasize that future investigations should focus on resolving complex and indirect surface–atmosphere interactions (e.g., reflectivity and impurities of glacier surfaces [93]) to more comprehensively understand physical processes driving glacier mass balance variability. In particular, the hypsometric evolution of glaciers with continued mass loss warrants considering dynamic modeling approaches to simulate the future state of glaciers as we observe glacier hypsometry and morphology playing considerable roles in defining mass balance trends. In an attempt to better understand the physical drivers of the Swiss glaciers’ mass balance, we recommend that future work provide a deeper evaluation of the spatiotemporal impacts of climatic teleconnections on the glaciers, including analyzing patterns of regional temperature advection along with wind speed and direction variations as it relates to changing glacial conditions.

Author Contributions

Conceptualization, S.G., J.L.R.J., and T.J.B.; methodology, S.G., A.B., and T.J.B.; software, S.G. and L.S.; validation, S.G.; formal analysis, S.G., A.B., and L.S.; investigation, S.G.; resources, S.G. and T.J.B.; data curation, S.G.; writing—original draft preparation, S.G. and A.B.; writing—review and editing, S.G., J.L.R.J., T.J.B., D.R.B., R.C.W., and A.B.; visualization, S.G., L.S., and A.B.; supervision, J.L.R.J., T.J.B., R.C.W., and D.R.B.; funding acquisition, S.G. and J.L.R.J. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Texas State University through the Doctoral Research Support Fellowship. Support for T.J.B. was provided by the University of Alaska Fairbanks Experimental Arctic Prediction Initiative.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The monthly CPC teleconnections (NAO, EA, SCA) are available at the CPC website. The NAO (Cropper daily version) can be obtained from Thomas Cropper, the NAO (Hurrell version) from Climate Data Guide, and the daily GBI [53] from NOAA. The AMO monthly data is also available at NOAA.

Acknowledgments

Weather data were purchased from MeteoSwiss and delivered by Gergely Rigo. We thank Huss for sharing annual and seasonal mass balance datasets, and Klimeš for revising the earlier manuscript draft.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

AMO	Atlantic Multidecadal Oscillation
AARs	accumulation area ratios
GBI	Greenland Blocking Index
EA	East Atlantic pattern
WGMS	world glacier monitoring service
GTNG	global terrestrial network for glaciers
NAOEA	North Atlantic Oscillation
GLOFs	glacial lake outburst floods
ASTER	advanced spaceborne thermal emission and reflection radiometer
GDM	global digital elevation model
SCA	Scandinavian pattern
VAW/ETHZ	laboratory of hydraulics, hydrology, and glaciology at ETH Zürich
ELA	equilibrium line altitude
Ba	annual mass balance
Bw	winter mass balance
Bs	summer mass balance
NBCN	Swiss national basic climatological network
NOAA	National Oceanic and Atmospheric Administration
CPC	NOAA Climatic Prediction Center
EOF	empirical orthogonal function
CUSUM	cumulated sum of standard residuals
efp	empirical fluctuation process
BIC	Bayesian information criterion
ACF	autocorrelation
PACF	partial autocorrelation
HI	hypsometric index
HIn	hypsometric integral
RMSE	root mean square error

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Overview maps of the research study area and the glaciers with mass balance measurements analyzed in this work (source: ESRI, USGS, NOAA, and DHM25/200). The map represents: (a) Silvretta (b) Clariden (c) Rhone (d) Gries (e) Grosser Altesch (f) Hohlaub, Allalin, Schwarzberg (g) Gietro.

Figure 2. Data coverage of annual mass balances (Ba) (a) and the time series of annual mass balance (Ba), summer mass balance (Bs), winter mass balance (Bw), and accumulation area ratios (AAR) for the selected glaciers from 1914 to 2017 (b) (Source: [33]).

Figure 3. A heatmap representation of the Spearman correlation matrix of detrended Ba and dendrogram of similar groups of glaciers (group 1: Silvretta (Silv.), Clariden (Clar.), Gries (Gri.) and Gietro (Gie.); group 2: Allalin (Alla.), Schwarzberg (Sch.) and Hohlaub (Hoh.)). The right bar shows the color code of the correlation coefficients (all significant at p < 0.01).

Figure 4. Number of days with the temperature above 0 °C over all glaciers per year. The blue and red shades show the dominant pattern of cold and warm periods below and above 130 days, respectively.

Figure 5. Winter air temperature anomalies are shown for (a) the southern Rhone (including Gietro, Hohlaub, Allalin, and Schwarzberg glaciers), (c) the northern Rhone (including Grosser Altesch, Gries, and Rhone glaciers), (e) Silvretta, (g) Clariden. Summer air temperature anomalies are presented in (b) the southern Rhone, (d) the northern Rhone, (f) Silvretta, (h) Clariden. The blue and red shades show significant changes.

Figure 6. Winter total precipitation anomalies in (a) the southern Rhone, (c) the northern Rhone, (e) Silvretta, (g) Clariden. May–Sept (summer) total precipitation anomalies in (b) the southern Rhone, (d) the northern Rhone, (f) Silvretta, (h) Clariden. The blue and red shades show significant changes.

Figure 7. The relationship between the annual East Atlantic (EA) (a) and summer Multidecadal Atlantic Ocean (AMO) (b) indices versus summer mass balance from 1970 to 2017, with strong negative correlation coefficients of −0.86 and −0.88 (p < 0.05), respectively, for the 38-year record.

Figure 8. Comparable hypsometric profiles for selected glaciers. The colored boxes for each profile mark the ELA of the respective glacier. The dashed curves represent the two glaciers outside the catchment. The inset plot represents the relationship between HI and HIn for the studied glaciers.

Figure 9. The elevation maps of the selected Swiss glaciers (Sources: swissALTI3D, and ASTER GDEM V3 for Clariden and Silvretta).

Figure 10. Glacier terminus changes during 1984–2014.

Table 1

Characteristics of the selected Swiss glaciers referring to the end of the observation period. Glaciers are ordered based on the time series length (Δt), elevation range (elev. range), equilibrium line altitude (ELA), area, accumulation area ratio (AAR), and the prevailing aspect. Rhone and Aletsch are the two discontinuous time series, indicating seasonal data are not available for the entire period (Source: [33]).

Glacier	Location (Alps)	Period	Δt (yr)	Elev. Range (m a.s.l.)	ELA (m a.s.l.)	Area (km2)	AAR (%)	Prevailing Aspect
Allalin	Valais	1956–2017	62	2693–4180	3575	9.64	27	NE
Hohlaub	Valais	1956–2017	62	2843–4030	3365	2.128	30	E
Schwarzberg	Valais	1956–2017	62	2680–3566	3155	5.10	33	NE
Gries	Valais	1962–2017	56	2432–3307	3285	4.34	0	NE
Gietro	Valais	1966–2017	52	2751–3817	3315	5.20	22	NW
Rhone	Bernese	1885–2017	39	1775–3620	3085	15.10	35	S
Aletsch	Bernese	1940–1999	59	1560–4085	3055	83.02	56	SE
Silvretta *	Central	1918–2017	99	2474–3071	3025	2.60	1	W
Clariden *	Glarus	1914–2017	103	2534–3251	2935	4.50	25	NE

(* two glaciers outside the study area).

Table 2

Overview of Landsat data used for termini delineation.

Year	Date of Capture	Landsat Product ID	Sensor	Glaciers Covered
1984, 1985	2 September 1984	LT05_L1TP_194028_19840902_20170220_01_T1	Landsat 5 Thematic Mapper (TM)	Aletsch, Allalin, Clariden, Gries, Hohlaub, Rhone, Schwarzberg
	2 September 1984	LT05_L1TP_194027_19840902_20170220_01_T1	Landsat 5 TM	Silvretta
	12 September 1985	LT05_L1TP_195028_19850912_20171212_01_T1	Landsat 5 TM	Gietro
1994	4 August 1994	LT05_L1TP_195027_19940804_20180215_01_T1	Landsat 5 TM	Clariden
	4 August 1994	LT05_L1TP_195028_19940804_20180215_01_T1	Landsat 5 TM	Aletsch, Allalin, Gietro, Gries, Hohlaub, Rhone, Schwarzberg
	6 August 1994	LT05_L1TP_193028_19940806_20180214_01_T1	Landsat 5 TM	Silvretta
2003, 2004	14 September 2003	LT05_L1TP_195028_20030914_20161204_01_T1	Landsat 5 TM	Gietro
	2 September 2004	LT05_L1TP_193028_20040902_20161129_01_T1	Landsat 5 TM	Silvretta
	9 September 2004	LT05_L1TP_194028_20040909_20161130_01_T1	Landsat 5 TM	Aletsch, Allalin, Clariden, Gries, Hohlaub, Rhone, Schwarzberg
2014	14 September 2014	LC08_L1TP_193028_20140914_20170419_01_T1	Landsat 8 Operational Land Imager (OLI)	Silvretta
	28 September 2014	LC08_L1TP_195027_20140928_20170419_01_T1		Clariden
	28 September 2014	LC08_L1TP_195028_20140928_20170419_01_T1	Landsat 8 OLI	Aletsch, Allalin, Gietro, Gries, Hohlaub, Rhone, Schwarzberg

Table 3

Structural break points of Ba and the pre/post-break change coefficients (mm w.e.).

Glacier	Breaks	Chow Test	Rate of Change (Seg1)	Rate of Change (Seg2)	Rate of Change (Seg3)
Allalin	1980	0.01	242.8	−505.4	-
Hohlaub	1983	0.02	140.03	−647.3	-
Schwarzberg	1971, 1980	0.40, 0.00	−80.2	686	−641.2
Gries	1980, 2001	0.04, 0.30	−91.6	890	−1561.4
Gietro	2001	0.20	−224.5	−874.8	-
Silvretta	2001	0.00	−253	−1076.6	-
Clariden	1988	0.01	17.3	−605.5	-

Table 4

Correlation coefficients of Ba versus Bw and Bs. Values in bold are significant at p < 0.01.

Glacier	Δt	Bw	Bs
Allalin	62	0.42	0.90
Hohlaub	62	0.48	0.89
Schwarzberg	62	0.33	0.89
Gries	56	0.12	0.88
Gietro	52	0.15	0.85
Rhone	39	0.44	0.90
Aletsch	59	0.39	0.92
Silvretta	99	0.30	0.87
Clariden	103	0.46	0.88

Table 5

Correlation coefficients of summer, winter, and annual precipitation versus Ba. Superscripts ^a–c indicate significance at p < 0.01, 0.05, and 0.10, respectively.

	Allalin	Gietro	Gries	Hohlaub	Schwarzberg	Silvretta	Clariden
Winter	0.15	0.15	0.28 c	0.05	0.11	0.22	0.27 b
Summer	0.01	−0.20	−0.16	−0.07	−0.14	−0.37 b	−0.43 a
Annual	0.16	0.01	0.13	−0.03	−0.01	0.00	0.00

Table 6

Correlation coefficients of summer, winter, and annual air temperature versus Ba. Superscripts ^a–c indicate significance at p < 0.01, 0.05, and 0.10, respectively.

	Allalin	Gietro	Gries	Hohlaub	Schwarzberg	Silvretta	Clariden
Winter	−0.16	−0.19	−0.43 a	−0.18	−0.29 b	−0.32 a	−0.35 b
Summer	−0.20	−0.23	−0.40 a	−0.22	−0.33 b	−0.21	−0.06
Annual	−0.25 c	−0.25 c	−0.51 a	−0.27 c	−0.40 a	−0.36 a	−0.27 c

Table 7

Multiple regression coefficients for winter precipitation and summer temperature explaining B_a variation of Swiss glaciers. Superscripts ^a–c indicate significance at p < 0.01, 0.05, and 0.10, respectively.

	Allalin	Gietro	Gries	Hohlaub	Schwarzberg	Silvretta	Clariden
Winter precipitation	96.69	130.90	116.80	62.70	102.34	131.95	167.76
Summer temperature	−174.30 b	−123.60	−325.32 a	−231.50 c	−317.78 a	−175.38	−58.68
Intercept	−318.00	−648.60	−409.60	−89.46	567.42	−500.80	−776.65
% of explained variability	10.00	6.00	23.00	9.00	19.00	9.00	5.00

Table 8

Spearman correlation coefficients of Bw, Bs, and Ba versus seasonal and annual indices. Superscripts ^a–b specify the significance level of correlations at p < 0.05 and p < 0.10, respectively.

Index	Bw			Bs			Ba
	Winter	Summer	Annual	Winter	Summer	Annual	Winter	Summer	Annual
AMO	−0.08	0.01	−0.05	−0.80 a	−0.86 a	−0.86 a	−0.79 a	−0.82 a	−0.85 a
NAOCPC	0.30 a	−0.16	0.26 b	−0.41 a	0.69 a	0.06	−0.26 b	0.61 a	0.17
NAOHurrell	0.35 a	−0.25 b	0.30 a	−0.15	0.39 a	−0.04	−0.02	0.31 a	0.06
GBI	−0.34 a	0.13	−0.15	0.33 a	−0.68 a	−0.51 a	−0.44 a	−0.61 a	−0.55 a
SCA	−0.03	−0.21	−0.17	0.33 a	0.66 a	0.56 a	0.25 b	0.59 a	0.47 a
EA	0.14	0.03	0.07	−0.86 a	−0.82 a	−0.88 a	−0.79 a	−0.77 a	−0.81 a

Table 9

HIn, HI and glacier categories.

Glaciers	HIn	HI	Glacier Geometry
Aletsch	0.58	−1.62	very top-heavy
Allalin	0.46	1.26	bottom-heavy
Gietro	0.51	−1.16	equidimensional
Gries	0.50	−1.19	equidimensional
Hohlaub	0.36	2.06	very bottom-heavy
Rhone	0.54	−1.24	top-heavy
Schwarzberg	0.44	1.41	bottom-heavy
Clariden	0.32	3.11	very bottom-heavy
Silvretta	0.51	−1.06	equidimensional

Table 10

Glacier retreat rates during 1984–2014.

Glacier	Period	Total Retreat (m)	Retreat Rate (m/yr)
Aletsch	1984–1994	433.09	43.31
	1994–2004	282.11	28.21
	2004–2014	451.03	45.10
Allalin	1984–1994	190.97	19.10
	1994–2004	461.70	46.17
	2004–2014	190.06	19.01
Gietro	1985–1994	102.59	11.40
	1994–2003	78.08	8.67
	2003–2014	287.81	26.16
Gries	1984–1994	85.06	8.51
	1994–2004	209.22	20.92
	2004–2014	228.60	22.86
Hohlaub	1984–1994	94.18	9.42
	1994–2004	90.58	9.06
	2004–2014	136.15	13.61
Rhone	1984–1994	98.81	9.88
	1994–2004	38.56	3.86
	2004–2014	202.17	20.22
Schwarzberg	1984–1994	15.43	1.54
	1994–2004	30.86	3.09
	2004–2014	230.67	23.07
Clariden	1984–1994	137.25	13.72
	1994–2004	170.66	17.07
	2004–2014	193.02	19.30
Silvretta	1984–1994	80.24	8.02
	1994–2004	61.11	6.11
	2004–2014	65.78	6.58