Darkening Swiss glacier ice?

The albedo feedback is an important driver of glacier melt over bare-ice surfaces. Light-absorbing impurities strongly enhance glacier melt rates but their abundance, composition and variations in space and time are subject to considerable uncertainties and on-going scientific debates. In this study, we assess the temporal evolution of shortwave broadband albedo derived from 19 end-of summer Landsat scenes for the bare-ice areas of 39 large glaciers in the western and southern Swiss Alps. Trends in bare-ice albedo crucially depend on the spatial scale considered. No significant negative temporal trend in bare5 ice albedo was found on a regional to glacier-wide scale. However, at higher spatial scales, certain areas of bare-ice including the lowermost elevations and margins of the ablation zones revealed significant darkening over the study period 1999 to 2016. A total glacier area of 16 km2 (equivalent to about 12% of the average end-of-summer bare-ice area in the study area) exhibited albedo trends significant at the 95% confidence level or higher. Most of this area was affected by a negative albedo trend of about −0.05 per decade. Generally, bare-ice albedo exhibits a strong interannual variability, caused by a complex interplay 10 of meteorological conditions prior to the acquisition of the data, local glacier characteristics and the date of the investigated satellite imagery. Although, a darkening of glacier ice was found to be present over only a limited region, we emphasise that due to the recent and projected growth of bare-ice areas and prolongation of the ablation season in the region, the albedo feedback will considerably enhance the rate of glacier mass loss in the Swiss Alps in the near future.


Introduction
Glaciers are known to be excellent indicators of climate change (IPCC, 2013).Increasing air temperatures and changing precipitation patterns provoke snowlines to rise to at higher altitudes and thus a spatially greater exposure of bare-ice surfaces.
In connection with a general prolongation of the ablation season, the increased climatic forcing causes an amplification of glacier melt (Kuhn, 1993;Oerlemans, 2001).However, these changing glacier characteristics trigger feedback mechanisms, in particular the positive albedo feedback which enhances bare-ice melting (Tedesco et al., 2011;Box et al., 2012).Hence, the strongly negative mass balances of many glaciers are not solely a direct signal of atmospheric warming but result from a complex interplay of changes in climate forcing and related surface-atmosphere feedback mechanisms.
Currently, there is an on-going debate about the occurrence and rate of glacier and ice sheet darkening worldwide.While studies like those of Oerlemans et al. (2009), Wang et al. (2014), Takeuchi (2001) or Mernild et al. (2015) observed a darkening for one or several glaciers, respectively in the European Alps, the Chinese and Nepalese Himalaya or in Greenland's peripheral 1 The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-18Manuscript under review for journal The Cryosphere Discussion started: 29 January 2018 c Author(s) 2018.CC BY 4.0 License.glacierised areas over varying time-scales, evidence of darkening from sectors of the Greenland Ice Sheet is less pronounced, leading to controversial discussions (cf.Box et al., 2012;Alexander et al., 2014;Dumont et al., 2014;Polashenski et al., 2015;Tedesco et al., 2016).Moreover, the emergence of legacy contaminants at lower elevations on Alpine glaciers (Bogdal et al., 2009;Pavlova et al., 2014;Steinlin et al., 2014Steinlin et al., , 2016) ) or outcropping ice in the ablation zone of the Greenland Ice Sheet that contain high dust concentrations potentially associated with paleo-climatic conditions (Wientjes and Oerlemans, 2010;Wientjes et al., 2012;Goelles and Bøggild, 2015) and the recognition of the potential role of biological impurities (Cook et al., 2017;Stibal et al., 2017;Tedstone et al., 2017) may emphasize and amplify the impact of light-absorbing impurities on ice melt.
To date, most long-term studies either used point data from automatic weather stations located in the ablation area of a glacier (e.g.Oerlemans et al., 2009), coarsely spaced satellite data from the Moderate Resolution Imaging Spectroradiometer (MODIS) (e.g.Stroeve et al., 2013;Mernild et al., 2015) or other remote sensing datasets (e.g.Wang et al., 2014) to infer trends in ice albedo.Mostly, studies validated the satellite-derived albedo values with in-situ data measured at one to several locations on the ground, which is not necessarily always ideal (Ryan et al., 2017).However, for the limited size of Alpine glaciers and the complex surrounding topography, the spatial resolution of MODIS data (500 m) is not suitable, and no appropriate highresolution albedo product is readily available.Thus, studies focusing on alpine glaciers often base their analysis on higher resolution datasets to obtain information about glacier surface albedo and related changes and processes (Dumont et al., 2011;Fugazza et al., 2016;Di Mauro et al., 2017).Furthermore, when debating the darkening of glaciers, a clear distinction between glacier-wide versus point-based investigations is necessary.Moreover, a separation between bare-ice and snow-covered albedo changes is required to correctly distinguish between differing processes and dependencies.
In this study, we use Landsat data to obtain spatially distributed bare-ice albedo for 39 glaciers with a total area of 480 km 2 (corresponding to about a quarter of the present glaciation of the European Alps) located in the western and southern Switzerland over the 17-year period 1999 to 2016.We focus on the bare-ice areas defined as neither snow-covered nor debris-covered glacier surfaces, as the latter affect glacier mass balance by different processes.We examine trends and their significance to better quantify and investigate a possible darkening of glacier ice in the western and southern Swiss Alps from the point to the regional scale.

