Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data

Less than 0.25 % of the 250 000 glaciers inventoried in the Randolph Glacier Inventory (RGI V.5) are currently monitored with in situ measurements of surface mass balance. Increasing this archive is very challenging, especially using time-consuming methods based on in situ measurements, and complementary methods are required to quantify the surface mass balance of unmonitored glaciers. The current study relies on the so-called albedo method, based on the analysis of albedo maps retrieved from optical satellite imagery acquired since 2000 by the MODIS sensor, on board the TERRA satellite. Recent studies revealed substantial relationships between summer minimum glacierwide surface albedo and annual surface mass balance, because this minimum surface albedo is directly related to the accumulation–area ratio and the equilibrium-line altitude. On the basis of 30 glaciers located in the French Alps where annual surface mass balance data are available, our study conducted on the period 2000–2015 confirms the robustness and reliability of the relationship between the summer minimum surface albedo and the annual surface mass balance. For the ablation season, the integrated summer surface albedo is significantly correlated with the summer surface mass balance of the six glaciers seasonally monitored. These results are promising to monitor both annual and summer glacier-wide surface mass balances of individual glaciers at a regional scale using optical satellite images. A sensitivity study on the computed cloud masks revealed a high confidence in the retrieved albedo maps, restricting the number of omission errors. Albedo retrieval artifacts have been detected for topographically incised glaciers, highlighting limitations in the shadow correction algorithm, although interannual comparisons are not affected by systematic errors.


Introduction
Mountain glaciers represent only 3 % of the ice volume on the Earth but contribute significantly to sea level rise (e.g.Church et al., 2013;Gardner et al., 2013;Jacob et al., 2012).In addition, millions of people partly rely on glaciers, either for drinking water or agriculture or due to related glacier hazards (Baraer et al., 2012;Chen and Ohmura, 1990;Immerzeel et al., 2010;Kaser et al., 2010;Sorg et al., 2012;Soruco et al., 2015).The surface mass balance (SMB) of glaciers is directly driven by the climate conditions; consequently, glaciers are among the most visible proxies of climate change (Dyurgerov and Meier, 2000;Haeberli and Beniston, 1998;Oerlemans, 2001;Stocker et al., 2013).Measuring and reconstructing glacier SMB therefore provides critical insights into climate change both at global and regional scales (Oerlemans, 1994).
Systematic SMB monitoring programmes began in the late 1940s and early 1950s in most of the European countries (e.g.France, Norway, Sweden, Switzerland).Gradually, more glaciers have become monitored, reaching the present worldwide figure of 440.However, this represents only a small sample of the nearly 250 000 inventoried glaciers worldwide (Pfeffer et al., 2014).Among the existing methods to quantify changes in glacier SMB, the well-established glaciological method has become a standard, widely used worldwide, yielding most of the reference datasets (World Glacier Monitoring Service, WGMS; Zemp et al., 2015).Based on repeated in situ measurements, this method requires intensive fieldwork.However, this method is unable to reconstruct the SMB of unmonitored glaciers.The Global Terrestrial Network for Glaciers (GTN-G) aims at increasing substantially the number of monitored glaciers to study regional climate signal through changes in SMB.To reach this objective, the development of methods complementary to the ground-based glaciological method is therefore required.Since the 1970s, several methods have taken advantage of satellite imaging to compute changes in glacier volume (Kääb et al., 2005;Rabatel et al., 2017;Racoviteanu et al., 2008).Several glacier surface properties have thus been used as proxies for volume fluctuations -changes in surface elevation from differencing digital elevation models (DEMs) (e.g.Belart et al., 2017;Berthier et al., 2016;Gardelle et al., 2013;Ragettli et al., 2016;Shean et al., 2016); end-of-summer snow line elevation from high spatial resolution optical images (Braithwaite, 1984;Chinn et al., 2005;Meier and Post, 1962;Mernild et al., 2013;Rabatel et al., 2005Rabatel et al., , 2008Rabatel et al., , 2016;;Shea et al., 2013); mean regional altitude of snow from low spatial resolution optical images (Chaponniere et al., 2005;Drolon et al., 2016); or changes in the glacier surface albedo from high temporal resolution images (Brun et al., 2015;Dumont et al., 2012;Greuell et al., 2007;Greuell and Knap, 2000;Shea et al., 2013;Sirguey et al., 2016).Widely used over icecaps or large ice masses, satellite derived DEMs cannot yet be confidently used to compute annual or seasonal SMB of mountain glaciers, although recent studies have revealed promising results for determining SMB changes of large mountainous glacierized areas (Belart et al., 2017;Ragettli et al., 2016).The method based on the correlation between the regional snow cover and glacier SMB has shown satisfying results to retrieve seasonal SMB, especially for the winter period.This method was used for the quantification of 55 glaciers SMB in the European Alps over the period 1998-2014 (Drolon et al., 2016).The method based on the identification on high spatial resolution optical images of the end-of-summer snow line altitude has shown encouraging results in the French Alps, multiplying by 6 the available long-term annual SMB time series (Rabatel et al., 2016), but needs to be automated to compute glacier SMB at regional scales.In addition, monitoring glacier surface properties on the daily or weekly basis and over large glacierized regions is still challenging with high spatial resolution images.The current study is based on the albedo method used in Dumont et al. (2012), Brun et al. (2015) and Sirguey et al. (2016).Images from the MODerate resolution Imaging Spectroradiometer (MODIS) are processed to compute daily albedo maps of 30 glaciers in the French Alps over the period 2000-2015.Then we rely on the methodological framework proposed by Sirguey et al. (2016) on Brewster Glacier (New Zealand), looking at the relationships between annual and seasonal SMB and the glacier-wide averaged sur-face albedo α.Our overall objective is to study the relationships between glacier SMB and albedo by (i) reconstructing the annual albedo cycle for 30 glaciers in the French Alps for the period 2000-2015, (ii) linking the albedo signal to the summer components of the SMB as well as to its annual values for 6 and 30 glaciers, respectively, and (iii) assessing the sensitivity of the retrieved albedo towards tuning parameters (cloud coverage threshold for images processing, reliability of detected shadows).Section 2 presents the available SMB datasets used for the comparison and describes briefly the in situ automatic weather stations (AWS) used to assess the quality of MODIS-retrieved albedo.The method to retrieve albedo maps is described in Sect.3. Results are presented and discussed in Sects.4 and 5.The conclusion gathers the main results of the study and provides perspectives for future works.
2 Study area and data

