Unravelling the evolution of Zmuttgletscher and its debris cover since the end of the Little Ice Age

Debris-covered glaciers generally exhibit large, gently sloping, slow-flowing tongues. At present, many of these glaciers show high thinning rates despite thick debris cover. Due to the lack of observations, most existing studies have neglected the dynamic interactions between debris cover and glacier evolution over longer time periods. The main aim of this study is to reveal such interactions by reconstructing changes of debris cover, glacier geometry, flow velocities, and surface features of Zmuttgletscher (Switzerland), based on historic maps, satellite images, aerial photographs, and field observations. We show that debris cover extent has increased from ∼13 % to ∼32 % of the total glacier surface since 1859 and that in 2017 the debris is sufficiently thick to reduce ablation compared to bare ice over much of the ablation area. Despite the debris cover, the glacier-wide mass balance of Zmuttgletscher is comparable to that of debris-free glaciers located in similar settings, whereas changes in length and area have been small and delayed by comparison. Increased ice mass input in the 1970s and 1980s resulted in a temporary velocity increase, which led to a local decrease in debris cover extent, a lowering of the upper boundary of the ice-cliff zone, and a strong reduction in ice-cliff area, indicating a dynamic link between flow velocities, debris cover, and surface morphology. Since 2005, the lowermost 1.5 km of the glacier has been quasi-stagnant, despite a slight increase in the surface slope of the glacier tongue. We conclude that the long-term glacier-wide mass balance is mainly governed by climate. The debris cover governs the spatial pattern of elevation change without changing its glacier-wide magnitude, which we explain by the extended ablation area and the enhanced thinning in regions with thin debris further up-glacier and in areas with abundant meltwater channels and ice cliffs. At the same time rising temperatures lead to increasing debris cover and decreasing ice flux, thereby attenuating length and area losses.

1 Motivation and objectives

Debris-covered glaciers have been observed to show a delayed adjustment of their length to climatic changes (e.g. Ogilvie, 1904; Scherler et al., 2011; Banerjee and Shankar, 2013). This behaviour can be explained by melt reduction due to insulation by the debris layer, which commonly increases in thickness towards the terminus (Östrem, 1959; Nakawo et al., 1986; Nicholson and Benn, 2006; Anderson and Anderson, 2018), and is expected to distinctly prolong the glacier's response time (Jóhannesson et al., 1989). Several studies showed that debris-covered glaciers exhibit similar thinning rates to debris-free glaciers in the Himalayas (e.g. Kääb et al., 2012; Gardelle et al., 2012; Nuimura et al., 2012). Explanations proposed for this behaviour include the emergence of surface features with exposed ice (e.g. ice cliffs, water flow channels, ponds), enhanced thinning further up-glacier that compensates for the debris-induced melt reduction, reduced mass flux from the accumulation areas and decreasing emergence velocities at the tongue, and englacial ablation with subsequent roof collapsing (e.g. Pellicciotti et al., 2015; Vincent et al., 2016; Ragettli et al., 2016; Banerjee, 2017; Brun et al., 2018). Ice cliffs and ponds can enhance ablation in comparison to debris-covered surfaces and even debris-free ice and are common features on debris-covered glacier tongues (Benn et al., 2012; Brun et al., 2016). Especially during periods of negative mass balance, the down-glacier increase in debris cover thickness reduces ablation through its insulating effect and can lead to a lower – and even reversed – mass balance gradient (Nakawo et al., 1999; Benn and Lehmkuhl, 2000; Benn et al., 2012; Ragettli et al., 2016; Rounce et al., 2018). Over time, this reduction in ablation can lead to a decrease in surface slope and, consequently, driving stress and ice flow velocity (Kääb, 2005; Bolch et al., 2008b; Quincey et al., 2009; Jouvet et al., 2011; Rowan et al., 2015; Dehecq et al., 2019). Furthermore, with an increase in equilibrium line altitude (ELA) the englacial debris melts out earlier, leading to an extended debris cover further up-glacier (Benn et al., 2012; Kirkbride and Deline, 2013; Carturan et al., 2013). As a consequence of reduced ablation and driving stress, heavily debris-covered glaciers often have long, gently sloping, low-lying tongues with low flow velocities or even stagnant parts (Benn et al., 2012). Current research on many debris-covered glaciers is mostly focussed on processes, such as ablation beneath the debris cover and in areas of ice cliffs and ponds, thinning of glacier tongues, surface flow velocities, or changes in debris cover over time (e.g. Hambrey et al., 2008; Bolch et al., 2012; Benn et al., 2012; Dobhal et al., 2013; Ragettli et al., 2016; Gibson et al., 2017). A few studies have started to integrate the existing understanding and interactions of some of these processes into numerical models for the evolution of such glaciers (Jouvet et al., 2011; Rowan et al., 2015; Anderson and Anderson, 2016; Banerjee, 2017).

However, most studies face difficulties leading to persistent shortcomings in our understanding of the development of debris-covered glaciers: (i) time series are often too short to inform on the full response time of larger glaciers; (ii) investigations are often local because repeated tongue-wide data are sparse, especially for longer periods; (iii) by considering only one or a few variables, the mutual influences of changing debris cover and glacier geometry cannot be assessed; (iv) long-term (> decade) glacier-scale studies have mostly been conducted in the Himalayas (e.g. Bolch et al., 2008b; Ragettli et al., 2016; Lamsal et al., 2017).