Study sites and data
Our study focuses on 39 glaciers located in the western and southern Swiss Alps (Figure 1).All of them are characterised by a surface area >5 km 2 and thus offer a large-enough spatial extent to study the evolution of bare-ice surfaces and related albedo changes.The investigated glaciers vary considerably in size, ranging from about 5 km 2 (Giétro, Schwarzberg) to almost 80 km 2 (Aletsch), and span an elevation range from about 1850 m above sea level (m a.s.l.) to over 4500 m a.s.l.(Table 1).According to the most recent Swiss Glacier Inventory, the 39 glaciers covered a total area of 483 km 2 in 2010 (Fischer et al., 2014).We used a Sentinel-2 scene (20 m spatial resolution) acquired on the 23 rd of August 2016 to manually adjust the glacier outlines and to obtain up-to-date glacier extents totalling 451 km 2 of high accuracy (Paul et al., 2016).For our analysis, we excluded heavily debris-covered parts, except for medial moraines, as we focus on the albedo of bare ice only.Hence, the area difference of 32 km 2 between 2010 and 2016 does not solely stem from glacier retreat, but is also due to our exclusion of debris-covered glacier tongues for this study (e.g.Zmutt, Unteraar, Zinal, Oberaletsch).The obtained glacier outlines for the year 2016 are used in all consecutive analysis, thus the glacier outline is kept constant over the study period and does not evolve with time.2).To obtain maximum information about the bare-ice area of the glaciers, only scenes acquired at the end-of summer (from months August or September) were chosen.Unfortunately, no good scenes are available for the years 2001, 2007 and 2010.The restriction to choose end-of summer scenes only hampers the investigation of seasonal changes, but favours an intercomparison over multiple years.On average, a scene comprised 132 km 2 of bare ice, of which usually less than 5% was cloud-covered (Table 2).

Bern
Our surface type retrieval approach based on the obtained broadband shortwave albedo (see Section 3.2) requires a digital elevation model.We used the DHM25 with an original spatial resolution of 25 m provided by Swisstopo [Swisstopo, 2005] and resampled it to 30 m spatial resolution to match the broadband shortwave albedo datasets derived from Landsat.
To contextualise our results, the lithology surrounding the individual glaciers based on the lithological-petrographic map of Switzerland (GK500) provided by Swiss Geotechnical Commission (SGTK) was used; the map is at 1:500000 scale, and shows the subsurface strata subdivided into 25 groups according to their formation, their mineralogical composition, their particle size and their crystallinity.Moreover, we also used observed annual mass balances for individual glaciers acquired in the frame of the Swiss glacier monitoring network (GLAMOS) to study the link between albedo and glacier mass loss.