Site description
The study focuses on 30 glaciers located in the French Alps (Fig. 1).Each glacier can be classified as a mountain glacier, extending over an altitudinal range from around 1600 m a.s.l.(Argentière and Mer de Glace glaciers) to 4028 m a.s.l.(Blanc Glacier), and located between the coordinates 44 • 51 to 46 • N and 6 • 09 to 7 • 08 E. The cumulative glacial coverage considered in the present study is 136 km 2 , i.e. half of the glacier surface area covered by 593 inventoried glaciers over the French Alps for the period 2006-2009 (Gardent et al., 2014).
Studied glaciers have been selected following four criteria related to the availability of field data and remote sensing constraints, namely (i) the annual glacier-wide SMB for the study period had to be available, (ii) the glacier surface area had to be wide enough to allow robust multi-pixel analysis, (iii) the glacier had to be predominantly free of debris to allow remotely sensed observations of the albedo of snow and ice surfaces, and (iv) summer SMB records had to be available to consider summer variability.Finally, 11 glaciers have been selected in the Ecrins Range, 14 in Vanoise and 5 in Mont Blanc (Fig. 1, and listed Table 1).

MODIS satellite images
The MODIS sensor, on board the TERRA-EOS/AM-1 satellite has been acquiring near-daily images of the Earth since 25 February 2000.With 36 spectral bands ranging from 0.459 to 14.385 µm, and spatial resolution ranging from 0.25 to 1 km depending on the spectral band, MODIS is nowadays one of the most used optical sensors for land surface observations.Because of its short temporal revisit time, its long acquisition period and its moderate resolution, images from MODIS are the most suitable for the present work.We there-  1).The four AWS used in the present study were set up on Saint-Sorlin Glacier (no.20).Adapted from Rabatel et al. (2016).

Surface mass balance data
In the French Alps, six glaciers allow both the summer and annual analyses to be conducted, due to the availability of summer SMB data (b s ) obtained from in situ measurements with the glaciological method (unpublished data, LGGE internal report, listed in Table 1).In addition, glacier-wide annual SMB values of the 30 studied glaciers were computed by Rabatel et al. (2016) using the end-of-summer snow line measured on optical remote-sensing images and the glacierwide mass change quantified from DEM differencing.
For the six glaciers where glacier-wide annual SMB is available from the two methods, i.e. in situ and satellite measurements, the average of the two estimates was used to calibrate and evaluate the albedo method, in order to derive for each glacier a single relationship SMB vs. computed albedo.We do not discuss here the differences between the consid-ered datasets because these differences have been investigated by Rabatel et al. (2016).