To better understand how a changing debris cover affects glacier geometry, flow velocities, and surface features, and how debris cover is in turn affected by these variables, it is necessary to consider the long-term development of glaciers beyond their potential response times. Few studies have investigated the temporal evolution of debris cover on glaciers (overview in Kirkbride and Deline, 2013), or the evolution of debris-covered glaciers over time (e.g. D'Agata and Zanutta, 2007; Capt et al., 2016). Several studies observed an increase in debris cover on glaciers during negative mass balance periods (e.g. Bolch et al., 2008a; Quincey and Glasser, 2009; Kirkbride and Deline, 2013). Because of the overall negative mass balance trend, glaciers are and will be increasingly affected by debris cover. It is important to understand the magnitude of this increase and how it affects the geometry and dynamics of the glaciers in order to simulate their future development (Jouvet et al., 2011; Anderson and Anderson, 2016).

Large debris-covered glaciers in the Himalayas and Karakoram are not ideal candidates for long-term investigations due to their long response times (>50 years because of thick ice and debris cover; Cogley, 2011; Banerjee and Shankar, 2013) and data scarcity before the 1960s. In contrast, the long history of length monitoring as well as the availability of topographic maps and aerial photos from the mid-19th century onwards make glaciers in the Alps and especially in Switzerland well suited to study long-term development. In Switzerland, the earliest maps (1850s, 1870s) already include debris symbols and contour lines, and stereo imagery throughout the 20th century allows for detailed 3-D surface investigations.

In this study we aim to understand the long-term geometric evolution and dynamics of debris-covered glaciers through the study of Zmuttgletscher in the Swiss Alps. This medium-sized valley glacier has been going through the transition from a mostly debris-free glacier in the late 1850s to one that is almost completely debris-covered in its ablation area (2017). We quantify the evolution of geometry (length, area, elevation, slope), mass balance, flow, and debris cover at a high spatial and temporal resolution since the end of the Little Ice Age (LIA) around 1850. We use these data to investigate interactions between geometry evolution, ice flow, and debris cover and the related drivers. We further analyse the occurrence of ice cliffs and their role in long-term glacier evolution.

2 Study site

Zmuttgletscher ($\mathrm{45}{}^{\circ }{\mathrm{59}}^{\prime }$ N, $\mathrm{7}{}^{\circ }{\mathrm{37}}^{\prime }$ E) is located in the southern end of the Matter Valley in the western Swiss Alps and ranges from ∼2240 to ∼4150 m a.s.l. It is surrounded by the Matterhorn (4478 m) and the Dent d'Hérens (4174 m) to the south, the Dent Blanche (4357 m) to the north, and high ridges (above 3400 m) in between (Fig. 1). In 2016, the glacier had a surface area of 15.74 km² and was substantially debris-covered in its ablation area. The debris mainly originates from the surrounding rock walls.

Zmuttgletscher is located in a relatively dry region at the main divide of the Alps; thus it receives precipitation from both northern and southern weather systems. There are no direct measurements at higher elevations, but model estimates suggest values between 0.8 and 1.5 m (MeteoSwiss, 2014). Glaciological mass balance measurements at the nearby debris-free Findelgletscher (15 km distant) suggest end-of-winter accumulation around 0.8–1.5 m water equivalent (w.e.) (Sold et al., 2016). However, at Zmuttgletscher, avalanching additionally contributes to accumulation. When including contributing rock walls and lateral moraines, the total area available for accumulation is up to 22 km² (restricted to areas >30^∘ slope angle, Fig. 1).

Zmuttgletscher has several independent and connected tributaries (Fig. 1). The major accumulation area in the south – Tiefmattengletscher – reaches almost up to the summit of Dent d'Hérens. The western accumulation area – Stockjigletscher – is a relatively flat area above two distinct icefalls between Tête Blanche (3710 m), Stockji (3092 m), and Wandfluehorn (3588 m). Schönbielgletscher is a tributary from the north reaching up to ∼3400 m, below the Dent Blanche summit wall. The upper part of Tiefmattengletscher is completely surrounded by rock walls, including the >1000 m high Dent d'Hérens north wall and the almost 1500 m high Matterhorn west face. At the end of the LIA, the main glacier tongue was nourished from all of these accumulation areas. In 2017, it is mainly fed by Tiefmattengletscher and to a lesser extent from Stockjigletscher. Around 2010 the central, almost debris-free Stockjigletscher branch (between point 2095 and Stockji) detached from the main tongue. Because Stockjigletscher is lacking high rock walls in its accumulation area, it does not deliver substantial amounts of debris down to the ablation area. In contrast, the relatively low-lying Schönbielgletscher receives a lot of accumulation through avalanching from the surrounding rock walls and consequently exhibits a continuous debris cover even above the icefall at ∼2900 m.

Figure 1 (a) Geographical location. (b) Zmuttgletscher, its topographical setting and different tributaries: SBG: Schönbielgletscher; STG: Stockjigletscher; TMG: Tiefmattengletscher. Background image: orthomosaic from aerial imagery by Swisstopo from 2013 (Swisstopo, 2010). (c) Glacier hypsography in the year 2010 (© SwissTopo).

3 Data and methods

Our analyses are based on topographical maps, aeroplane-based and UAV-based aerial images, and satellite images (Table 1), in addition to various field observations, and long-term air temperature measurements. The use of these data is discussed in the respective sections below. Areal images were available from (i) post-war mapping flights by the American military, (ii) national flight campaigns for topographic map production, (iii) specific flight campaigns for glacier monitoring purposes conducted by Swisstopo, (iv) fixed-wing UAV flights (using a Sensefly eBee Classic; Sensefly SA) in 2016 and 2017, and (v) Pleiades tri-stereo images from 2017.

Table 1 Topographic maps and satellite and aerial images used in this study. Abbreviations are as follows: dH: elevation change; DTM: digital terrain model; OP: orthophoto; obl. aer.: oblique aerial; Ple.: Pleiades. Entries denoted with an asterisk are used as a final product. All data except 1894, 2016, and 2017 are from Swisstopo (2018a). Oblique photos from 1894 were taken by Reid (1894).

3.1 Generation of digital terrain models (DTMs) and orthophotos

Glacier surface information was extracted from stereo-photogrammetric DTMs generated from aerial images and the Siegfried map (Table 1). To obtain a DTM representing the glacier tongue in 1879 the Siegfried map from 1880 (Swisstopo, 2010) was georeferenced using ground control points (GCPs). Subsequently, the contour lines were semi-automatically digitised by separating the differently coloured (blue is on-glacier; brown is off-glacier) contour lines from other symbols (Siedler, 2011). Their elevation information was extracted and interpolated using the Topo-to-Raster tool in ArcGIS. The point of origin for elevation measurements, situated at the shore of Lake Geneva, changed in the 1930s from 376.2 to 373.6 m (Swisstopo, 2018a) and the elevations derived from the Siegfried map were corrected accordingly.

The DTMs from 1946 to 2017 (apart from 2010 and the Pleiades DTM from 2017) were created with photogrammetric methods using structure-from-motion software (Agisoft LLC, 2016; Pix4D, 2016). For each date, a point cloud was produced from the available number of input images (4–29) and then georeferenced using a set of 10–15 GCPs, i.e. reference points in the images that could be referred to points on stable ground, taken from Swissimage 2013. The quality of the DTMs is comparable to DTMs from traditional photogrammetric software (see Mölg and Bolch, 2017). Their uncertainty in elevation – defined as the standard deviation over stable terrain around the glacier tongue – mostly lies within 2 m depending on the resolution of the aerial images, the number of images, and image quality factors (somewhat higher values were obtained for the DTMs from 1879 and 1946; Table 1). The uncertainties were derived for terrain with comparable steepness to the glacier surface ($<\mathrm{30}{}^{\circ }$) and are assumed constant in space, although DTM quality decreases in steep areas (e.g. rock walls). The DTM from 2010 (SwissAlti3D) was taken as a final product from Swisstopo and had also been produced from aerial stereo images; it has a nominal resolution and a vertical accuracy of 2 m (Swisstopo, 2018b). In addition, we generated a glacier-wide DTM for 2017 (spatial resolution 1 m) using Pleiades tri-stereo, high-resolution satellite imagery, based on photogrammetric principles using PCI Geomatics (Geomatica, 2016). All DTMs were co-registered to the reference DTM (2010) before further analysis by using the analytical approach by Nuth and Kääb (2011), followed by a second-order trend surface correction to eliminate remaining elevation differences (Pieczonka et al., 2013).

Orthophotos were generated by rectification of the stereo images using the respective DTM.

3.2 Glacier area and length

Glacier areas were measured from the Dufour map (1859), the Siegfried map (1879), all available orthophotos, and the Swissimage (2005, 2007, 2010, 2013; Swisstopo, 2010; Table 1) by manual digitisation. The Dufour map (map sheet 22, section 8, number 485) from 1859 was the first map containing elevation information (in the form of contour lines) and distances acquired with modern methods (Wolf, 1879; Graf, 1896; Rastner et al., 2016). The extent of the glacier and the supraglacial debris could be extracted from the map, whereas its elevation information was disregarded due to strong, non-linear, horizontal distortions.

The time series of maps and orthophotos resulted in glacier area values for each corresponding date since 1859. The mapping quality is often lower in debris-covered areas compared to debris-free areas (Paul et al., 2013), but the high resolution of these images allowed correct interpretation of the glacier margin. The glacier boundary in the accumulation area to the west was taken from the 2010 ice divide and was kept constant over time. The hanging glaciers at the north face of Dent d'Hérens have been included in the glacier area since they contribute mass to the main glacier through frequent ice avalanches.

Front variation measurements were conducted by using the glacier outlines for each date. Along the ice front – perpendicular to the flow – the change was measured at distances of 100 m and then averaged (Koblet et al., 2010; Bhambri et al., 2012). For the comparison to other glaciers we used GLAMOS length variations (GLAMOS, 2018), which were acquired in the field with the same concept of several parallel measurements, equidistant by 50 m. At Zmuttgletscher, GLAMOS measurements have regularly been conducted between 1892 and 1997, but with serious data issues before 1946. The GLAMOS data for Zmuttgletscher since 1946 are almost identical to our own relative length change record (Fig. S1 in the Supplement) and are therefore not further used. To interpret the observed length changes and the potential influence of debris cover, we compared the variations from our measurements to the GLAMOS data of 10 other Swiss glaciers.

The uncertainty of both area and length results is estimated to lie within $\pm \mathrm{1}/\mathrm{2}$ pixel (Table 1) along the glacier boundary and the start and end of the length profile, respectively (Hall et al., 2003; Granshaw and Fountain, 2006; Bolch et al., 2010; Bajracharya et al., 2015). For the calculation of the changes, the uncertainties of the two respective dates are combined by error propagation, analogous to the mass balance data (Hall et al., 2003; Zemp et al., 2013).

3.3 Debris cover, on-site ablation, and on-site air temperature

Debris cover extent for the orthophotos and historic maps was manually digitised (Fig. 2a, b). We mapped only areas of continuous debris cover, thus avoiding the inclusion of sparse debris cover. The debris extent of the Siegfried map (1879) was verified using two photographs taken in 1894 (Fig. 2c). This information was valuable to limit the debris extent at and up-glacier of the confluence of Tiefmattengletscher and Schönbielgletscher, which was the region of the strongest changes (Fig. 2d).

Figure 2 Example of mapped debris cover on (a) Dufour map, 1859 (© SwissTopo), and (b) Siegfried map, 1879 (© SwissTopo). (c) Parts of the tongue in 1894 towards the south-west (Reid, 1894; National Snow and Ice Data Center, 2002, updated 2015). (d) Similar extent as in (c) in 2013 (Swisstopo, 2010). Panels (c) and (d) show the main tongue and Tiefmattengletscher.

Ablation was measured at seven points on the glacier tongue (Fig. 1) during summer 2017 to better understand the influence of debris on melt and to complement the long-term elevation change data. These point measurements were conducted on bare ice and on surfaces with variable debris thickness at elevations between 2300 and 2600 m using 2 m long PVC stakes that were connected by zip ties and drilled into the ice. Debris cover was removed for the drilling and repositioned after inserting the stakes. To get a rough estimate of ablation on ice cliffs, the horizontal backwasting of a south-facing and north-facing ice cliff was measured using horizontally inserted ablation stakes. All stakes were measured in intervals of several weeks over the course of the 2017 ablation season.

Debris thickness data were collected in the field by manual excavation along and in between three transects perpendicular to the glacier flow direction in September 2017 (for transect locations see Fig. 1). Each data point represents the average of three measurements ∼1 m apart, and the standard deviation is used as an uncertainty measure. However, for debris thicknesses greater than 20 cm only one measurement was taken.

Homogenised time series of air temperature measurements by MeteoSwiss (Füllemann et al., 2011) were used to interpret the observed glacier evolution. The selected climate stations had to be as close to the study site as possible, lie at similar elevations, and cover a long period. Since no single station fulfilled all requirements, we used the stations on Col du Grand St. Bernard (2472 m, ∼40 km, 1864–2018) and in Sion (484 m, ∼35 km, 1865–2017).

3.4 Surface flow velocities

Surface flow velocities provide important information about the glaciers' dynamical state and its change over time. Automatic feature tracking methods were not feasible for imagery acquired before 2005 because the time differences and, thus, displacements of the features were too large. Therefore we manually tracked boulders to infer flow velocities along the debris-covered part of the main tongue as well as on the lower debris-covered part of Schönbielgletscher. The tongue was divided into four sections and 11 subsections according to differences in dynamic state and ice flow units (Fig. 8). The individual measurements were averaged for every time period and subsection to achieve a comprehensive picture of the dynamic changes. For the periods 2005–2007 and 2016–2017 we extracted flow velocity fields using the feature tracking module IMCORR in SAGA GIS (Fahnestock et al., 1992; Scambos et al., 1992; Conrad et al., 2015). The results were filtered using a visually defined threshold of correlation quality (“strength” =10.0). Outliers were manually removed, e.g. in the area of strong ice-cliff change or proglacial water surfaces.

Between 22 and 24 August 2017 the tongue of Zmuttgletscher was observed with a terrestrial radar interferometer (GPRI) developed by Gamma Remote Sensing AG (Werner et al., 2008). The GPRI was installed on a hill about 3 km away from the terminus in 2017 (Fig. 1). Measurements were acquired every minute for 1.5 d with a final range resolution of 3.75 m and an azimuth resolution of 7 m at a slant range of 1 km. The interferograms were determined with a standard workflow following Caduff et al. (2015) using the Gamma software, were stacked over a window of 8 h to reduce noise, and were afterwards unwrapped using a stationary point on the ground. The unwrapped phases were then converted to line-of-sight displacements according to Werner et al. (2008), whereby negative velocities were considered to be noise and filtered out. The results were georeferenced by rotating the displacement map to best match with the DEM25 (Swisstopo, 2005). Afterwards, the data were resampled into the new grid using nearest-neighbour interpolation. To assess the uncertainties in the velocities we looked at the difference from zero in measured displacement of 10 stationary points. This results in an uncertainty of the stacked velocity maps of ±0.03 m d⁻¹.

3.5 Surface features

Ice cliffs, exposed ice at supraglacial meltwater channels, and lakes were extracted in a semi-automatic process using an object-based approach with Trimble eCognition (eCognition Essentials 1.3, 2016); see Kraaijenbrink et al. (2016). The location and area of these surface features was determined on all orthophotos from 2013. A primary segmentation divides the image into polygons based on pixel intensity and image texture. Ice cliffs often consist of a lower, steeper section of bare ice and a flatter section of ice covered with sand and pebble-sized debris particles in the upper part, where the maximum solar angle seems to define the slope of pole-facing cliffs (Sakai et al., 2002; Buri et al., 2016). These changing slope areas were often separated by segmentation, and were in a second step manually selected and combined into one polygon per ice cliff or lake. Supraglacial channels cover only small areas and are thus incorporated into the term “ice cliffs”. The described approach is effort efficient and assures a low uncertainty level, which is on the order of $\pm \mathrm{1}/\mathrm{2}$ pixel. In order to assess the potential influence of ice cliffs on surface elevation change, we investigated their changes in area on the glacier tongue for different time periods. A 35 m buffer was generated around the ice-cliff polygon for both date 1 and date 2, and the overlapping area of these two buffer zones defined the category of ice-cliff areas. For each period the average surface elevation change in this category was compared to the surface elevation change at the rest of the tongue (yellow area in Fig. 9b).

3.6 Surface elevation changes and geodetic mass balance

The time series of surface elevation change over the glacier tongue was restricted to the overlapping area of all available DTMs from 1879 to 2017 (11 dates), which reaches up to ∼2750 m a.s.l. and covers >50 % of the ablation area of Tiefmattengletscher (2017) as well as the lowest parts of Stockjigletscher and Schönbielgletscher. The surface elevation change maps were generated and investigated for each individual period as well as for the entire time span. A glacier's elevation change is typically stronger towards lower elevations, whereas the insulation effect of a continuous debris layer can reverse this gradient. To estimate whether such an effect can be discerned at Zmuttgletscher over time, we analysed the average elevation change of 50 m elevation bins (starting from 2125±25 m up to 2675±25 m) for seven periods since 1879.

Glacier-wide geodetic mass balance estimates were calculated between dates for which DTMs covered large parts of the glacier surface (1879, 1946, 1977, 1988, 2001, 2010, 2017). In areas of data voids or artefacts in the DTM, especially at higher elevations, no surface elevation change values could be calculated. In order to reduce the sensitivity to data voids and outliers, 100 m elevation bins (starting from 2100 m) were used for the elevation change calculation. In elevation bins without DTM coverage, the average of three extrapolation methods was used as the final mass balance value. The following three methods were used to fill the data voids: (1) a linear relationship between elevation and elevation change (e.g. Kohler et al., 2007; Kääb, 2008), which is based on the strong respective correlation up to R²=0.93; (2) the same relationship but additionally setting the elevation change values to 0 when they become positive in the highest elevations; (3) the mean elevation change value of the uppermost elevations that contain data and applied to the glacierised area above. The area-weighted average elevation change of all bins was summed to calculate the average elevation change of the total glacier surface. We assumed an average density of 850±60 kg m⁻³ (Sapiano et al., 1998; Huss, 2013) to convert the volume to mass changes.

3.7 Uncertainties

The total uncertainty of the surface elevation estimates and the mass balance consists in (i) the accuracy of the individual DTMs, (ii) the filling of the empty elevation zones, (iii) the glacier volume density, (iv) the debris volume changes, and (v) the DTM co-registration. Ice density changes are on average assumed to be negligible over a longer time span but we included an uncertainty of ±60 kg m⁻³ for the density to mass conversion (σ_dens); see Huss (2013). Debris volume changes were not accounted for due to a lack of such data, but they would have a negligible effect on the total uncertainty in volume change. The error from the filling of the elevation zones is taken into account as the standard deviation between the total mass balance after applying the three different interpolation methods (σ_fill). This uncertainty measure was combined with the uncertainties of the two DTMs framing the respective period (σ_DTM1 and σ_DTM2, standard deviation over stable terrain) to be used as the total uncertainty for the resulting mass balance by applying the law of error propagation (e.g. Hall et al., 2003; Zemp et al., 2013) and hence given by

$$\begin{array}{}\text{(1)}& {\mathit{\sigma }}_{\mathrm{total}}=\sqrt{{\mathit{\sigma }}_{\mathrm{DTM}\mathrm{1}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{DTM}\mathrm{2}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{fill}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{dens}}^{\mathrm{2}}}.\end{array}$$

Uncertainties range from ±0.12 to ±0.2 m w.e. yr⁻¹, whereas the uncertainty for the total period 1879–2017 is ±0.03 m yr⁻¹ (Table 2).

4 Results

4.1 Geodetic mass balance

Zmuttgletscher's long-term mass balance from 1879 to 2017 was $-\mathrm{0.31}\pm \mathrm{0.03}$ m w.e. yr⁻¹ (Table 2). The mass balance was least negative between 1879 and 1946 (Fig. 3a). Also during the climatically favourable period in the 1970s and 1980s Zmuttgletscher had less mass loss, even below the Swiss-wide mean (Fig. 3a). Between 1977 and 1988 the glacier-wide mass balance stayed negative ($-\mathrm{0.27}\pm \mathrm{0.18}$ m w.e. yr⁻¹), although certain areas on the tongue showed an increase in elevation (Fig. 10c, d, e). After 1988, the mass balance became more negative again, even though the lowest areas on the tongue still showed stable or even increasing surface elevation (Fig. 10f). The negative trend continuously intensified and after 2001 Zmuttgletscher's mass balance became more negative than the Swiss-wide mean (Table 2).

Table 2 Glacier-wide surface elevation changes and geodetic mass balance of Zmuttgletscher.

Figure 3 (a) Geodetic mass balance of Zmuttgletscher and Findelgletscher (data from Rastner et al., 2016) in comparison to all Swiss glaciers. Mass balance data from WGMS (2017). (b) Evolution of air temperature anomaly (T. anom.) at the climate stations Col du Grand St. Bernard (black lines, 2472 m, ∼40 km) and Sion (blue lines, 484 m, ∼35 km). Note that the Sion climate station was slightly repositioned in the 1960s (dashed vs. solid blue). The air temperature data were smoothed with a 5-year running mean.

4.2 Glacier area and length changes

From close to the end of the LIA (∼1850s) to 2017, Zmuttgletscher retreated by 1907±12 m (i.e. 12.1±0.09 m yr⁻¹; Fig. 4). The retreat rate peaked between 1961 and 1977 ($-\mathrm{21.7}\pm \mathrm{0.04}$ m yr⁻¹). Between 1977 and 2001, parts of the glacier terminus stagnated or even slightly advanced. Between 2001 and 2005 Zmuttgletscher accelerated its retreat, reaching $-\mathrm{19}\pm \mathrm{0.15}$ m yr⁻¹, and slightly decelerated to $-\mathrm{12.1}\pm \mathrm{0.06}$ m yr⁻¹ after 2005. The total glacier area decreased from 21.24±0.34 km² in 1859 to 15.74±0.02 km² in 2016 (Fig. 5). Similar to the length change, the area decreased only slightly between 1977 and 2001 ($-\mathrm{0.01}\pm \mathrm{0.006}$ km² yr⁻¹), after which the loss rate increased again to $-\mathrm{0.06}\pm \mathrm{0.004}$ km² yr⁻¹ (see also Figs. S1, S2, S3). Overall, Zmuttgletscher shows a relatively modest retreat while other debris-covered glaciers – Unteraargletscher and Glacier de Zinal – have retreated more rapidly (Fig. 4).

Figure 4 Cumulative length change of Zmuttgletscher with data points (vertical red lines) from this study, and several other Swiss glaciers (GLAMOS, 2018). For a list of total length and area changes see Fig. S1.

4.3 Debris cover

4.3.1 Debris cover evolution

At the end of the LIA 2.75±0.2 km² of Zmuttgletscher was debris-covered, which increased to 5.03±0.01 km² in 2013. During this time the total glacier area decreased from 21.24±0.4 to 15.82±0.3 km², resulting in an increase in percentage of debris cover from ∼12.9 % to 31.8±0.06 % of the glacier area from ∼1850 to 2013 (Fig. 5).

Figure 5 Evolution of total glacier area, debris-covered glacier area, and the share of the debris-covered glacier area.

Generally, the debris cover extent has expanded up-glacier in Tiefmattengletscher and Schönbielgletscher. The extension was pronounced along the surface of Schönbielgletscher and in the central part of the main glacier tongue (Fig. 6). In both areas the debris extent has expanded to above the icefall into the former accumulation areas and now is close to the foot of the contributing rock walls of Dent d'Hérens (Tiefmattengletscher) and Dent Blanche (Schönbielgletscher). The debris extent also strongly expanded at the glacier margins below Stockji, due to further input from moraines and the disconnection of a contributing tributary. By 2010, the icefall at Stockjigletscher had thinned sufficiently to disconnect and expose the rock wall beneath, from which a small rock fall detached between 2010 and 2013 (Fig. 6). This debris mound covered an area of approximately 170×300 m and is now being slowly transported down-glacier and spreading laterally.

Figure 6 Current extent of debris cover and its evolution since 1859 (∼ the end of the LIA) together with the glacier outlines for LIA and 2016. Even though the glacier has retreated in the debris-covered region, debris covers more and more of the total glacier area.

4.3.2 Debris cover characteristics

Field investigations revealed a homogeneous debris cover in some regions with stone sizes in the topmost layer mainly between 10 and 30 cm in diameter, and a much more heterogeneous debris cover in other regions, with pebble-sized stones to metre-scale boulders (Fig. S4). The typical base layer of the debris consists of fine-grained material (sand and even silt), and is overlain by a few centimetres of pebbles. The thickness of the debris layer along each transect varied from less than 5 cm to greater than 70 cm. The thickest areas were found on the elongated ridge, especially on the steep, southern slope (Fig. 7a, transect length 400–500 m). The average thickness along the upper transect was 16.3±1.3 cm.

Figure 7 (a) Debris thickness along the upper transect (Fig. 1). Uncertainty bars represent the standard deviation of three measurements per location for depths <20 cm. (b) Ablation of seven stakes over the time period 5 July–22 August 2017 (for locations see Fig. 1). For comparison the curves of Rakhiot Glacier, Barpu Glacier, Kaskawalsh Glacier, and Isfallsglaciären are displayed (adapted from Mattson et al., 1993).

Figure 8 Evolution of surface flow since 1961 for the different sections of the glacier tongue as labelled in the inset map. After an increase in the 1970s and 1980s a rapid slowdown occurred in all observed regions. SBG: Schönbielgletscher; TMG: Tiefmattengletscher.

Ablation measurements from seven locations and over a 7-week period in summer 2017 show an “Östrem-like” behaviour with respect to debris thickness (Fig. 7b; Östrem, 1959). These data indicate that melt is strongly dependent on debris thickness and much less on elevation. Compared to the reference stake at 2606 m on debris-free ice, field measurements in summer 2017 showed a reduction of ∼15 % and ∼50 % for debris thicknesses of 10 and 25 cm, respectively. These values are somewhat smaller than in the literature, where melt reductions of ∼40 %–60 % and 70 %–80 % for ∼10 and ∼25 cm, respectively, were found (Östrem, 1959; Mattson et al., 1993; Nicholson and Benn, 2006; Brock et al., 2007). Note that our debris-free reference stake is at the highest elevation of all stakes, which may partly explain our underestimation in melt reduction by debris. The backwasting rate of ice cliffs was 1.3 times higher than at the debris-free reference stake. The uncertainty of the ablation measurements is ∼3 cm, which results in an ablation rate uncertainty in the order of ±0.09 cm d⁻¹.

4.4 Surface flow velocities

Overall, there is a trend of decreasing velocities since the first measurement period of 1961–1977 until the most recent one (2010–2017), but with a clear phase of acceleration in the late 1970s and 1980s. Velocities in the lowest section of the glacier tongue (yellow) as well as the lower Schönbielgletscher have decreased from 10 to 20 m yr⁻¹ in the 1980s to almost zero (<3 m yr⁻¹ since 2005) and can be considered close to stagnant. The middle section (green) showed values of 30–40 m yr⁻¹ until the 1980s; afterwards velocities strongly decreased to ∼5–10 m yr⁻¹. Velocities in the highest section (blue) have decreased over the same period from 50 to 60 to ∼15 m yr⁻¹. In the periods 1977–1983 and 1983–1988 a velocity increase of close to 50 % was observed in all areas, with a slightly delayed signal down-glacier. It was followed by a strong velocity decrease in the 1990s that became weaker thereafter. In all time periods central glacier parts moved faster than the margins, and velocities strongly decreased down-glacier (Figs. 8, S5). Uncertainties are on the order of ±0.1 m yr⁻¹.

Figure 9 Influence of ice cliffs on total elevation change rate over the lowest area of the glacier tongue for eight time periods from 1961 to 2017. (a) Elevation change rate periods; (b) area that was considered “tongue” for elevation change analysis; (c) examples of ice cliffs of two different years (2005 in red; 2001 in blue), the corresponding borders of the 35 m buffer, and the overlapping area of this buffer (yellow, grey dotted area).

The displacements from radar interferometry from the 1.5 d period in summer 2017 yield similar results and confirm the quasi-stagnation of the lower tongue with line-of-sight flow speeds <2 m yr⁻¹ (Fig. S5). The more active part of the glacier with velocities >10 m yr⁻¹ extends up-glacier from just above the confluence of Tiefmattengletscher and Schönbielgletscher. During recent years, the flow in this mid-tongue section has strongly decreased and velocities have halved between 2005–2007 and 2016–2017 (Figs. 8 and S5).

4.5 Ice cliffs and related elevation change

Even though the glacier has become more debris-covered, we did not find any clear long-term trend in area and location of ice cliffs (Table 3). They are almost exclusively located in the lowest part of the glacier tongue (yellow area in Fig. 9). The total ice-cliff area decreased by over 40 % in the 1970s and early 1980s but increased again afterwards. The topographically corrected ice-cliff area amounts to less than 0.5 % of the total glacier area and less than 1.8 % of the debris-covered glacier area (Table 3). As only a few small ponds have been detected (<5 per time step) they have not been further analysed.

Figure 10 Evolution of elevation change over the glacier tongue since 1879. Ice-cliff areas are shown in red. Note how the “wave” of mass flux moves down-glacier from (c)–(f). High values in the glacier forefield (e.g. 1983–1988) are linked to artefacts from water surfaces and the activities of a gravel mining company. See Figs. S7 to S16 in the Supplement for enlarged panels.

Table 3 Evolution of ice-cliff area from 1946 to 2013; dc.: debris-covered.

Over the entire period, the average elevation change of ice-cliff areas (cliffs + buffer overlap) was 1.2 times higher than over the rest of the tongue (Fig. 9). During the 1970s and 1980s the thinning rates were almost the same and close to zero whereas since 1988 the ratio of average elevation change ranged between 1.5 and 1.7. When considering average thinning rates over the entire glacier tongue from 1995 to 2017, the inclusion of ice-cliff areas enhances the surface lowering only by less than 5 % compared to average thinning of the glacier tongue without ice cliffs. Thus the presence of ice cliffs does not seem to affect the mass loss substantially.

Figure 11 Average annual elevation change as a function of surface elevation (in 50 m elevation bins) for seven time periods since 1879.

4.6 Tongue-wide surface elevation change and patterns

Since the end of the LIA, the surface elevation of the Zmuttgletscher tongue below 2750 m a.s.l. has almost continuously decreased (Fig. 10). Between 1879 and 1946, the average elevation change was negative at $-\mathrm{0.63}\pm \mathrm{0.13}$ m yr⁻¹, and substantially decreased to $-\mathrm{1.75}\pm \mathrm{0.34}$ m yr⁻¹ between 1946 and 1961. In the 1970s and 1980s the tongue partially thickened, which is most pronounced at the tip of the Stockjigletscher tongue ($+\mathrm{1.2}\pm \mathrm{0.29}$ m yr⁻¹ from 1977 to 1983) and in the main trunk of the tongue, while the terminus area was still substantially thinning. Between 1977 and 1983 the tongue's average surface elevation barely changed, at a rate of $-\mathrm{0.11}\pm \mathrm{0.29}$ m yr⁻¹. Between 1983 and 1995 the upper part of the tongue already started to thin again whereas in the lowest 1 km a mixed signal of thickening and thinning occurred (Fig. 10f), consistent with a slight advance at the centre of the terminus after 1995. After 1988 the tongue began to thin substantially at a rate of −1 m yr⁻¹, becoming even more negative with $<-\mathrm{2}$ m yr⁻¹ after 2005 (see Fig. S6 for exact values of all periods). Since 1879, the average thinning over the ablation area has been 104.7 ±0.05 m (∼0.76 m yr⁻¹), with maximum values of >190 m in the area of the terminus position in 2017, i.e. the former central part of the tongue. Before 1946, surface lowering was greatest near the terminus (Fig. 10a–d), whereas later it became more heterogeneous throughout the tongue, particularly since 2001. Locally, surface lowering was most pronounced in ice-cliff areas and in the debris-free upper regions (panels f, h, i in Fig. 10).

The relation between surface elevation change and surface elevation has strongly shifted over time (Fig. 11). During the periods from 1879 to 1988 thinning was stronger in lower elevations, with an almost stable surface above 2500 m in the period of 1977–1988. After 1988 the thinning has in general become stronger but also more homogeneous along the glacier tongue.

5 Discussion and interpretation

5.1 Method suitability and shortcomings

The DTMs used for calculating mass balance and surface elevation changes have uncertainties of ∼1–2 m (i.e. ±0.2–0.4 m yr⁻¹; Table 1), which are derived as the standard deviation over stable terrain and therefore represent the upper bound of the uncertainty (Magnússon et al., 2016). Except for the periods 1879–1946 and 1977–1988, these uncertainties are significantly lower than the glacier changes and thus have little impact on the final decadal elevation change rates. High elevation change values outside of the glacier surface are found in very steep terrain and in the immediate forefield, and in the latter case stem from artefacts on water surfaces (e.g. Fig. 10e), or the activities of a gravel mining company (e.g. Fig. 10g, h, i), and do not affect the results on glacier elevation change.

The time series of glacier surface flow velocities along the lowest ∼4 km is over 50 years long (1961–2017) and includes periods of glacier acceleration and deceleration. Even though it covers less than half of the total study period, it is one of the longest time series of a debris-covered glacier, comparable to the velocity series from Glacier de Tsarmine from 1967 to 2005 (Capt et al., 2016). Despite difficulties in tracking boulders in more active areas on Zmuttgletscher due to large displacements (long periods between image dates) and crevasses (e.g. massive increase in crevasses in 1980s), several representative data points could be extracted per glacier subsection and time period, which provide a clear picture of the dynamic evolution of the entire glacier tongue since 1961. Unfortunately, the time period 1946–1961 yielded too few data points (due to image quality and snow cover) and had to be excluded from the time series. In the lowest areas the velocity field derived from the 2016–2017 orthophotos (Fig. S14) shows slightly higher values than the interferometer (e.g. 4±0.25 m vs. 1±0.57 m yr⁻¹) and a reduced extent in low flow, but confirms the stagnation of the lower tongue. This difference can be attributed to the short duration and time of observation in late summer, when flow velocities are likely below the annual average due to the establishment of an effective basal drainage system (e.g. Nienow et al., 1998).

The uncertainty of the mapped debris extent is equal to that of the glacier and ice-cliff mapping at $\pm \mathrm{1}/\mathrm{2}$ pixel. For 1961, the debris cover extent might be slightly underestimated as the orthophoto does not fully cover the ablation area and a thin snow-cover was present in higher areas. This could partly explain the slight reduction in debris extent from 1946 (3.98±0.01 km²) to 1961 (3.79±0.01 km²).

The automatic detection of ice cliffs using topographic and geometric variables as applied by others (e.g. Brun et al., 2016; Kraaijenbrink et al., 2016) showed no satisfying results. Many steep areas (>40^∘) were incorrectly identified as debris-covered, and slope values extracted over ice cliffs were often below those typically observed in situ elsewhere (e.g. Reid and Brock, 2014). Our object-based mapping approach is efficient to map larger numbers of cliffs for several orthophotos. However, the polygons are based on a somewhat subjective choice of separation parameters and the manual correction of the polygons is time-consuming. When considering image quality (texture, contrast, snow cover, cast shadow) the results include both omission and commission errors, which also lie in the range of $\pm \mathrm{1}/\mathrm{2}$ pixel along the polygon margins.

5.2 Chronological evolution and process interactions

On the basis of their dynamic behaviour we distinguish four distinct phases over the entire observation period since 1879. In the first phase, between 1879 and 1977, as a response to the atmospheric warming in the mid-19th century, the mass balance of Zmuttgletscher turned negative (Fig. 3), first only slightly (1879–1946: $-\mathrm{0.09}\pm \mathrm{0.12}$ m. w.e. yr⁻¹), then more strongly (1946–1977: $-\mathrm{0.67}\pm \mathrm{0.16}$ m. w.e. yr⁻¹), leading to glacier thinning and retreat. At the same time, the debris-covered area increased from 2.75±0.2 km² in 1859 to 3.67±0.01 km² in 1961. Rising air temperatures are probably the main reason for the increasing debris extent as they lead to a rise of the rain–snow transition altitude as well as higher ice ablation. Therefore, debris emergence accelerates and moves further up-glacier (Stokes et al., 2007). Decreasing ice flow further supports the development of a continuous debris cover of substantial thickness as the emerging debris is evacuated more slowly, and the debris has more time to melt out and thicken (Kirkbride and Deline, 2013). In the mid-20th century, a continuous debris cover existed on the lower part of the tongue, likely already responsible for a flattening of the mass balance gradient. Ice cliffs and supraglacial meltwater channels already existed in the lowest, debris-covered part of Zmuttgletscher in 1930 (Fig. S17) and thereafter increased in area until 1977. These surface features are partly responsible for increased thinning in these areas close to the terminus compared to further up-glacier.

In a second phase between 1977 and 1988, atmospheric cooling and increased precipitation brought a shift towards less negative mass balance ($-\mathrm{0.27}\pm \mathrm{0.18}$ m w.e. yr⁻¹, Fig. 3). This is reflected in a partial thickening of the tongue (Fig. 10d, e) as well as in increased flow velocities (Figs. 8, 12b). At the same time, the upper boundary of debris cover in these areas migrated slightly down-glacier (Fig. 12a), but because the debris cover continued to expand in other areas, the total debris-covered area was roughly stable during this period. The upper boundary of the zone with most ice cliffs was also moved down-glacier during this period and the total area of ice cliffs dropped by a factor of 2 (Table 3). The adjustment of the glacier length and area was less direct, but the retreat and area loss gradually slowed down.

In the third phase between 1988 and 2001 the mass balance became more negative again – along with an abrupt rise in air temperature (Fig. 3a, b) – and flow velocity reductions were pervasive across the glacier tongue, first higher up, then also in the lowest regions (Figs. 8, 12b). The debris-covered area substantially increased – especially during the 1990s – while the glacier area was stable and the central part of the terminus even advanced for some metres, leaving behind a small moraine (Fig. S18). The delayed response to the preceding period with less negative mass balance is, however, still detectable in the slight elevation gain right at the terminus in the period 1995–2001, whereas it was already highly negative on the rest of the tongue (Fig. 10f).

In the most recent period (2001 to 2017), the glacier has continued to develop a more negative mass balance (approx. −0.8 to −1.0 m w.e. yr⁻¹). Thinning has been most pronounced in debris-free areas – especially on the tongue of Stockjigletscher – and in the ice-cliff zone, and rates remained high (>1 m yr⁻¹) across the tongue. Towards the terminus these thinning rates are, however, not comparable to values typically observed on debris-free glacier tongues at similar elevations (e.g Findelgletscher: >5 m yr⁻¹; Sold et al., 2016), which might explain the relatively small length and area changes given the high overall mass loss (Figs. 3, 4). Flow velocities have continuously decreased after 1988, and since ∼2001 the lowest 1.5 km of the tongue has become almost stagnant (mostly below 3 m yr⁻¹).

Figure 12 (a) Central elevation profiles over the glacier tongue from the 19th century until 2017; the red line (2017) marks the former glacier bed from 0 to ∼2000 m along the profile; the rugged lower part of the lines is shaped by ice cliffs and flow channels. (b) Surface flow velocities from the 1960s until 2017 along the same profile; velocities were highest during the positive mass balance period in the 1970s and 1980s; the lowest part of the tongue is almost stagnant (2017). Vertical brown and blue bars show the upper boundaries of the debris-covered area and the ice-cliff zone, respectively; both boundaries have interrupted their up-glacier trend between 1977 and 1988.

5.3 Unprecedented debris cover increase

Zmuttgletscher's debris extent has increased by ∼19 % points since the end of the LIA without evidence of particularly large rock falls. Here the headwalls surrounding the accumulation area are the most relevant input sources, but lateral moraines also have a minor influence in marginal areas (van Woerkom et al., 2019). Increases in debris cover extent (for periods from several years to decades) have been shown for other Alpine glaciers (Deline, 2005; Kellerer-Pirklbauer et al., 2008), the Caucasus (Stokes et al., 2007), and the Himalayas (Bolch et al., 2008a; Bhambri et al., 2011), but have so far not been quantified for a century-long period. These studies found debris cover extent increases of 2 %–10 % points for glaciers with absolute debris cover extents of 2 %–42 %, which is well below the rates we found for Zmuttgletscher, even on a decadal scale. At Zmuttgletscher we could also identify strong temporal variations in the rate of debris cover extent change (e.g. stable extent in 1970s and 1980s and high rates of change of >4 % points per decade in the 1990s and early 2000s), which coincide with distinct changes in climate and related dynamic responses. This strong link to climate may also help to explain the strong debris cover extent increase on Zmuttgletscher, which we suggest is the result of (i) a combined effect of elevated temperatures and decreasing velocities, (ii) the large debris-free ablation area that still existed a few decades ago, and likely also the (iii) the comparatively small glacier size, i.e. the shorter response time.

5.4 The influence of debris cover on geometry

To understand the potential long-term effects of a debris cover on glacier evolution, regional comparisons provide a useful context. Nevertheless, one has to be careful in comparing length and area changes of glaciers for a given climate history because of differences in geometry and hence response times (Jóhannesson et al., 1989). The retreat of Zmuttgletscher is relatively modest compared to that of other large Swiss glaciers (Fig. 4) and it has shown little terminus fluctuation, even during the climatically favourable period in the 1970s and 1980s. Other glaciers with similarly subdued fluctuations are Unteraargletscher and Glacier de Zinal, which are also debris covered in their lower reaches. Unlike smaller and debris-free glaciers, neither Unteraargletscher nor Glacier de Zinal advanced in the 1980s and 1990s (Fig. 4). Aletschgletscher (>80 km²) is too large and thick to react to short periods of positive mass balance and shows a smooth and accelerating long-term retreat trend since the end of the LIA (i.e. −24 m yr⁻¹; GLAMOS, 2018). Conversely, Findelgletscher, a nearby debris-free glacier of similar size, thickness, and elevation range as Zmuttgletscher, showed much stronger fluctuations in length and mass balance. Findelgletscher experienced periods of balanced and even positive mass balances, accompanied by a sharp terminus advance in the 1980s (+41 m yr⁻¹), followed by a strong retreat since 1985 (−44 m yr⁻¹). During these recent periods of strongly negative mass balance, the retreat rate at Zmuttgletscher ($-\mathrm{7.2}\pm \mathrm{0.01}$ m yr⁻¹) was lower than that of all other glaciers in Fig. 4.

Area changes from the beginning of the 1970s until ∼2012 reveal that Zmuttgletscher has lost 7 % of its area (−0.17 % yr⁻¹), which is comparatively low (Table S1). In comparison, Findelgletscher has lost more than 11 % of its area (−0.32 % yr⁻¹) during a similar period. Thus, in terms of length and area changes, Zmuttgletscher shows a distinctly subdued and delayed response to climatic perturbations compared to debris-free glaciers in the region.

The relatively uniform and reduced thinning over much of Zmuttgletscher's tongue can be attributed to the debris cover and is observed at several debris-covered glaciers in Switzerland. This suppression of thinning by debris cover on the lower ablation area is best illustrated when plotting the elevation change for one period (1980–2010) for 11 larger Swiss glaciers along 10 vertical elevation bins of equal size (Fig. 13). For all debris-covered glaciers, including Zmuttgletscher, the thinning becomes almost independent of elevation for the lower elevation bins, whereas for debris-free glaciers, it continues to increase towards the terminus. The non-linear relationship between thinning and elevation causes the slower adaptation of debris-covered glacier length and area and hence more extended tongues. A similar flattening or even reversal of the elevation change gradient towards debris-covered glacier termini has also been observed in other regions, e.g. in the central Himalayas (Inoue, 1977; Bolch et al., 2008a; Benn et al., 2012; Ye et al., 2015; Ragettli et al., 2016) and the Tien Shan (Pieczonka and Bolch, 2015).

Figure 13 Elevation change between 1980 and 2010 against normalised elevation for a selection of Swiss glaciers. Solid lines: mostly debris-free glaciers; dashed lines: debris-covered glaciers. Data from Fischer et al. (2015). The lowest section of Aletschgletscher (yellow line) is actually also slightly debris covered.

5.5 Climate-driven glacier mass balance and flow dynamics

The evolution of glacier-wide mass balance of Zmuttgletscher since 1879 is comparable to centennial trends in other debris-free glaciers in the Swiss Alps, which are closely related to variations in climate (Fig. 3; Schmidli, 2000; Bauder et al., 2007; Huss et al., 2010; Hirschi et al., 2013; World Glacier Monitoring Service, 2017). For Zmuttgletscher, the strong connection of mass balance to climate is exemplified by the almost balanced conditions during colder phases before 1920 and in the 1970s and 1980s as well as by the strongly negative mass balance in the 1940s and 1950s and in recent decades (Fig. 3; Bauder et al., 2007; Escher-Vetter et al., 2009). This clearly demonstrates that the glacier-wide mass balance of Zmuttgletscher has foremost been governed by climatic changes and has not been strongly influenced by debris cover and changes thereof. This is in contrast to the case of Glacier de Miage in Italy, for which Thomson et al. (2000) related the positive mass balance on the tongue over the 20th century partly to debris cover. The debris cover of Glacier de Miage is, however, thicker than at Zmuttgletscher.

The relatively thin debris at Zmuttgletscher may on the one hand limit the decoupling from climatic influences, and on the other hand still be sufficient to reduce glacier thinning and terminus retreat. This results in an extended and stagnating glacier tongue that increases the area of ablation (relative to a debris-free glacier tongue), which in turn enhances the sensitivity of the glacier-wide mass balance to climatic changes.

The observed temporal variations in ice flow and thickness on the main glacier tongue of Zmuttgletscher also closely reflect variations in climate as exemplified by the acceleration and thickening in the period of more positive mass balance in the 1970s and 1980s. Positive mass balances impacting flow velocities were also observed at debris-covered Glacier de Tsarmine (Capt et al., 2016) and a number of other debris-free glaciers in Switzerland (Glacier de Giétro, Glacier de Corbassière, Mattmark; Bauder, 2017), Austria (Hintereisferner, Kesselwandferner; Stocker-Waldhuber et al., 2018), and France (Glacier d' Argentière; Vincent et al., 2009). Also at Glacier de Miage, a kinematic wave migrating down-glacier was observed between the 1960s and 1980s that led to a small terminus advance in the late 1980s (Thomson et al., 2000). Also the long-term (1960s until 2017) reduction in flow velocities on the tongue of Zmuttgletscher is in line with observations from all of the above-mentioned glaciers.

All these examples show a direct reaction of flow velocities and hence ice flux to climatic changes. This link is consistent with the principle of mass conservation, regardless of debris coverage, but breaks down when the glacier tongue starts to stagnate as a result of the debris cover slowing down the retreat. Such a dynamical stagnation and decoupling has in the last decade also been observed at the tongue of Zmuttgletscher and is characteristic of strongly debris-covered glaciers (e.g. Bolch et al., 2012).

On large debris-covered glaciers, lower flow velocities are often the result of a decreased driving stress due to flat tongues that result from a sustained reversal of the mass balance gradient (Bolch et al., 2008a; Anderson and Anderson, 2016). At Zmuttgletscher, the elevation change gradient has indeed decreased to almost zero during recent decades. Nevertheless, the glacier surface slope has slightly increased over time; therefore we attribute the decrease in dynamic activity to the reduction in ice thickness and, hence, mass flux and climate. This conclusion is coherent with findings by Dehecq et al. (2019) from glaciers in south and central Asia.

Figure 14 Examples of superficial and englacial ablation and their consequences. (a) Meander of a large supraglacial water channel with confining, under-cut ice cliff. (b) Almost circular surface lowering as a consequence of en- and subglacial ablation. (c) Subglacial water flow channel. (d) Supraglacial pond and ambient cliffs.

5.6 The role of debris-related surface features

Whilst the presence of a debris mantle may not substantially influence the overall mass loss rate of Zmuttgletscher, we find that the main reason for the observed heterogeneous thinning pattern is the increasingly extensive and also likely thicker debris cover that is heterogeneously distributed over the tongue. The heterogeneity is further increased by the presence of ice cliffs. Field visits and patterns of surface lowering suggest a close association between ice-cliff formation and the presence of supraglacial meltwater channels. Surface meltwater often runs off superficially over substantial distances and water flow channels are abundant all over the glacier tongue outside of crevassed areas. Inside and alongside these channels bare ice becomes exposed when water washes away the debris or laterally cuts into the ice (Fig. 14a) and the debris slides off the oversteepened channel walls. The location of supraglacial meltwater channels on Zmuttgletscher seems closely related to areas of compressional flow – in flat and stagnating areas – and has been most pronounced in the lowest ∼1.5 km, where most ice cliffs also occur. This is consistent with the cessation in up-glacier migration of the ice-cliff boundary in 2005, just below a step in the surface with extensional flow (Fig. 12a at a profile distance ∼3500 m). Ice cliffs also develop when en- and subglacial voids collapse and the shear planes of the fractured ice become exposed as also observed at Zmuttgletscher (Fig. 14b). We suggest that englacial ablation is effectively creating and enlarging such voids (Benn et al., 2012; Immerzeel et al., 2014; Fischer, 2011; Fig. 14c). Supraglacial ponds are known to have a high potential to increase ablation (Sakai et al., 2000; Röhl, 2008; Miles et al., 2016) and are also potential origins of ice cliffs (Fig. 14d) but these features are not prevalent on Zmuttgletscher.

Large ice cliffs also exist at the terminus and are responsible for exacerbated retreat compared to where the terminus is gently sloping and debris mantled. The situation of a terminal cliff at the glacier mouth combined with terminus retreat is also found at Gangotri glacier for example (Bhambri et al., 2012; Bhattacharya et al., 2016), whereas some glaciers with a stable terminus (e.g. Khumbu, Miage) do not show such terminal cliffs (Bolch et al., 2008a; Diolaiuti et al., 2009). Remarkably, depressions and irregular surface topography near the terminus were already indicated in the Siegfried map from 1879. Certain ice cliffs on Zmuttgletscher have reached >25 m in height and have persisted for over 2 decades. The consequences of ice cliffs at Zmuttgletscher are (i) a more chaotic pattern of surface elevation changes due strong local ablation, (ii) debris redistribution through cliff backwasting and fluvial transport, (iii) and stronger surface elevation changes at the lower part of the tongue (especially downstream of topographic steps). As a result, the elevation change is still stronger towards the terminus than without ice cliffs.

The strong expansion of debris cover might suggest a high and increasing importance of ice cliffs. However, we found that ice cliffs are only responsible for approximately 5 % of glacier-wide volume loss because (i) their total area is small (Table 3; also compared to other glaciers: e.g. Sakai et al., 2000; Juen et al., 2014; Ragettli et al., 2016) and (ii) the debris is relatively thin on Zmuttgletscher.

5.7 Dynamic interaction between flow velocities, debris cover, and ice cliffs

Even though high flow velocities may be expected to subdue the emergence of ice cliffs, we could not find a clear velocity threshold linked to their occurrence, e.g. to designate areas of potential ice-cliff formation. However, the upper boundary of the ice-cliff zone was moved downstream by increased flow velocities in the 1970s and 1980s (Fig. 12). Further, during this period the total ice-cliff area was clearly reduced compared to both periods with lower flow speeds before and after (Table 3). This suggests a clear link between dynamic state and the occurrence of ice cliffs and would imply expanding ice-cliff areas on stagnating tongues, consistent with the general interpretation of observations (e.g. Pellicciotti et al., 2015).

Increased flow velocities also led to a local down-glacier movement of the upper boundary of the debris extent (Fig. 12). Similarly, Deline (2005) documented a decrease in debris-covered area for Glacier de Miage between ∼1890 and ∼1920, caused by an increase in surface elevation and flow velocity (Thomson et al., 2000), which confirms that changes in glacier dynamics can directly impact supraglacial debris extent.

6 Conclusions

This study presents a >150-year-long reconstruction of the geometry, surface topography, and debris cover of Zmuttgletscher, as well as observations of surface flow velocities over 5 decades. We observed the strongest debris cover extent increase on record, from 12.9±0.9 % of the total glacier surface in 1859 to 31.8±0.06 % in 2013. The heterogeneous debris thickness distribution reduces ablation from 15 % to 80 % (compared to bare ice) over large parts of the debris-covered area, causing persistent debris-related thinning patterns. Even though this debris cover is relatively thin (mostly ∼15 cm), it has efficiently attenuated and delayed the response to climatic changes in terms of glacier length and area adjustment compared to debris-free glaciers in the Alps, during periods of both mass gain and mass loss. Hence, elevation change has become (almost) independent from elevation. Despite the expansive debris cover, our observations demonstrate that the evolution of Zmuttgletscher has been mainly driven by climatic changes, on both decadal and centennial timescales. Glacier-wide mass balances are throughout the last century comparable to other debris-covered but also debris-free Swiss glaciers, corroborating a direct link to climate. The prevalence of localised ablative processes such as ice-cliff expansion caused surface elevation change in otherwise heavily debris-covered areas, which leads to increased ice loss in the lowest glacier parts, thereby preventing a reversal of the elevation change gradient.

In general, flow velocities have strongly decreased in the last 2 decades similar to other glaciers in the Alps, and the lower 1.5 km of the tongue is almost stagnant (2017), even though the glacier tongue has slightly steepened over time. The low flow velocities are due to thinning and reduced ice flux, and to a minor degree influenced by debris cover. Higher flow velocities between 1977 and 1988 were triggered by a more positive surface mass balance and the related increased mass flux from the accumulation area, which also caused local glacier thickening. This increase led to a slight down-glacier migration of the debris cover boundary, followed by a strong debris cover increase when velocities dropped in the 1990s and rising air temperature led to increasingly negative mass balance. Higher flow velocities also moved the upper boundary of the ice-cliff zone down-glacier, temporarily reduced the ice-cliff area, and eventually led to a slight advance in the 1990s.

These findings suggest a clear and direct influence of flow dynamics on debris extent and ice-cliff formation. The above-described processes and feedbacks are likely valid and relevant for other debris-covered glaciers in the Alps and elsewhere at potentially different rates and magnitudes. In the context of global warming, it is therefore crucial to include these findings in models for glacier projections.