Pre-processing
All reflectance data was downloaded through earthexplorer.usgs.gov(Table 2).The final selection of all scenes is based on a visual check.To detect and delineate clouds obscuring the glacier surfaces, we used a semi-automatic classification approach based on the Spectral Angle Mapper (SAM, Kruse et al. (1993)) algorithm implemented in ENVI.For each sensor (TM, ETM, OLI) a spectral library of cloud signatures was manually compiled, which served as reference library for the respective sensor.
Hence, for each scene we obtained a cloud mask that was used to exclude cloud-affected pixels from all consecutive analyses.
Except for the scene taken on the 09 th of September 2013, when about 20% of the bare ice was cloud-covered (north-east part of the study area), the cloud coverage was generally smaller than 5% over the bare-ice area (Table 2).

Albedo retrieval
We applied the narrow-to-broadband conversion by Liang (2000) to obtain shortwave broadband albedo α short from the surface reflectance data.The conversion is based on five of the seven individual bands, and is formulated as follows: where b n represents the spectral band number of Landsat TM/ETM.For Landsat OLI, the band numbers were adjusted accordingly.This conversion was developed based on a large empirical data set and the band configurations of Landsat TM/ETM.
As shown by Naegeli et al. (2017) this albedo retrieval approach can be applied also to the most recent mission Landsat 8 and is suitable for mountain glaciers.It provides albedo products that are of very high accuracy and deviate by less than < 0.001 on average from a more sophisticated albedo retrieval approach.

Surface type evaluation
The delineation of bare-ice area versus snow-covered surfaces is based on a multi-step classification scheme of the surface albedo values (Figure 2).The classification is thus based on a physical parameter specific for both snow and ice.In a first step, two threshold values for certainly snow (α > 0.55) and certainly ice (α < 0.25) are defined (primary surface type evaluation, Figure 2) based on recommendations in the literature (Cuffey and Paterson, 2010).This results in a critical albedo range (0.25 > α < 0.55), where an unambiguous assignment of the surface type, i.e. snow or ice, is not possible without considering other parameters.We, therefore, take advantage of a digital elevation model available for all glaciers to evaluate the average albedo in elevation bands of 20 m within this critical albedo range.The transition between ice and snow is typically characterized by a distinct change in albedo (e.g.Hall et al., 1987;Winther, 1993;Zeng et al., 1983).We thus derive an estimate of the mean snowline altitude (SLA) for each glacier and scene based on the greatest slope of the albedo-elevation profile.The albedo for this altitude is considered to be the site-and scene-specific albedo threshold discriminating snow and ice and is henceforth termed α crit .
In a second step, we use the SLA and α crit as reference to evaluate the surface type within the range of critical albedo values, where there is ambiguity between snow and ice (secondary surface type evaluation, Figure 2).Finally, all grid cells are evaluated regarding their relative position compared to the SLA (probability test to eliminate extreme outliers, Figure 2).
Grid cells located clearly above the SLA are more likely to be snow than ice, and vice versa.An increasing positive/negative  vertical distance from the SLA thus results in penalties for the likelihood of the cell within the critical albedo range of being either snow or ice.As an example, a grid cell near the glacier terminus with an albedo of 0.42, i.e. a rather high albedo for Alpine glacier ice, will be classified as ice.An albedo of e.g.0.35 observed for the highest regions of the glacier, in contrast, will be classified as snow, as the low albedo is more likely to be explained by an erroneous albedo determination (e.g.shadows)

Probability test to eliminate extreme outliers
than by actually snow-free conditions.In summary, our procedure to distinguish between snow and bare-ice surfaces relies on 5 remotely-determined surface albedo and merges this information with surface elevation with a probability-based approach to detect outliers and to automatically adapt the classification to the site-and scene-specific conditions.