In situ albedo measurements
Albedo measurements acquired punctually using an AWS on Saint-Sorlin Glacier have been used to evaluate the MODISretrieved albedo.In situ albedo measurements were available for three periods in the ablation zone (July-August 2006; June-August 2008; June-September 2009) and for one period in the accumulation zone (June-September 2008).Albedo data from these AWS have been calculated as the ratio of the reflected to incident shortwave radiation (0.3-2.8 µm) using two Kipp and Zonen pyranometers.With a potential tilt of the instrument with respect to surface melting and the intrinsic sensor accuracy (±3 %, Six et al., 2009), the calculated albedo at the AWS shows a ±10 % accuracy (Kipp and Zonen, 2009;Dumont et al., 2012 3 Methods

MODImLab products
MODIS L1B images were processed using the MOD-ImLab toolbox (Sirguey, 2009).Image fusion between MOD02QKM bands 1 and 2 at 250 m resolution and MOD02HKM bands 3-7 at 500 m resolution allows seven spectral bands at 250 m resolution to be produced (Sirguey et al., 2008).Then, atmospheric and topographic corrections are applied that include multiple reflections due to steep surrounding topography (Sirguey, 2009).Various products are derived from the corrected ground reflectance including snow and ice surface albedo (Dumont et al., 2012).As recommended by Dumont et al. (2012), the white-sky albedo (estimated value of the surface albedo under only diffuse illumination) is considered.The use of an anisotropic re-flection model for snow and ice has been preferred to the isotropic case, due to its closer agreement with in situ measurements (Dumont et al., 2012).The MODImLab toolbox also produces sensor geometrical characteristics at the acquisition time such as the solar zenith angle (SZA) and the observation zenith angle (OZA) used for post-processing the images (Sect.3.4).The MODImLab cloud detection algorithm is more conservative than the original MODIS product (MOD35), and has been preferred as recommended in Brun et al. (2015).
According to Dumont et al. (2012) and further assessed by Sirguey et al. (2016) the overall accuracy of MODIm-Lab albedo product under clear-sky conditions is estimated at ±10 %.
To mitigate the impact of shadows over the glaciers, MODImLab uses a DEM from the Shuttle Radar Topography Mission (SRTM -90 m resolution -acquired in 2000) to estimate the sky obstruction by the surrounding topography and to correct the impact of shadows (see Sirguey et al., 2009).The algorithm implemented in MODImLab is fully described by Sirguey et al. (2009) and was inspired by Dozier et al. (1981) and Dozier and Frew (1990) for the sky obstruction factor processing (Horizon and Vsky in Sirguey et al., 2016), and Richter (1998) for the correction of shadows.It is first computed at 125 m resolution, providing Booleantype products of self and cast shadows per pixel.Results are then averaged and aggregated to 250 m resolution, producing a sub-pixel fraction of shadow (further detailed in Sirguey et al., 2009).Finally, MODIS data processed with MOD-ImLab provide, among others, near-daily maps of white-sky albedo at 250 m resolution together with cloud masks and cast and projected shadows.
Albedo maps have been processed for 5068 images for the Ecrins range, 4973 for Mont Blanc and 5082 for Vanoise over the period 2000-2015.Only images acquired between 09:50 and 11:10 UTC (+2 h in summer for local time conversion) were selected to get minimum SZA and limit projected shadows of surrounding reliefs.

Glacier masks
Following Dumont et al. ( 2012) and Brun et al. (2015), we manually created raster masks of the 30 glaciers, based on the glaciers' outlines from the 1985-1987 (Rabatel et al., 2013) and high spatial resolution (6 m) SPOT-6 images from 2014.All debris-covered areas, together with mixed pixels (rocksnow/ice) have been removed to capture only the snow/ice albedo signal.The resulting number of pixels per glacier is listed in Table 1.

Basis of the method
For one glacier in the Alps (Dumont et al., 2012), two in the Himalayas (Brun et al., 2015) and one in the Southern Alps of New Zealand (Sirguey et al., 2016), the summer minimum glacier-wide averaged albedo (α min a ) has been significantly correlated with the glacier-wide annual SMB.The relationship between α min a and glacier-wide SMB results from the fact that solar radiation is the main source of energy for melting snow and ice, both at the surface and within the first centimetres below the surface (Van As, 2011).But this is not sufficient to explain why averaged surface albedo is suitable for monitoring glacier SMB.
If we consider a temperate glacier in the mid-latitudes, its surface is fully covered by snow in winter, leading to high and uniform surface albedo (α min ≈ 0.8 in Cuffey and Paterson, 2010).During the ablation season, the accumulation area is still covered with snow conversely to the ablation area where the ice is exposed and sometimes covered by debris.The overall albedo of the glacier surface is therefore decreasing over the course of the ablation season, providing information on the ratio of these two areas.The ratio between the size of the accumulation zone and the entire glacier, called the accumulation-area ratio (AAR) has often been used as a predictor of SMB both qualitatively (LaChapelle, 1962;Meier and Post, 1962;Mercer, 1961) or quantitatively (Dyurgerov et al., 2009).Therefore, assessing α min a provides insight into the relative share of the exposed ice and the snow-covered areas at the end of the ablation season, also quantified by the AAR.

From annual to summer SMB
In this study, α min a has been computed for the 30 glaciers in order to validate the method at a regional scale.Only the α min a occurring in summer have been considered because minimum values out of the summer period are artifacts.Then, α min a has been directly correlated with available annual SMB data (listed in Table 1).
Following the work by Sirguey et al. (2016) on Brewster Glacier, a similar approach has been used in order to validate the method at a summer scale but only on six glaciers (within our sample of 30) for which the summer SMB are available.Conversely to Sirguey et al. (2016), the summer SMB b s has been compared to the integrated albedo signal α int s during the entire ablation season (1 May to 30 September) computed as follows and illustrated in Fig. 2 Integrated summer albedos enable us to account for snowfall events that can occur during the ablation period (punctual high albedos).As an example, a strong summer snowfall event leading to a rather persistent snow coverage of the glacier will "feed" the integrated albedo, and physically reduces the glacier melting, which has an impact on the SMB (Oerlemans and Klok, 2004).The method therefore accounts for snowfall events to retrieve the glacier summer SMB.To compare each year together and remove the impact of the variable integration time period for each glacier, α int s has been divided by the number of integrated days.

Data filtering
MODIS offers the opportunity to get daily images, but retrieving daily maps of Earth surface albedo remains challenging.Indeed, various sources of error require filtering of the available images in order to only capture physical changes of the observed surface and not artifacts.Clouds are known to be a major problem in optical remote sensing of the Earth surface especially in the case of ice-and snow-covered surfaces.Even though some algorithms exist to differentiate clouds and snow-covered areas (e.g.Ackerman et al., 1998;Sirguey et al., 2009), omission errors are difficult to avoid, leading to erroneous albedo of the surface.
In this study, all images with a presence of cloud greater than 30 % of the total glacier surface area have been discarded.This threshold is higher than that chosen in Brun et al. (2015) on the Chhota Shigri Glacier (20 %), and we thus discuss in Sect.5.1 the impact of the computed cloud threshold on the derived albedo results.When determining α min a , 0 % of cloud cover has been imposed as a condition, and a visual check for each year and each glacier has been performed.Snapshots from the fusion of MODIS bands 1-3 and from bands 4-6 (Sirguey et al., 2009) have been used to visually check the images, together with images from other satellites (mostly from the Landsat archive) and pictures and comments from mountaineering forums.This last step, although laborious when studying 30 glaciers, allowed the identification of the summer minimum to be improved.A visual check of these images also confirms that projected shadows of clouds are not affecting the albedo maps of summer minimum.Another source of error is the impact of the OZA.As mentioned in Sirguey et al. (2016), accuracy of the MODIS retrieved albedo strongly decreases for viewing angles above 45 • as pixel size increases from 2-to 5-fold from OZA = 45 • to 66 • (Wolfe et al., 1998).This phenomenon is accentuated when observing steep-sided snow/ice surfaces, surrounded by contrasted surfaces (rocks, forests, lakes).This distortion could lead to capturing the mean albedo of a glacier plus its surroundings.As a result of this, we decided to filter the images according to their OZA angle, as further described Sect.4.1.

Retrieved albedo assessment
A quantitative evaluation of the retrieved albedo has been performed with AWS deployed on Saint-Sorlin Glacier.Measurements have been synchronized between punctual albedo for MODIS and a 2 h averaged albedo around MODIS acquisition time for the AWS.It is worth recalling some differences between the in situ measured albedo data and the one retrieved using MODIS.The downward-facing pyranometer stands at around 1 m above the surface, corresponding to a monitored footprint of ca.300 m 2 (theoretical value for a flat terrain) while the pixel area of MODImLab products matches 62 500 m 2 .Quantified albedos from each method are therefore not representative of the same area.On the other hand, incoming radiation data are extremely sensitive to a tilt of the sensor located on the AWS, and maintaining a constant angle throughout the monitoring period remains challenging, especially during the ablation season.For instance, a tilt of 5 • of the pyranometer at the summer solstice can increase by 5 % the error on the irradiance measurement (Bogren et al., 2016).No sensor tilt was deployed on the AWS, thus preventing the application of tilt correction methods (e.g.Wang et al., 2016).Nonetheless, regular visits allowed us to maintain the sensor horizontal and to limit errors in the irradiance measurements.
Figure 3 illustrates the comparison between the retrieved and measured albedos at the AWS locations for various OZA classes.One can note minor differences between the data plotted in Fig. 3 and those presented in Dumont et al. (2012, Fig. 2).These differences are related to changes in the MOD-ImLab algorithm and different computation of the in situ albedo, integrated over a 2 h period in the current study.
In Fig. 3, the spread between MODIS and AWS albedos is higher for low albedos (i.e.ablation area).This is related to the footprint difference as described earlier, accentuating the albedo differences when monitoring heterogeneous surfaces (snow patches, melt pounds. . .), which are even more pronounced in summer.One can also note that MODIS albedo often overestimates the AWS albedo value.This overestimation could be explained by (1) the MOD-ImLab albedo retrieval algorithm.(Under-estimation of the incoming radiation computed in the MODImLab algorithm would lead to overestimated retrieved albedo values, and in addition the atmospheric corrections used to compute the incident radiation could be hypothesized as source of errore.g.modelled transmittance through a simplified computed atmosphere; see Sirguey et al. (2009) for further description), or (2) the AWS albedo measurements.Indeed, view angles of AWS pyranometers (170 • ) could influence the retrieved albedo by monitoring out-of-glacier features (e.g.moraines, rock walls,), resulting in underestimated albedo values.However, it is worth noting that most of the points are within the combined uncertainty of both sensors and these differences .The continuous grey line illustrates the 1 : 1 relationship between AWS and MODIS retrieved albedo.Thin and dotted lines represent the combined uncertainties on both AWS-and MODIS-retrieved albedo (absolute value of 10 % for each), only accounting for intrinsic sensor accuracy and not for errors related to the acquisition context, e.g.size of the footprint.in albedo retrieved from MODIS and the AWS are thus hard to interpret.Finally, Fig. 3 shows substantial differences between OZA < 10 • and other OZA classes.For OZA < 10 • , MODIS albedos better agree with AWS albedos than for the three other classes.Integrating MODIS images with OZA > 10 • substantially deteriorates the agreement with AWS albedos (in term of r 2 , RMSE and the slope P MODIS 1 ), especially on "narrow" targets as alpine mountain glaciers.We therefore chose to prioritize images acquired with low OZA to avoid detection of non-glacierized surfaces.Therefore, four classes of images have been selected following the criteria presented in Table 2.
For the rest of the computation, the absolute ±10 % accuracy per pixel estimated in Dumont et al. ( 2012) has been considered.We determined the uncertainty on α by accounting for the spatial variability of the albedo signal within the glacier and considering that our sets of pixels are independent from each other (Eq.2): where σ stands for the standard deviation (SD) of the pixels albedo with N the number of pixels.

Temporal variability of the albedo signal
Using the "step-by-step" filtering procedure explained in Sect.3.4, the ∼ 16-year albedo cycle of each of the  The green dots spot for each summer the minimum average albedo, and have been manually checked for all years and glaciers.Dashed red and blue lines stand for the beginning of the defined ablation and accumulation seasons (May and 1 October respectively).
30 glaciers was obtained.Figure 4 illustrates the entire albedo time series for Saint-Sorlin Glacier over the period 2000-2015.We observed that the albedo decreases from the beginning of summer (dashed red line), reaching α min a in August/September and rising again at the end of September.This cyclicality is a proxy of surface processes.The snow cover decreases at the beginning of summer until reaching its lowest extent, and finally increases again with the first snowfall in late summer to reach its maximum extent in winter/spring.
The periodicity of the albedo signal is however not so well defined for some of the studied glaciers.For instance, Argentière Glacier exhibits a severe drop of α in winter, reaching values as low as summer minima (α ∼ 0.4).The observed drop of albedo in winter occurs during more than 1 month centred on the winter solstice (21 December) and is observed for nine glaciers (Argentière, Baounet, Casset, Blanc, Girose, Pilatte, Vallon Pilatte, Tour and Sélé glaciers).These glaciers are located within the three studied mountain ranges and have the common characteristic of being very incised, with steep and high surrounding faces.We studied the albedo series as a function of the SZA to reveal possible shadowing on the observed surfaces.Figure 5 displays the same cycle as Fig. 4 for Argentière Glacier but providing information about SZA.As a reminder, the MODImLab white-sky albedo is independent of the illumination geometry but the computed albedo for each pixel can be subject to shadowing from the surrounding topography.
Two main observations stand out from the winter part of the cycle in Fig. 5: (i) most of MODIS α severely decrease under α = 0.6 for SZA greater than 60 • corresponding to November to January images, and (ii) these drops are not systematic and we rather observe a dispersion cone than a welldefined bias.As there are no physical meanings to systematic change of the surface albedo during a part of the winter period and owing to the fact that this dispersion is only observed for topographically incised glaciers, these decreases in albedo have been considered as artifacts.These observations led us to perform a sensitivity study on the validity of the shadow mask produced by MODImLab, and to study the impact of these shadows on the retrieved glacier-wide albedo (see Sect. 5.2).

α min a and annual SMB
The summer minimum average albedo for each year and each glacier has been linearly correlated with the glacier-wide annual SMB. Figure 6 illustrates the relationship between α min a and b a for Blanc Glacier.Error bars show the dispersion of the SMB dataset for each year, and from the glacier intrinsic variability of the albedo signal on the day of α min a acquisition.For the glaciers where the glacier-wide annual SMB is available from the SLA method, the uncertainty is about ±0.22 m w.e. on average (ranging from 0.19 to 0.40 m w.e.depending on the glacier; Rabatel et al., 2016).
Twenty-seven glaciers show significant correlations (see Table 1 for full results) if considering a risk of error of 5 % (according to a Student's t test), which confirms the robust correlation between α min a and b a .However, the linear correlation has no statistical significance for three glaciers with r 2 < 0.25.A possible explanation is the high number of removed images in summer due to manually checked thin overlying clouds not detected by the MODImLab cloud algorithm.
Looking at the 27 glaciers for which significant relationships have been found, 2001 is regularly identified as an outlier.According to existing SMB datasets, 2001 is the only year of the period 2000-2015 for which the annual SMB has been positive for all the studied glaciers (+0.80 m w.e.yr −1 on average).
To predict correctly the SMB values for the year 2001 using the albedo method, monitored minimum glacier-wide average albedo would need to be extremely high (often greater than 0.7, i.e. 0.83 and 0.95 for Rochemelon and Vallonnet glaciers, respectively), to match the regression line derived from other years of the time series (Table 1).Taking into consideration snow metamorphism during the summer pe- The Cryosphere, 12, 271-286, 2018 riod, melting at the surface and possible deposition of debris or dusts, monitoring such high albedo values averaged at the glacier scale is unrealistic.As removing 2001 from the time series does not increase the number of glaciers for which the correlation is significant, 2001 has been conserved in the time series.However, this observation reveals a limitation of the albedo method by underestimating the annual SMB value for years with very positive annual SMB.

α int s and summer SMB
Studying the integral of the albedo signal during the ablation season can provide insight into the intensity of the ablation season and thus into the summer SMB b s .As described in Sect.3.3.2,α int s has been computed and connected to the in situ b s .Figure 7 illustrates the results for Saint-Sorlin Glacier.
Saint-Sorlin Glacier, together with the five other seasonally surveyed glaciers, showed a significant correlation between the two observed variables (from r 2 = 0.46 to   r 2 = 0.94 with an error risk < 5 %, all statistics detailed in Table 1).Conversely to α min a , the α int s is slightly more robust to the presence of undetected clouds as its value does not rely on a single image.The lowest correlation has been found for Talèfre Glacier.The latter accounts for a relatively large debris-covered tongue that has been excluded when delineating the glacier mask (see Supplement).Consequently, the low correlation could be partly explained by this missing area, considered in the glaciological method but not remotely sensed.To conclude, α int s has been significantly correlated with b s and is therefore a reliable proxy to record the ablation season.

Discussion
In this section, we first discuss the impact of the threshold applied to the cloud cover fraction on the obtained results.Then, a sensitivity study focused on the algorithm correcting the shadows is presented.We finally express the main limitations and assessments of the albedo method.

Cloud coverage threshold
As stated in Sect.3.4, a value of 30 % of cloud coverage over the glacier mask has been defined as the acceptable maximum value for considering the albedo map of the day.We computed a sensitivity study on the impact of this threshold on the value of the obtained correlations between the integrated summer albedo and the in situ summer SMB.The summer period has been chosen as it represents the period when the albedo of the glacier is the most contrasted, between bare ice and snow/firn.The glacier-wide average albedo in this period is therefore more sensitive to possible shading of a part of the glacier.Figure 8 illustrates the results for the six seasonally surveyed glaciers.The used value of the allowed cloud coverage appears not to have a substantial impact on the correlation.This observation implies that the MODImLab cloud product is reliable enough to only compute surface albedo and to avoid too frequent misclassification between the clouds and the surface.It also suggests that removing too many images because of partial cloud cover removes information about the glacier-wide average albedo variability.However, allowing all images, even when the glacier-wide average albedo is computed on only 10 % of the glacier (90 % of detected cloud coverage), does not reduce significantly the correlation for most of the six glaciers.
Nevertheless, hypothesizing that the glacier-wide average albedo of a small fraction of the glacier (e.g. 10 %) is suitable to represent the entire glacierized surface is questionable.It therefore depends on the size of the observed glacier, where 10 % of a glacier of 3 and 30 km 2 do not have the same meaning, but also on the delineated mask (ablation area not entirely considered because of debris coverage).The summer-integrated albedo is also highly dependent on the time gap between valid seful images.In other words, if an image has an "anomalous" glacier-wide average albedo because of high cloud coverage, the impact on the integrated value will be smaller if "normal condition" albedos are monitored at nearby dates.
The average number of available images per year does not largely differ between the various computed cloud coverage thresholds.It varies on average from 95 to 123 images per summer period for respectively 0 and 100 % cloud coverage threshold.Intermediate values are 106, 111 and 116 images per summer for 30, 50 and 75 % cloud coverage threshold, respectively.The difference in significance of r 2 (according to a Student's t test) between opting for 0 and 100 % is almost negligible, and choosing the best cloud threshold value is rather a compromise between the number of used images and the resulting correlation with glacier-wide SMB.We finally concluded that selecting a cloud coverage threshold of 30 % presents the best determination coefficients between the integrated summer albedo and the summer balance for most of the six glaciers without losing too much temporal resolution.

Assessment of the impact of shadows on retrieved albedos
In light of the documented dispersion on α during some of the winter months on several studied glaciers (Sect.4.2), sensitivity of the MODIS retrieved albedo against correction of shadows had been assessed.This work has only been conducted on the 250 m resolution raster products and specifically on the cast shadow product, because self-shadow corrections can be considered as reliable enough as they are only related to the DEM accuracy.We thus defined a pixel as "corrected" when at least one of its sub-pixels was classified as shadowed.From then on, two glacier-wide albedos α have been defined: (i) α non-cor computed on non-corrected pixels only, classified as non-shadowed, and (ii) α of both corrected and non-corrected pixels, equal to the glacier-wide average albedo.Figure 9 illustrates the difference between α non-cor and α as a function of the percentage of corrected pixels over the entire glacier.The study was performed on Argentière Glacier (111 pixels) that exhibited large α artifacts in winter (Fig. 5).The inner diagram allows us to emphasize the annual "cycle" of modelled shadows, contrasted between nearly no cast shadows in summer and an almost fully shadowed surface in winter.We represent the 1 SD of α, averaged by classes of 5 % corrected pixels.In other words, it illustrates the mean variability of the glacier-wide surface albedo.Therefore, for images with α non-cor − α within the interval defined by 1 SD of α, errors resulting from the correction algorithm are smaller than the spatial variability of the glacierwide albedo glacier.We also selected only significant values, following a normal distribution of the averaged α.Consequently, only values at ±1σ (68.2 %) in terms of percentage of corrected pixels have been retained (i.e. when the relative share of corrected pixels ranged from 15.9 to 84.1 %).Between 0 and 15.9 %, α non-cor and α are not sufficiently independent because of the low number of corrected pixels, and beyond 84.1 %, α non-cor is computed over a too small number of pixels.As a consequence, even if the albedo correction in the shadowed parts of the glacier could be improved, most of the errors related to this correction do not depreciate the results.Above 80 % of corrected pixels (December to early February), differences between α non-cor and α exceed the monitored spatial variability of α.These anomalies are at the root of the observed artifacts in Fig. 5 caused by the severe drops of albedos and described Sect.4.2.
In addition, a seasonality in the albedo signal can be observed -with α non-cor − α > 0 in early spring (February to April) and α non-cor − α < 0 in summer and autumn (June to November).This could be explained by different localizations of shadowed area for a given ratio of corrected pixel.As an example, a glacier could have in October a snowand shadow-free snout and a fresh snow-covered and shadowed upper section.This configuration would induce a negative difference, as we observe from June to November.Conversely, this glacier could present in March (same ratio of corrected pixels as October) a complete snow coverage, leading to a smaller difference between α non-cor and α (< 0.1), that could even result in a positive difference, as we observe from February to April.
Finally, observed albedo artifacts in winter are most likely due to the correction of shadows.On the other hand, correcting shadows accurately and consistently is extremely challenging.As illustrated by Fig. 9, a way to confidently consider the albedo signal is to exclude values with too large a share of corrected pixels.However, because of the interannual approach carried out in this study, such a systematic artifact is not depreciating the results but would be a major issue in studies focused on albedo values themselves (e.g.maps of snow extent).

Limits of the albedo method
In agreement with Dumont et al. (2012) and Brun et al. (2015), retrieving the glacier annual SMB from albedo summer minima proves to be an efficient method.Low correlations often result from high and persistent cloud coverage during summer, reducing the chance of spotting the albedo summer minimum.For SMB reconstruction purpose, a future line of research could rest upon linking morpho-topographic features of the glacier, such as glacier surface area, mean altitude or slope, to the regression coefficients of both annual and seasonal SMB vs. albedo relationships, giving the opportunity to establish analogy between monitored and unmonitored glaciers.Tests have been carried out, but no significant and satisfying results have been obtained, due to a presumably too heterogeneous dataset, where large glaciers (> 10 km 2 ) and/or south-facing glaciers are largely under-represented.Larger-scale studies and multi-variable correlations in between morpho-topographic features could be for instance envisaged.Rabatel et al. (2017) recently proposed an alternative approach to reconstruct the annual mass balance of unmonitored glacier on the basis of the albedo method.This approach relies on the ELA method (Rabatel et al., 2005), but using the remotely sensed monitored α min a together with the AAR, the glacier hypsometry, and the regional SMB elevation gradient (which is the annual SMB gradient in the vicinity of the glacier ELA).For an exhaustive description of this approach, see Rabatel et al. (2017).
Using the albedo method for the summer period has shown promising results, with significant correlations found for the six seasonally monitored glaciers.There is still in this approach a step to retrieve the summer SMB of an unmonitored glacier with high confidence.
The winter period has also been considered in the framework of this study, but has not been presented in the main body of this publication because of underwhelming results.The albedo signal between 1 October and 30 April has been computed similarly to Sirguey et al. (2016) by integrating the (3) According to Sirguey et al. (2016), the use of α T allows for detection of all snowfall events on the glacier by monitoring abrupt rises of α.One of the main conclusions of the above study was the ability of the computed α int w to monitor the frequency of snowfall events, themselves proxy of the accumulation of snow on the glacier, known to be one of the main components of the winter SMB.
α T has been chosen to maximize the correlation between the retrieved cumulative winter albedo α int w and the winter SMB.Threshold values have been computed independently for each of the six seasonally monitored glaciers.To evaluate the impact of this threshold, α int w has also been computed without threshold over winter months (equivalent to α T = 0).Table 3 gathers all the coefficients obtained from the relationship α int w vs. b w , with and without the use of an albedo threshold α T .
For Argentière and Mer de Glace glaciers, a significant correlation is found whatever the value of the albedo threshold α T .For the four other glaciers, using α T largely improves the correlation.However, α T is far from being uniform on the six glaciers (0.53 ≥ α T ≥ 0.76).In addition, for most of the considered glaciers, correlation coefficients abruptly deteriorate when changing this threshold, which does not allow us to use a "regional" threshold for all considered glaciers.On the other hand, Argentière and Mer de Glace without the use of α T provide the best correlation coefficients compared to the other four glaciers; it is noteworthy that they are by far the largest glaciers of our monitoring set (14.59 and 23.45 km 2 for Argentière and Mer de Glace glaciers, respectively).With a glacier snout reaching 1600 m a.s.l., the tongue of these glaciers can experience melting events (resulting in contrasted pixels in terms of albedo value), even during the winter season.Another difference between our study and Sirguey et al. (2016) is that their work focused only on Brewster Glacier, defined as a maritime glacier.These types of glaciers, even during the accumulation period, can experience strongly varying albedos in their lower reaches, which leads to similar behaviours in winter as for Argentière and Mer de Glace glaciers.We therefore reconsider the idea of Sirguey et al. (2016) to use a threshold as a representative value of fresh snowfall, as there is no physical reason that this threshold varies, at least within the same region.However, an interesting perspective would be to apply the method without threshold, on a set of other maritime or large glaciers (> 10 km 2 ).
An additional approach has been carried out, aiming at retrieving b w by deduction from the reconstructed b a and b s from the albedo signal.This approach, not using the winter albedo signal, is poorly correlated (r 2 < 0.16) with in situ b w for the six seasonally monitored glaciers.Indeed, the result extremely depends on the quality of the correlations between b a , b s and the albedo signals.Saint-Sorlin Glacier is a good example, being one of the glaciers with the highest correlations for the annual (r 2 = 0.86) and summer (r 2 = 0.94) SMB.Subtracting b s from b a to computed b w leads to an average difference between computed and measured b w of ±0.41 m w.e. for the 10 simulated years.As a consequence, in the case of low correlations between SMB and albedo, errors in the computed winter SMB become exacerbated.

Conclusion
In this study, we used the so-called albedo method to correlate annual and summer SMB to glacier-wide average albedos obtained from MODIS images.This method has been applied to 30 glaciers located in the French Alps, over the period 2000-2015.Image processing has been performed using the MODImLab algorithm, and filters on the images have been applied, removing images with more than 30 % cloud coverage, and excluding images with satellite observation angles greater than 30 • .Quality assessment has been performed and close agreement has been found between albedos from AWS installed on Saint-Sorlin Glacier and MODIS retrieved albedo values.Annual SMB has been significantly correlated with the summer minimum albedo for 27 of the 30 selected glaciers, confirming this variable as a good proxy of the glacier-wide annual SMB.For the six seasonally monitored glaciers, summer SMBs obtained from the glaciological method have been significantly linked to the integral of the summer albedo.However, calculating the integral of the winter albedo to quantify the winter SMB as done by Sirguey et al. (2016) has shown underwhelming results.Monitoring winter glacier surface albedo may provide good insight into the frequency of snow accumulation at the surface of the glacier but is poor at quantifying the amount of accumulation.Glaciers that experience complete snow coverage during most of the winter season showed the lowest correlation (r 2 ≤ 0.33) while the two glaciers showing the best correlations are subject to some events of surface melting in their www.the-cryosphere.net/12/271/2018/The Cryosphere, 12, 271-286, 2018 lower reaches.This approach should not be definitively forsaken, but it requires improvements in order to confidently retrieve winter SMB.
Sensitivity study on the impact of the considered cloud coverage has revealed a high confidence in the MODIm-Lab cloud algorithm, limiting pixel misclassifications, and a rather high tolerance of the integrated signal to the number of partly cloud-covered images.This confidence on cloud filters is very promising to document unmonitored glaciers.Correction of shadows by the MODImLab algorithm has however revealed some limitations when a large share of the glacier is shadowed by the surrounding topography (around winter solstice).Despite this, severe and artificial drops of albedo in winter have not been identified as an obstacle for monitoring both summer and winter SMB.Such systematic errors are not an issue for inter-annual studies, but would be a serious issue on studies focused on albedo values themselves.For future works, the MODIS archive together with albedo maps, cloud and shadow masks processed with MODImLab, together with validation data from AWS, offer a unique dataset to monitor the temporal and spatial evolution of the surface albedo of glaciers at a regional scale.For instance, computing the absorbed solar radiation (Bair et al., 2016) by date and for each glacier would be an appropriate protocol to estimate the impact of a changing glacier surface albedo in terms of snow or ice melt.Quantifying albedo changes and resulting mass losses with such an approach would be of major interest to better understand the potential effects of possibly increasing dust content, glacier orientation or snow grain growth on glacier surface melt processes.
To conclude, the use of optical satellite images to estimate glacier surface processes and quantify annual and summer SMB from the albedo cycle is very promising and should be expanded to further regions.Using images from different satellites, combining high spatial and temporal resolution instruments, could substantially reduce uncertainties, especially for spotting the albedo summer minimum with more confidence, but also to improve the temporal resolution.This method could then in the short term become reliable for retrieving SMB values of monitored and unmonitored glaciers.
Data availability.Processing of the albedo images has been performed using the open-source MODImLab algorithm.This algorithm can be accessed by contacting its administrator, P. Sirguey.MODIS images have been downloaded using the permanent ftp server: n5eil01u.ecs.nsidc.orgmaintained by the National Aeronautics and Space Administration (NASA) Distributed Active Archive Center (DAAC).
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Map of the region of interest with the studied glaciers shown in red (numbers refer to Table1).The four AWS used in the present study were set up on Saint-Sorlin Glacier (no.20).Adapted fromRabatel et al. (2016).

Figure 2 .
Figure2.Schematic of a typical albedo cycle over one summer, displaying parameters which have been linked to annual and summer (between 1 May and 30 September in the Northern Hemisphere) SMB.α int s is retrieved using Eq.(1).The summer minimum value of albedo is represented by α min a .

Figure 4 .
Figure 4.The ∼ 16-year albedo course for Saint-Sorlin Glacier.Glacier-wide averaged albedo is represented with the continuous black line.The green dots spot for each summer the minimum average albedo, and have been manually checked for all years and glaciers.Dashed red and blue lines stand for the beginning of the defined ablation and accumulation seasons (May and 1 October respectively).

Figure 5 .Figure 6 .
Figure 5. Albedo cycle for Argentière Glacier as a function of the SZA.Each point corresponds to glacier-wide averaged albedo for each available image.The 16 years are displayed.Colour scale gives indication on the date of the used image.The thick grey line describes the weekly albedo averaged over the entire study period.For readability purpose, the averaged albedo has been smoothed, using a seven-point running average.

Figure 7 .
Figure7.Summer SMB b s expressed as a function of the integrated albedo over the entire ablation season for Saint-Sorlin Glacier.Error bars result from the uncertainties related to the glaciological method (measurements and interpolation at the glacier scale of the punctual measurements, ±0.20 m w.e. in total), and on the quadratic sum of the systematic errors made on each albedo measurement.The thin dashed grey line represents the linear regression showing the best correlation between the two variables, together with correlation coefficients.

Figure 8 .
Figure8.The r 2 for the six seasonally surveyed glaciers for the albedo summer integral vs. summer SMB relationship against the cloud threshold above which images have been discarded during the summer season.For the computation, 100 thresholds have been tested between 0 and 100 %.The inner histogram illustrates the number of considered images per summer and averaged on the six glaciers.

Figure 9 .
Figure9.Impact of the ratio of corrected pixels toward the difference between non-corrected and glacier-wide albedo for Argentière Glacier.Each point corresponds to one acquisition and the 16 years are therefore displayed on this graph.Colour scale gives some indication of the date of the acquired image.Grey shaded areas correspond to ratios of corrected pixels for which α non-cor − α has low statistical robustness (refer to the main text).Thin grey lines represent 1σ SD of α, averaged by classes of 5 % corrected pixels.The inner graph illustrates the amount of corrected pixel, as function of the selected month.

Table 2 .
Filtering the images from OZA values.

Table 3 .
Coefficients of determination for the relationship between the winter SMB b w and the integrated winter albedo, computed with and without the albedo threshold α T .