Trend analysis
Over the study period 1999 to 2016, at least one good Landsat end-of-summer snapshots was available for 15 years (cf.Table 2).In three years, two or three scenes at the end-of summer were available.Thus, in an ideal case the albedo trend of an of snow-covered areas in the scenes and/or sensor artefacts, less scenes were usually available to evaluate a temporal bare-ice albedo trend for single grid cells.We arbitrarily set the necessary number of scenes to 50%, thus at least ten albedo values are required for calculating the albedo trend of one individual grid cell.We used the non-parametric Mann-Kendall (MK) test (Mann, 1945;Kendall, 1975) to evaluate the confidence level of the trends (significant at the 95% / 90% / 80% level, or not significant).For grid cells with significant trends, the magnitude of the trend was determined based on linear regression through all available data points.Trends are given as albedo change per decade.For some analyses, we only use trends significant at the 95% confidence level for further interpretation and exclude less significant trends.

Spatially distributed shortwave broadband albedo
Figure 3 shows the spatio-temporal evolution of glacier-wide shortwave broadband albedo for Findelengletscher.The retrieval of meaningful albedo values is restricted by the quality of the surface reflectance data and, thus, the availability of realistic values in the individual bands needed for the narrow-to-broadband conversion.For Landsat TM/ETM, a saturation problem over snow-covered areas exists, resulting in missing values for these regions (years 1999-2012 in Figure 3).This problem is not present in the Landsat 8 data (years 2013-2016 in Figure 3).The striping in the Landsat ETM data also occurs in our albedo retrievals (e.g.09.09.2004 in Figure 3) and produces regions with missing data.Although, we applied a cloud removal algorithm, our results are still impacted by cloud shadows that are harder to detect without manual effort (e.g.18.08.2002 in Figure 3).However, the bare-ice area is almost always well represented and inferred albedo is realistic, hence allowing for a monitoring through time.
Generally, the average albedo values for the bare-ice surfaces are rather low, ranging from 0.17 to 0.30 for individual glaciers as a mean over the entire study period.For all 39 glaciers and over the entire study period we obtained a mean bare-ice surface albedo of 0.22.Extreme years with generally very high snowline altitude (2003,2011,2015) or very low snowline altitude (2013,2014) are linked to summers with exceptionally long and warm or rather cold and humid weather situations, and thus strong or weak ablation, respectively (Glaciological Reports, 2017).

Regional and glacier-wide trend in bare-ice albedo
We averaged mean albedo over the entire bare-ice area for each year and glacier to obtain a 39 individual time-series for the study period 1999 to 2016.In addition, the overall average series was evaluated as the mean bare-ice albedo of all glaciers and each year (Figure 4a).
Individual glaciers show considerable variations (up to 0.32 difference between minimum/maximum values) of mean bareice albedo between years; however, some glaciers show only minor interannual variability of about 0.06 such as for Grosser Aletsch.On average, the glaciers exhibit a range of 0.18 in minimum and maximum values in-between individual years.Due to these large interannual variations, no significant trends in average glacier-wide bare-ice albedo between 1999 and 2016 for The Cryosphere Discuss., https://doi.org/10.37 out of the 39 glaciers were found.Only two (Saleina, Trift) show slightly positive trends that are significant at the 95% confidence level according to the MK test.
The yearly values of the 39-glacier average albedo time series range from 0.18 to 0.27, with a mean of 0.22.There is no clear link between annual glacier mass balances and average bare-ice albedo for the 39 investigated glaciers (Figure 4b).As for the individual glaciers, no significant trend was found for the averaged time series over the period 1999 to 2016.As trends in bare-ice albedo for the entire ablation area of glaciers might be diluted by averaging over larger areas, or be affected by data uncertainty, we also evaluated the trend in albedo for all grid cells individually.For 133 km 2 (30%) of the entire surface area of all glaciers, trends were significant at the 80% level according to the MK test (Table 3).Thereof, 16 km 2 (12%) showed trends significant at the 95% confidence level or higher.Trends were classified according to their magnitude 10 for interpretation.Our classification is shown in Table 3. Classes with clear negative trends (class 1-3) are more abundant compared to classes with no clear or positive trends (class 4-7) at very high confidence level (95%) (Table 3).Thus, significant albedo trends at a confidence level of 95% in the bare-ice areas of the studied glaciers were only detected for grid cells with 9 The Cryosphere Discuss.a rather strong reduction of albedo over the 17 years.Surprisingly, more than 80% of all grid cells with significant albedo changes at the 95% confidence level showed negative trends: 20% exhibited changes of around −0.02 per decade, but 60% of cells showed trends more negative than −0.03 per decade.For some grid cells, also positive albedo trends significant at the 95% confidence level were detected however (about 15% or 2 km 2 ).
For most of the bare-ice area, the derived trends in albedo were only significant at low levels.Compared to the glaciers' overall ablation area only relatively few grid cells with trends significant at the 95% confidence level or higher (dark blue areas in Figure 5) are present.The cells with significant trends at high confidence levels are usually situated at the termini or along the lower margins of the glaciers and trends are mostly negative (cf.Table 3, Figures 5-7).The darkening can be attributed to different causes.At the glacier termini, an accumulation of fine debris due to the deposition of allochthonous material and/or melt-out of englacial debris is most likely.These materials, together with the presence of organic material, usually dark and humic substances, decrease local albedo values considerably and foster the growth of algae and bacteria (Hodson et al., 2010;Yallop et al., 2012;Takeuchi, 2013;Stibal et al., 2017).However, many of these effects and interactions are still unclear.Along the glacier margins an increase in debris cover due to small collapses or input of morainic material and, hence, a deposition of rather thick debris on the bare-ice is possible.Moreover, the appearance of debris-rich basal ice alongside the lower glacier margins due to the general glacier recession poses a further cause of local darkening (Hubbard and Sharp, 1995;Hubbard et al., 2009).Along the central area of the glacier tongue, particularly in the vicinity of medial moraines (e.g. in the case of Gornergletscher, Figure 7), a strongly negative albedo trend indicates an expanding medial moraine, changing the local area  from clean to (partly) debris-covered ice.In contrast, we also find significant positive albedo trends for some locations on the glacier tongues (see Figure 7).moraine, hence leading to a transition from debris-covered to clean ice with a higher albedo for certain grid cells.Lateral shifts of the position of medial moraines are possible for retreating glaciers (Anderson, 2000).

Uncertainty analysis
Our results are subject to uncertainties arising from errors in the input data, the albedo retrieval approach, the general data processing and the availability of data.An assessment of the uncertainties stemming from the input data, such as saturation or stripping issues in Landsat TM and ETM data, as well as the performance of the narrow-to-broadband formula (Equation 1) is beyond the scope of this study.Moreover, it became clear, that the most recent Landsat 8 data provide better results for snow-covered areas compared to Landsat TM and ETM (cf. Figure 3).However, as this study only focuses on bare-ice areas, this input data quality issue is not influencing the analysis of temporal albedo evolution.these areas.We minimized this effect by using glacier outlines updated to 2016 in order to exclude grid cells from the analysis that become ice-free towards the end of the study period.
To account for the uncertainty introduced by the temporal availability of remote sensing data we performed a comprehensive uncertainty analysis based on ten end-of summer Landsat 8 scenes acquired between 2013 and 2016 (Table 4).The analysis was performed for one glacier, Findelen, as most scenes were available for this glacier due to the overlapping coverage of this glacier by two different Landsat scenes (path/row 194/28 and 195/28).For the same grid cell and multiple satellite scenes acquired during the same year (1-5 weeks apart at maximum) we found an average variability in inferred albedo of 0.026 over all four investigated years (2013-2016) (Table 4).Assuming that bare-ice albedo remains constant over this short time period in reality, this value provides a direct uncertainty estimate for local satellite-retrieved albedo.
To assess the impact of local albedo uncertainty on the determination and the robustness of potential temporal trends, we  the computed average uncertainty of local albedo of 0.026 (average pixel number of 3500).The re-evaluation of the longterm albedo trends significant at the 80% level according to the MK test revealed that they were not affected by the random perturbation of the albedo values.Both a very similar area of the glaciers' bare-ice surfaces and distribution of trend magnitude was found in the perturbed datasets.However, for trends significant at the 95% confidence level or higher a slightly smaller area (11 km 2 ) was detected (c.f.Table 3).Within this area, the majority (77%) of all pixels is affected by negative trends, which is highly similar as obtained by the original albedo datasets (cf.Table 3).Thus, the inferred trends in local bare-ice albedo are considered to be robust despite the uncertainty in the albedo retrieval.

Snap-shot uncertainty and dependencies
In contrast to the quasi-continuous measurement setup of an automatic weather station, which is however only representative for a limited spatial extent (Ryan et al., 2017), airborne and spaceborne remote sensing datasets only represent a snap-shot in time.Hence, the temporal variability is only included to a certain degree and thus provokes a snap-shot uncertainty in surface albedo for evolution analyses.The meteorological conditions prior to the acquisition of the remote sensing imagery are highly important for the snap-shot uncertainty (Fugazza et al., 2016).Naegeli et al. (2017)  pared to a dataset acquired at the beginning of the melting period.However, this is only true, if meteorological characteristics between the individual acquisition dates are relatively constant.Snowfall or heavy rainfall events might significantly alter the ice surface conditions and the associated albedo values.While fresh snow increases the albedo strongly (e.g.Brock, 2004) and decreases the extent of the bare-ice area (Naegeli et al., 2017), rain can have a two-sided effect.A heavy precipitation event can lead to a short-term (between 1 to 4 days (Azzoni et al., 2016)) increase in albedo due to decreasing surface roughness and/or wash-out of fine debris present on the ice surface (between 5 to 20% according to Brock (2004) and Azzoni et al. (2016)), whereas light rainfall can cause the presence of a thin waterfilm on the glacier ice surface that absorbs radiation much stronger than the underlying ice and thus result in a decreased albedo.Similarly, a long-lasting phase with high air temperatures or intense shortwave radiation input during mid-day can lead to a permanent or temporary waterfilm on the ice surface that reduces reflectivity and thus shortwave broadband albedo considerably (Cutler and Munro, 1996;Jonsell et al., 2003;Paul et al., 2005).
Moreover, a remaining thin snow cover might cause slightly increased albedo values in the ablation area (still being in the typical range of glacier ice) that is difficult to be recognized with remote sensing data sets only (Naegeli et al., 2017).
The strong spatio-temporal variability of shortwave broadband albedo is dependent on various factors modulating the local glacier surface characteristics.Besides influencing the state of glacier surface albedo, these dependencies also impact on the changes in bare-ice albedo.For example, the surrounding lithology of a glacier determines (at least partially) the availability of fine debris material that can be transported by wind and water, and be deposited on the glacier ice, reducing its albedo considerably (Di Mauro et al., 2015, 2017;Azzoni et al., 2016).The investigation of the lithology surrounding the 39 individual glaciers and their overall albedo trend observed for the study period (Table 1) revealed that glaciers predominantly surrounded by less abrasive rocks (calcareous phyllites, limestones and marly shales, CERCHAR Abrasivity Index (CAI) 0-2 after Käsling and Thuro ( 2010)) exhibited a stronger negative albedo change of −0.05 per decade compared to glaciers that are located in an area of very to extremely abrasive rocks (−0.03 albedo change per decade; amphibolites, basic rocks, gneiss, granites, mica shists and syenites, CAI 2-6 after (Käsling and Thuro, 2010)).Thus, easily erodible rock-types provide more loose material that might be transported by wind and water on to the glacier surface and hence impact the bare-ice albedo.No relation between the albedo of the surrounding geology and the magnitude of the ice albedo change was evident however.These findings indicate the importance of the surrounding rocks as possible debris input source on a glacier, in particular as lateral moraines tend to become steeper and more instable due to general glacier recession in times of global atmospheric warming (Fischer et al., 2013), as well as their influence on the energy balance of the nearby glacier ice and snow surfaces.The discussion of these uncertainties and dependencies, highlights only parts of the complex spatio-temporal evolution of glacier surface albedo.While some influential factors mediating bare-ice albedo are obvious but challenging to quantify (e.g.meteorological conditions prior to the acquisition of data, micro-topography of the surface, etc.) others despite being quantifiable, are more ambiguous (glacier geometry, surface slope and aspect, surrounding lithology, etc.).Based on the presented results we therefore emphasize the need for further investigations of temporal and spatial dependencies of bare-ice albedo changes regarding various meteorological or geomorphological conditions and their interactions.

Temporal evolution of shortwave broadband albedo
Throughout the study period of 17 years, the spatial pattern of bare-ice albedo remained relatively stable for the 39 glaciers.
However, the extent of the bare-ice area exhibits a strong interannual variability (Table 2).This is mainly determined by local, temporary meteorological conditions varying strongly from year to year.The meteorological conditions prior to the acquisition dates are crucial as they considerably alter the surface characteristics and thus the observed broadband shortwave albedo.
Moreover, a prolonged ablation period has a strong impact on surface properties such as surface roughness and, hence, also impacts on glacier surface albedo (Cathles et al., 2011;Rippin et al., 2015).On smaller spatial and temporal scales, variations in glacier surface albedo are further evoked by meltwater redistribution of impurities (Hodson et al., 2007;Irvine-Fynn et al., 2011).
However, these complex surface-atmosphere interactions are still rather poorly constrained, in particular the temporal dimension, and further research in this area is needed.Nevertheless, as this study focused on end-of summer (August and September) scenes only, the relative variations between the individual years is comparable and robust.

Spatial scales of trends
A clear distinction between regional/glacier-wide bare-ice and local albedo changes is necessary if temporal trends are investigated.A negative trend in glacier-wide albedo (i.e.including both the ablation and the accumulation area) does not necessarily indicate a darkening of the glacier surface but rather a shift in snowline, or in other words, an enlargement of the bare-ice area relative to the total glacier surface.This effect is particularly pronounced in times of rising air temperatures and prolonged ablation periods.In contrast, a negative trend in bare-ice albedo can be an indicator of a darkening phenomenon due to an increased abundance of light-absorbing impurities (mineral dust, organic matter, algae, soot, etc.).Similarly, a lack of trends in bare-ice albedo change at the regional scale does not necessarily exclude the presence of significant trends at the local scale for individual grid cells.
In the frame of this study, we were unable to detect a spatially wide-spread, regional trend in bare-ice glacier albedo.
However, for certain regions of the glaciers, such as the lowermost glacier tongues or along the lower margins, significant negative trends were found.Hence, a clear darkening was observed at the local scale for a limited number of grid cells rather than for entire ablation areas.These findings are in agreement with published literature, but also show that findings of previous studies conducted at the local scale cannot be generalized for an entire glacier or the regional scale.For example, Oerlemans et al. ( 2009) observed a strongly negative albedo trend at a fixed location close to the terminus of Vadret dal Morteratsch, Switzerland, based on weather station data.They found an albedo reduction of 0.17, from 0.35 to 0.15, between 1995 and 2007.This trend is substantially higher than those detected in the present study over the period 1999 to 2016 (see also Table 3).
Other studies investigated glacier-wide albedo trends and found negative albedo trends of around 0.1 over the period 2000 to 2013 for Mittivakkat Gletscher, Greenland (Mernild et al., 2015) or up to 0.06 during the period 2000 to 2011 for nine glaciers in western China (Wang et al., 2014).However, the differing study periods, the varying observation scales and the impact of local characteristics on albedo changes make a direct comparison of albedo trends susceptible to misinterpretations.Based on 19 Landsat scenes over a 17-year study period, we assessed the spatio-temporal evolution of bare-ice glacier surface albedo for 39 glaciers in the western and southern Swiss Alps.Our results indicate that the spatial scale is crucial for the investigation of albedo trends and the detection of a potential darkening effect that is often referred to in recent literature (Takeuchi, 2001;Oerlemans et al., 2009;Dumont et al., 2014;Wang et al., 2014;Mernild et al., 2015;Tedesco et al., 2016).
While we did not find a darkening of bare-ice glacier areas at the regional scale or averaged for the ablation areas of individual glaciers, significant negative albedo trends were revealed at the local scale.These individual grid cells or small areas were mainly located at the glacier termini or along the lower glacier margins.
We emphasize the importance of the snap-short uncertainty -limited numbers of end-of-summer scenes demand recognition.Specifically, the meteorological conditions preceding the acquisition of the satellite data can influence albedo and so need to be considered.Although, only snap-shots of glacier surface albedo are available, the almost two-decade long time-series indicate a highly significant darkening trend for about 10% of the ablation areas which typically ranges from about −0.05 to −0.03 per decade.For these areas, the positive ice-albedo feedback enhanced melt rates and is expected to be enforced in the near future.Even though the darkening of glacier ice has been found to occur over only a limited area of the investigated glaciers, the projected enlargement of bare-ice areas characterised by low albedo coupled with the predicted prolongation of the melt season will most likely strongly impact on the glacier surface energy balance and substantially enhance ice mass loss.

Figure 1 .
Figure 1.Overview of studied glaciers located in the western and southern Swiss Alps.The black frames mark the extent of detailed figures.The inset indicates the position of the study site in Switzerland.

Figure 2 .
Figure 2. Flow chart of the methodology applied to evaluate the surface type (snow/ice) of every glacier grid cell based on the derived shortwave broadband albedo.For further explanations please see text.

Figure 4 .
Figure 4. (a) Time series of mean bare-ice albedo of all 39 glaciers (grey dots) and their overall average (black dots with dashed line).(b) Annual mass balance anomalies relative to 2006-2016 for the Rhone catchment (Glaciological Reports, 2017).

Figure 5 .
Figure 5. Confidence levels of bare-ice albedo trends over the study period 1999 to 2016 according to the MK test.See Figure 1 for the location of the different panels in the Swiss Alps.

Figure 6 .
Figure6.Classified albedo trends per decade for all grid cells with trends significant at the 80% confidence level or higher according to the MK test.Averages for areas with bare-ice albedo trends significant at the 95% confidence level are given for selected large glaciers.See Figure1for the location of the different panels in the Swiss Alps.

10Figure 7 .
Figure 7. (a) Close-up of bare-ice albedo trends per decade significant at the 95% confidence level or higher for the tongue of Aletsch, and (b) time-series of bare-ice albedo between 1999 and 2016 for ten randomly selected points on the terminus (crosses in (a)) including a linear fit (dashed purple, r = −0.6).

10
randomly perturbed the distributed bare-ice albedo values of every grid cell and scene, and for all 39 individual glaciers with 13 The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-18Manuscript under review for journal The Cryosphere Discussion started: 29 January 2018 c Author(s) 2018.CC BY 4.0 License.
highlighted this fact by cross-comparing albedo products from three different sensors with acquisition times within one week.If glacier-wide albedo is compared, a dataset acquired later in the ablation season is expected to show a larger bare-ice area characterised by low albedo values com-The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-18Manuscript under review for journal The Cryosphere Discussion started: 29 January 2018 c Author(s) 2018.CC BY 4.0 License.

Table 2 .
Overview of Landsat scenes used.The bare-ice area is given in km 2 and relative to the total study area of 451 km 2 .Cloud coverage is given for bare-ice areas only, and is also evaluated relative to the overall bare-ice area of the respective scene.

Table 3 .
Overview of classes of bare-ice albedo trends for individual grid cells between 1999 and 2016 corresponding to the confidence levels of 80% and 95% according to the MK test.Numbers refer to the sum of all bare-ice grid cells of all 39 study glaciers.

Table 4 .
Overview of scenes used in the uncertainty analysis.Px refers to the number of pixels that were used to derive uncertainty.Mean (αmean), minimum (αmin) and maximum (αmax) albedo, as well as the mean (σmean) standard deviation of point-based bare-ice albedo for each individual scene pair or triple per year are given.