Seasonal changes in surface albedo of Himalayan glaciers from MODIS data and links with the annual mass balance

Few glaciological field data are available on glaciers in the Hindu Kush–Karakoram–Himalayan (HKH) region, and remote sensing data are thus critical for glacier studies in this region. The main objectives of this study are to document, using satellite images, the seasonal changes of surface albedo for two Himalayan glaciers, Chhota Shigri Glacier (Himachal Pradesh, India) and Mera Glacier (Everest region, Nepal), and to reconstruct the annual mass balance of these glaciers based on the albedo data. Albedo is retrieved from Moderate Resolution Imaging Spectroradiometer (MODIS) images, and evaluated using ground based measurements. At both sites, we find high coefficients of determination between annual minimum albedo averaged over the glacier (AMAAG) and glacier-wide annual mass balance (Ba) measured with the glaciological method (R 2 = 0.75). At Chhota Shigri Glacier, the relation between AMAAG found at the end of the ablation season and Ba suggests that AMAAG can be used as a proxy for the maximum snow line altitude or equilibrium line altitude (ELA) on winteraccumulation-type glaciers in the Himalayas. However, for the summer-accumulation-type Mera Glacier, our approach relied on the hypothesis that ELA information is preserved during the monsoon. At Mera Glacier, cloud obscuration and snow accumulation limits the detection of albedo during the monsoon, but snow redistribution and sublimation in the post-monsoon period allows for the calculation of AMAAG. Reconstructed Ba at Chhota Shigri Glacier agrees with mass balances previously reconstructed using a positive degreeday method. Reconstructed Ba at Mera Glacier is affected by heavy cloud cover during the monsoon, which systematically limited our ability to observe AMAAG at the end of the melting period. In addition, the relation between AMAAG and Ba is constrained over a shorter time period for Mera Glacier (6 years) than for Chhota Shigri Glacier (11 years). Thus the mass balance reconstruction is less robust for Mera Glacier than for Chhota Shigri Glacier. However our method shows promising results and may be used to reconstruct the annual mass balance of glaciers with contrasted seasonal cycles in the western part of the HKH mountain range since the early 2000s when MODIS images became available.


Introduction
Approximately 800 million people rely directly on water originating from the high mountains of Asia for fresh water supply and hydropower (Immerzeel et al., 2010).The runoff of the large Asian river systems such as Indus, Ganges or Brahmaputra rivers is partly controlled by the melting of the glaciers located in the Hindu Kush-Karakoram-Himalayan Published by Copernicus Publications on behalf of the European Geosciences Union.
F. Brun et al.: Surface albedo of Himalayan glaciers and links with mass balance (HKH) region (e.g., Immerzeel et al., 2013).These glaciers represent the largest glacierized area in the lower latitudes and are significant contributors to current and future sea level rise (Gardner et al., 2013;Radić et al., 2014).It is therefore important to assess how these glaciers would respond to climate change, and what would be the consequences of their evolution in terms of glacial hazards and modification of local and regional hydrology (e.g., Richardson and Reynolds, 2000;Kaser et al., 2010).
The most popular method to measure glacier mass balance with remote sensing data is the geodetic method.The method is based on the difference between two or more digital elevation models (DEMs), or repeat laser altimetry measurements, acquired at different dates (e.g., Berthier et al., 2007;Kääb et al., 2012;Nuimura et al., 2012).The geodetic method allows a region-wide glacier mass balance to be calculated, and has revealed the heterogeneous response of HKH glaciers to climate change (Gardelle et al., 2013;Kääb et al., 2012).However, the geodetic method provides estimates of height and volume change averaged over typical periods of 5 years or more (e.g., Bolch et al., 2011;Kääb et al., 2012;Nuimura et al., 2012;Gardelle et al., 2013), and therefore fails to capture the inter-annual or seasonal variability of glacier mass balance.Furthermore, although this method is accurate over large regions, it is less efficient when applied to a single glacier for which no high resolution topographic data are available (e.g., Vincent et al., 2013).
The equilibrium-line altitude (ELA) and minimum mean bi-hemispherical broadband albedo (referred to hereafter as albedo) of a glacier have the potential to be good proxies of the annual mass balance in some regions (e.g., Dumont et al., 2012;Rabatel et al., 2013;Shea et al., 2013).On temperate glaciers, the average albedo of a whole glacier surface reaches a minimum at the end of the ablation season, when the elevation of the transient snow line elevation reaches a maximum (Rabatel et al., 2005;Dumont et al., 2012).The maximum elevation of the transient snow line is often referred to as a proxy of the ELA (Rabatel et al., 2005).The ELA, and thus the annual minimum albedo averaged on the whole glacier (AMAAG), is strongly correlated with the annual glacier mass balance (B a ; Dumont et al., 2012).Here we present the first test of this method on two Himalayan glaciers: Chhota Shigri Glacier and Mera Glacier.
The goal of this study is to examine the links between albedo and B a at two Himalayan glaciers.Specifically, our objectives are the following: (1) to use ground measurements of albedo to evaluate the accuracy of albedo retrievals algorithms for Moderate Resolution Imaging Spectroradiometer (MODIS) data, (2) to assess the relation between albedo and B a , (3) to characterize the seasonal variability of the albedo; and (4) to use albedo records from MODIS to reconstruct the mass balance of both glaciers since 2000.

Chhota Shigri Glacier
Chhota Shigri Glacier (32.3 • N, 77.6 • E) is a valley-type glacier located in the Chandrabhaga River basin, Pir Panjal Range, India (Fig. 1).This glacier extends from 6263 to 4050 m a.s.l. and covers a total area of 15.7 km 2 .It is mostly free of debris except in the lower ablation area (≤ 4500 m a.s.l.).Approximately 3.4 % of the total glacierized area is debris-covered (Vincent et al., 2013).A central moraine separates the glacier into two branches below 4800 m a.s.l., and the glacier is fed by numerous tributaries exhibiting various aspects (Fig. 2).Its glacier-wide annual mass balance has been monitored since 2002 using the ground-based glaciological method (Azam et al., 2012(Azam et al., , 2014b)).The equilibrium-line altitude for a zero annual mass balance (ELA 0 ) was found to be close to 4900 m a.s.l.(Wagnon et al., 2007).

Mera Glacier
Mera Glacier (27.7 • N, 86.9 • E) is a 5.1 km 2 debris free glacier situated in the Dudh Koshi basin, Everest region, Nepal (Fig. 1).The main glacier flows from 6420 m a.s.l. and divides into two branches around 5800 m a.s.l.The main branch (referred hereafter as Mera branch, located on the western part of the glacier) flows down to 4940 m a.s.l., while the smaller eastern branch (Naulek branch) flows down to 5260 m a.s.l.The eastern part of the Naulek branch is also part of the Naulek glacierized complex that flows down from Naulek Peak (Fig. 3).Glacier-wide mass balance has been monitored at Mera Glacier since 2007 using the glaciological method, and ELA 0 is close to 5550 m a.s.l.(Wagnon et al., 2013).

Climate setting
The southern flank of the Himalayas is dominated by two contrasting climate regimes: the westerlies and the Indian summer monsoon.The influence of the Indian monsoon decreases from east to west and dominates the climatic signal for longitudes eastern from approximately 77 • E. Conversely, the influence of the westerlies decreases from west to east (Bookhagen andBurbank, 2006, 2010;Yao et al., 2012;Maussion et al., 2013).Monsoon precipitation occurs from June to September as moist air masses from the Bay of Bengal sweep toward the Himalayas.These orographic rainfalls are highly related to the local relief (Bookhagen and Burbank, 2006), with a strong south-north gradient due to the presence of the Himalayan range.At high altitudes, monsoon-season precipitation is characterized by high frequency and low intensity (Shrestha et al., 2012).Westerly circulation patterns bring precipitation mainly during the winter months (November to April).They are produced by synopticscale upper-tropospheric waves amplified by the topography with a capacity to generate heavy snowfalls (Hatwar et al., 2005).
Because of both climatic regimes, the seasonality of rainfall varies greatly with the location in the range (Yao et al., 2012).According to Burbank et al. (2012) and local precipitation measurements (e.g., Wagnon et al., 2013), precipitation amount at Mera Glacier is mostly governed by the Indian monsoon, with 70 to 80 % of precipitation occurring during summer time (Wagnon et al., 2013).Chhota Shigri Glacier is influenced by both monsoon and westerly circulation systems, and receives precipitation throughout the year (Azam et al., 2014a, b).These observations are consistent with the glaciological data: Mera Glacier is a summer-accumulation type glacier (Maussion et al., 2013;Wagnon et al., 2013), while Chhota Shigri Glacier is a winter-accumulation-type glacier despite being located in the transition zone and influenced by both regimes (Wagnon et al., 2007;Maussion et al., 2013).
The accuracy of the albedo retrieval from satellite data was evaluated using shortwave radiation observations from automatic weather stations (AWSs) installed in the ablation zones of the two glaciers (Figs. 2 and 3) at 4670 and 5360 m a.s.l. for Chhota Shigri and Mera glaciers, respectively.Incident and reflected shortwave radiation (SW in and SW out , respectively) were measured with Kipp & Zonen CM1 sensors (0.3 ≤ λ ≤ 2.8 µm) every 30 s and averaged every half-hour.According to the manufacturer the sensors have an intrinsic accuracy of ±3 %.In turn, the albedo, calculated as the ratio between SW out and SW in , has a compounded accuracy of about 4 %.Nevertheless, the local slope, the tilt of the de- vice, and frost deposition on the sensor can increase the uncertainty.At Mera Glacier the AWS collected data from 9 November 2009 to 10 November 2010 and from 27 November 2012 to 16 November 2013.At Chhota Shigri Glacier the AWS collected data from 13 August 2012 to 2 February 2013 and from 8 July 2013 to 3 October 2013 (Fig. 4).
Snow and ice albedo are highly sensitive to cloud cover, which affects the ratio of direct and diffuse radiation (e.g., Warren, 1982;Gardner and Sharp, 2010).Therefore, we filtered the data to retain only those collected in clear-sky conditions to be able to compare them to satellite-retrieved albedo.To do so, we calculated the local clear-sky broadband atmospheric transmissivity χ as the ratio between the value defining the upper envelope of measured SW in at ground level and that defining the upper envelope of SW in at the top of the atmosphere (TOA) (Garnier and Ohmura, 1968).We found clear-sky χ equal to 0.99 and 0.85 on Mera and Chhota Shigri glaciers, respectively.Clear-sky days were selected as the days for which χ was in the range 0.99 ± 4 % for Mera Glacier and 0.85 ± 10 % for Chhota Shigri Glacier at local solar noon.The lower estimated atmospheric transmissivity at Chhota Shigri Glacier is likely due to its lower elevation, its sky-view factor, and the vicinity of important aerosol sources.The very high clear-sky χ estimated above Mera Glacier is not realistic.It is probably due to reflections on surrounding snowy slopes and/or tilting of the sensor which artificially increase the incoming shortwave radiation measured by the AWS.Snow and ice albedo also de- pend strongly on the solar zenith angle (Warren, 1982).In order to allow a direct comparison between field and satellite measurements, we averaged the AWS albedo over a 2 h time span centered on the MODIS satellite passing time, which is coincident with a high zenithal angle.

MOD10A1 products
Daily MODIS/Terra Snow Cover Daily L3 Global 500 m SIN Grid (MOD10A1, referred hereafter as MOD10), processed and distributed by the National Snow and Ice Data Center, are available at http://nsidc.org(Hall et al., 2006).The products provide an estimate of daily blue-sky (BS) albedo, which corresponds to the broadband albedo for the actual direct and diffuse illumination, daily snow cover, and a quality assessment (QA) at 500 m resolution (Klein and Stroeve, 2002).The daily albedo is calculated from the daily MODIS reflectances corrected from the atmospheric effects (MODIS/Terra Surface Reflectance Daily L2G Global 500m SIN Grid; MOD09GHK).The adjustment of snow and ice anisotropic scattering effects is done using a DIScrete Ordinates Radiative Transfer model (DISORT).Through additional processing of MOD09 products, we filtered the albedo products to keep only those calculated from images acquired with a sensor zenith angle lower than 40 and Mera (bottom) glaciers.The values reported are the daily maxima of the radiative flux.The daily albedo is calculated as the average of the albedo on a 2 h time span centered on the satellite passing time.Note the difference in x axis scale.The vertical grey lines show the days detected as cloudy, using a threshold on incoming irradiance (see text for details).

MCD43 products
MODIS MCD43A3 products provide white-sky spectral albedo (WS albedo, which corresponds to the albedo of the surface illuminated only by diffuse irradiance), black-sky spectral albedo (which corresponds to the albedo of the surface illuminated only by direct irradiance), and broadband albedo (called WS shortwave albedo), for the seven visible and near infrared bands of MODIS at local solar noon.MCD43A3 products are produced at 500 m resolution every 8 days using 16 days of MODIS Terra and Aqua acquisitions.The separate products MCD43A2 are used to provide quality assessment for the MCD43A3 products.These data are available through the online data pool at the NASA Land Processes Distributed Active Archive Center (https: //lpdaac.usgs.gov/).

MODImLab products
MODImLab is an algorithm developed to provide snow and ice albedo with a 250 m resolution (Sirguey et al., 2009;Dumont et al., 2012).It retrieves albedo in mountainous regions from MODIS TOA calibrated radiances corrected with topographic and atmospheric data (see Table 1 for a complete list of inputs and their availability).MODImLab supports both BS albedo and WS albedo, as well as a cloud mask derived from MODIS reflective and emissive bands.We filtered the albedo products to keep only those calculated from images acquired with a sensor zenith angle lower than 40 • .For a complete description of the albedo retrieval method, the reader is referred to Sirguey et al. (2009), Sirguey (2009), Dumont et al. (2011), andDumont et al. (2012).

Glacier masks
In order to compute glacier albedo from georeferenced satellite imagery provided at 250 and 500 m resolutions, we created raster masks of the glaciers with the same resolutions following Dumont et al. (2012).Glacier masks of 250 m resolution were derived from glacier boundaries used in Wagnon et al. (2007) and Wagnon et al. (2013), and the 500 m raster masks derive from the 250 m masks.These masks were manually defined using SPOT5 (2.5 m resolution) and Pléiades (1 m resolution) orthoimages acquired on 21 September 2005 and 25 November 2012, respectively.They slightly differ from the glacier outlines of the Randolph Glacier Inventory (Pfeffer et al., 2014) but are believed to be more accurate because they have been determined based on visual inspection of high resolution images and field expertise.
For the raster glacier mask at Chhota Shigri Glacier (Fig. 2), we removed all pixels corresponding to debriscovered zones (all pixels below 4500 m a.s.l.corresponding to the main tongue), to rocky areas (southeastern tributaries), to narrow valleys (eastern tributary) and to steep slopes with seracs (northwestern tributary), which might affect the satellite-derived albedo.The resulting mask contains 89 pixels.As this mask is large compared to other studies (Dumont et al., 2012 used a 29 pixel mask), a mean albedo can be assessed for the glacier even when it is partially cloud covered (applied coverage threshold is 20 %).In contrast, for Mera Glacier we defined two separate masks to conduct a parallel analysis.Mask 1 (red line in Fig. 3) includes the whole glacier (44 pixels).Conversely, we retained only the 16 pixels corresponding to elevations lower than the ELA 0 (5550 m a.s.l.) for mask 2 (black line in Fig. 3).We proceeded in this way because Mera Glacier is almost always covered by clouds during the second half of the ablation sea- son (end of July to September).Further, for such a summeraccumulation type glacier, the satellite images obtained immediately after the monsoon already exhibit snow accumulation at higher elevations, hiding the glacier surface representative of the end of ablation season.Moreover, calculating the albedo on a small part of the glacier (for mask 2) increases the availability of good quality images, especially when the higher part of the glacier is covered by clouds, while the lower part remains cloud-free, which is often observed in spring.

Comparison with field measurements
We compared the broadband albedo measured under clear sky conditions at the AWSs with the BS albedo retrieved by MODImLab and MOD10 algorithm on the corresponding pixel for both glaciers (Fig. 5).While ground-based and satellite-retrieved measurements of albedo are not directly comparable, as their spatial and temporal resolutions differ, they will be coherent in time.Indeed, the albedo measured from the AWS comes from a 80 m 2 surface, whereas the albedo derived from satellite has been calculated over a 62 500 or 250 000 m 2 surface.Additionally, the AWS albedo is averaged over a 2 h time span centered on the satellite passing time.As the MCD43 products provide an albedo averaged over 16 days we did not compare it against instantaneous field measurements.The limited number of albedo samples (n = 15 for MODImLab and n = 16 for MOD10) and the relatively small variance (σ = 0.05 for the AWS albedo used to compare with MOD10) on Chhota Shigri Glacier prevented a reliable conclusion about the quality of the albedo retrieval techniques (Fig. 5a).As the MOD10 albedo estimates correspond to snow-covered periods at the Chhota Shigri AWS, the regression is poorly constrained.At Mera Glacier the measurement period is long enough to thoroughly assess the quality of the albedo retrieval (Fig. 5b).The agreement between albedo retrieval from MOD10 and AWS measurements is slightly better than between MODImLab and AWS measurements (RMSE are 0.08 and 0.10, respectively).The mean difference between the AWS measurements and MODImLab retrieval is 0.03 vs. 0.15 for MOD10.Thus the MOD10 albedo retrieval appears to be biased towards lower values at the Mera site.This is in contradiction with other studies that found a bias toward higher values in MOD10 data (e.g., Şorman et al., 2007).The mean absolute errors are 0.10 for MODImLab and 0.16 for MOD10 data.
In conclusion, MODImLab and MOD10 retrievals seem to have the same quality in the context of this study.Nevertheless, given the higher spatial resolution of MODImLab and the ability of its cloud detection algorithm to better detect clouds, we choose to mostly rely on MODImLab products for our analyses.

Comparison of cloud masks
The performances of the cloud masks calculated by MOD-ImLab and that provided by the MOD35 Cloud product (and used in MOD10 products) were compared through a visual inspection of 13 ASTER images with variable cloud coverage (Table 2) covering the Mera Glacier area and corresponding to a total of 13 312 and 3328 pixels for the 250 and 500 m resolutions, respectively (Fig. 6).MODIS and ASTER images were acquired at the same time (no more than 5 minutes apart), and were compared using four categories, taking the ASTER mask determined by visual inspection as a reference: - Reference data were obtained via visual interpretation of clouds in simultaneous ASTER images.Category 1 corresponds to cloudy pixels detected as cloudy pixels, category 2 to cloudy pixels detected as cloud-free pixels, category 3 to cloudfree pixels detected as cloudy pixels and category 4 to cloud-free pixels detected as cloud-free pixels.Categories 1 and 4 correspond to the correctly classified pixels (green background) and categories 2 and 3 correspond to the misclassified pixels (red background).
spectively.The MODImLab cloud mask correctly detected 96 % of cloud-obscured pixels, whereas MOD10 detected only 67 % (percentages calculated from the frequencies of categories 1 and 2 on Fig. 6).MODImLab misclassified 46 % of the cloud-free pixels, whereas MOD10 misclassified only 21 % of them (percentages calculated from the frequencies of categories 3 and 4 on Fig. 6).Overall, the MODImLab cloud detection algorithm is more conservative than the MOD35 algorithm, as cloud coverage is overestimated.To avoid unrealistic values of albedo retrieved from MODIS data, we relied on the MODImLab cloud mask at the cost of reducing the number of available images assumed to be free of clouds.

Multi-annual trend of the albedo signal
Following Ming et al. (2012) we calculated the 2000-2013 trend in MOD10 albedo for the highest elevation pixels of both glaciers (located at 5191 m a.s.l. for Chhota Shigri Glacier and 6008 m a.s.l. for Mera Glacier).The linear trends are +4.3×10−4 ±7×10 −5 and −3.6×10 −3 ±3×10 −5 yr −1 , respectively (Fig. 7).As the MODImLab cloud mask is more conservative than the MOD35 cloud mask (i.e., the cloud cover is overestimated by MODImLab), the albedo values obtained with MODImLab are not uniformly distributed along the year, and therefore introduce a sampling bias in the trend calculation.Consequently, we did not use MODImLab to compute a multi-annual trend of the albedo signal.
We also compared the MOD10 albedo averaged over 16 days to MCD43 albedo (produced from 16 days of acquisition).MCD43 quality flags guaranteed only 4 % of high quality retrieval (with full bidirectional reflectance function inversion) on Mera Glacier and 21 % on Chhota Shigri Glacier.Such quality assurance was considered too low to use MCD43 albedo products for trend analysis and for the rest of the study.

Seasonal variations of albedo
We linearly interpolated between the dates with information the mean albedo values at a daily time step to calculate an average albedo cycle for each glacier (Fig. 8).Although such a linear interpolation may sometimes hide unobserved changes  ± 7 × 10 −5 yr −1 for Chhota Shigri Glacier and −3.6 × 10 −3 ± 3 × 10 −5 yr −1 for Mera Glacier).The red dots are the MCD43 8 day albedo constructed from multiangle acquisitions over 16 days.Only 22 % of MCD43 data points correspond to good quality products on Chhota Shigri Glacier and only 4 % on Mera Glacier.Note that many data points acquired during monsoon time above Mera Glacier, or during winter time above Chhota Shigri Glacier, probably correspond to undetected cloudy pixels. in albedo, we believe that the seasonal cycle is well captured on Chhota Shigri Glacier where observations are numerous.
On the other hand, for Mera Glacier, the cycle presented in Fig 8 has to be considered with caution.On Chhota Shigri Glacier, the annual albedo signal exhibits a marked seasonality.Maximum albedos are observed in winter (November to February), and decrease through the spring and sum-mer to a minimum in August (Fig. 8 and Table 3).In contrast, the cycle observed for the lower part of Mera Glacier (mask 2) exhibits a less pronounced seasonality compared with Chhota Shigri due in part to cloud cover during the monsoon.The glacier surface is seldom apparent from space in visible and near infrared channels during summer time, thus potentially concealing some variability of the annual cycle.Consequently, when the glacier can be observed again after the monsoon, it is mostly snow covered, even in the ablation area, and has a relatively high albedo ( 0.7).
For Mera Glacier, the 16-day average curve of MOD10 products shown in Fig. 7 exhibits a stronger seasonality.Given that the detection rate of clouds by MOD35 cloud mask is relatively low (67 %, Fig. 6) and that the cloud coverage during monsoon season is heavy, the albedo minima visible on Fig. 7 are potentially artifacts due to pixel misclassification.Consequently, the albedo annual cycle observed on Fig. 7 cannot be fully considered as a true cycle.

Albedo and annual glacier-wide mass balance
The mean albedo of a glacier has been recognized as being a potentially good proxy for the transient snow line altitude because of the contrast in albedo between snow and ice (Dumont et al., 2012).Among selected Alpine glaciers regularly surveyed, the altitude of the transient snow line reaches a maximum at the end of the ablation season which corresponds to the ELA (Rabatel et al., 2005).In turn, the mean albedo of the glacier surface reaches a minimum and can effectively be used as a proxy to estimate the ELA.
We processed daily MODIS data from May 2000 to December 2013 for the two glaciers.Daily maps of WS albedo were calculated for the two glaciers and averaged over cloudfree glacierized pixels for images with a cloud coverage over the glacier lower than 20 % for Chhota Shigri Glacier and 0 % for Mera Glacier (Figs. 2 and 3).WS albedo does not depend on solar zenith angle, and is thus a more relevant variable than BS albedo to compare albedo maps generated at different dates, and therefore with different solar illumination.Variations in WS albedo can be interpreted more readily as a variation in characteristics of the surface.Following Dumont et al. (2012) we extracted the minimum value of the albedo signal for the period corresponding to the end of the ablation season of Chhota Shigri Glacier.From field observations and precipitation regime, the snow line was expected to be at its highest elevation between late July and late September.Consequently, we processed images systematically from July to December and only partially for the rest of the year to reduce the amount of L1B data to be downloaded and processed.Minimum albedo values range from 0.22 to 0.46, and the AMAAG correlates well with the fieldbased measurements of B a (Fig. 9a).This is supported by a determination coefficient (R 2 ) equal to 0.75 (n = 11 years; significant at 95 % confidence with a Student's t test).The agreement is better for negative mass balance years.All the dates of the observed minimum albedo are between 25 July and 13 September (Table 3).It is also noteworthy that the AMAAG values are scattered when B a is 0 or slightly positive for Chhota Shigri Glacier.
On Mera Glacier the snow line was expected to be the highest at the end of the monsoon season (between August and September) during which most of the ablation and accumulation occur (Wagnon et al., 2013).Unfortunately, the minimum albedo value during these months could hardly be observed from space because of the persistent cloud coverage during this period.On the first images captured after the monsoon, the glacier is usually covered by snow that likely conceals the ELA information.Nevertheless, Wagnon et al. (2013) observed that between October and April, the snow deposited above Mera La (5400 m a.s.l., Fig. 3) is systematically re-mobilized or sublimated by strong winds.As a consequence, the surface representative of the glacier at the end of the previous monsoon season, unfortunately not visible from space at this period, is likely to be exposed again in the www.the-cryosphere.net/9/341/2015/The Cryosphere, 9, 341-355, 2015 w.e. and ±0.28 m.w.e. for Chhota Shigri and Mera glaciers, respectively.The uncertainty on the mean albedo is calculated as the accuracy of the MODImLab albedo established on Mera Glacier (RMSE = 0.10) divided by the root mean square of the number of pixels in the average (n = 89, n = 44 and n = 16 for Chhota Shigri and Mera glaciers masks 1 and 2, respectively).At Mera Glacier, the 2012-2013 mass balance has been calculated excluding the heavy snow fall that occurred 13-15 October 2013 to be consistent with the end of the hydrological year ending with the ablation season.It is noteworthy that this mass balance also has a larger uncertainty (not shown on the figure) because of ablation measurements were partly compromised by stakes buried in the snow.Note that the scales are different on the two figures.
course of the following post-monsoon or winter seasons.We thus hypothesized that the annual ELA information could be retrieved after the monsoon season during the rest of the hydrological year, before melting starts again.Albedo minima for Mera Glacier were thus searched for between August and June of the following year.
When albedo is averaged over the entire glacier (mask 1), no significant relation is found between AMAAG and B a (Fig. 9b).However, when averaged over the ablation area (mask 2), the minimum albedo is highly correlated with B a (determination coefficient of 0.75; n = 6 years; significant at 95 % confidence with a Student's t test, Fig. 9c).Con-sequently for the rest of the study, we will use mask 2 to average albedo over the ablation area only, and therefore AMAAG refers to "annual minimum albedo averaged over the glacier ablation zone".However, the occurrences of these minima are scattered over the entire year, from 10 August to 27 May of the following hydrological year (Table 3), while the albedo values are in a higher and narrower range than observed on Chhota Shigri Glacier (0.47 to 0.55).2013) and the annual mean of the reconstructed mass balance, respectively.The purple line corresponds to the annual mass balance values calculated with a degree-day model (Azam et al., 2014a).The error bars for the reconstructed mass balance correspond to the 90 % confidence interval.For Mera Glacier some reconstructed mass balance values (for the years 2001 and 2007) correspond to albedo values outside the range of calibration and are therefore believed to be questionable.Note that the scaling of mass balance is different on both panels.

Mass balance reconstruction, 1999-2013
Using the significant relations obtained between albedo and glacier-wide mass balance, we reconstructed B a for years with available MODIS images.For each year, the AMAAG was retrieved between July and September on Chhota Shigri Glacier, and between August and June the next year on Mera Glacier.These minima were then inverted to estimate B a , with the linear regression parameters of mass balance expressed as a function of AMAAG (Fig. 10).This method is referred to hereafter as the albedo method.
We calculated the uncertainty associated with these reconstructed values by performing two Monte Carlo simulations.First, we did a simulation to obtain a distribution of the regression parameters assuming that the errors on both mass balance and albedo were normally distributed (Fig. 11).Then, using this distribution, for every albedo minimum, we randomly chose 100 000 slope-intercept couples and calculated the corresponding B a .We kept the values included between the 5th and 95th percentiles of this distribution to obtain the 90 % confidence intervals for the reconstructed mass balance values.
Both the observed and reconstructed B a values were averaged over their corresponding period of record to allow comparison with the glacier-wide geodetic mass balance reported in Gardelle et al. (2013) for the period 1999-2011.Although a systematic bias may exist between the geodetic mass balance and B a obtained with the glaciological method used to calibrate the albedo method (e.g., Zemp et al., 2013), this comparison allows us to check the order of magnitude of our results.The average reconstructed mass balance on Chhota Shigri Glacier was found to be −0.68 ± 0.10 m.w.e.yr −1 , which is significantly more negative than the value derived from the geodetic method (−0.37 ± 0.19 m.w.e.yr −1 ).Nonetheless, the reconstructed B a values for the period October 1999-October 2002 are in good agreement with the values reconstructed using a degree-day model (Azam et al., 2014a).For Mera Glacier the reconstructed mass balance is on average positive (0.13 ± 0.19 m.w.e.yr −1 ) as was the geodetic mass balance estimate from Gardelle et al. (2013) (0.18 ± 0.31 m.w.e.yr −1 ).However, the uncertainties on mass balance values retrieved with the albedo method are large, and some retrieved mass balance values are estimated from albedo values out-of-calibration range (years 2000-2001 and 2006-2007).This is likely to lead to a questionable mass balance average.

Limitations of the optical remote sensing in the study area
In this study we first assessed the quality of remotely sensed albedo retrieval techniques using ground-based AWS data, which is rarely possible in the HKH region due to the difficulty in glacier access.Nevertheless, we stress that a direct comparison between satellite products and contemporaneous AWS data is complicated by the inconsistent footwww.the-cryosphere.net/9/341/2015/The Cryosphere, 9, 341-355, 2015 The relatively high rate of misclassification of cloudy pixels (33 %) in MOD10 products raises questions about the significance of the calculated decadal trend.In addition to the sampling problem associated with the persistent cloud coverage during winter time on Chhota Shigri Glacier and monsoon season on Mera Glacier, the observed decadal trends might not be related to changes in the glacier surface, but could also be related to changes occurring in the atmosphere or due to a drift in MODIS sensors calibration (e.g., Zhang and Reid, 2010).Nevertheless these trends are consistent with those observed by Ming et al. (2012), at glaciers PT1 and PT4, which are located relatively close to Chhota Shigri and Mera glaciers.Albedo trends at PT1 and PT4, respectively, are +6.4 × 10 −3 and −1.1 × 10 −3 yr −1 .Albedo trends at Chhota Shigri Glacier and Mera Glacier are +4.3×10−4 ±7×10 −5 and −3.6×10 −3 ±3×10 −5 yr −1 , respectively.These trends are both significant for a confidence level of 85 % with a Student's t test.Ming et al. (2012) observed only negative albedo trends except for the most western glacier of their study PT1.The positive trend we obtained at Chhota Shigri Glacier, close to PT1, seems to confirm this finding.
Special cloud detection algorithms are required in mountainous and high altitude regions (Hall and Riggs, 2007;Sirguey et al., 2009).In particular, thresholds on the different bands, band ratios, and band combinations need adjustments when applied to different regions.For instance, the relatively low water vapor content of the atmosphere above Mera Glacier (clear sky χ of 0.99) yields frequent misclassifications of cloud free pixels as high altitude clouds in MOD-ImLab cloud detection algorithm.Our work suggests that independent tuning of cloud detection algorithms is required, and this requires extensive comparison with other interpreted images, a process that is both time consuming and operator dependent.

Advantages, drawbacks, and robustness of the albedo method applied in the Himalayas
An attractive application of the albedo method is to reconstruct B a for periods when surface observations are not available.The mass balance reconstruction is robust for Chhota Shigri Glacier but appeared questionable for Mera Glacier.We suggest three main reasons for this difference.First, the choice of pixels included in Mera Glacier mask is questionable.In our case selecting only the lower part of the glacier lead to a strong improvement in the correlation between B a and AMAAG (Fig. 9b and c).This approach is only suitable for glaciers where ELA 0 is known.As mask 2 includes only pixels located below the ELA 0 (Fig. 3), the mean albedo cannot be related to ELA detection in the cases when ELA is above ELA 0 (i.e., negative B a observed for years 2008-2009, 2009-2010 and 2011-2012).
Second, persistent cloud coverage at Mera Glacier does not allow the true AMAAG to be resolved in some years.In order to test this hypothesis we also considered the second lowest value of the mean albedo observed between August and June on Mera Glacier.This revealed that the second lowest albedo value is often observed more than two months apart from the AMAAG and often exceeds the actual observed minimum by 0.05 to 0.3.In contrast, the second minimum mean albedo measured on Chhota Shigri Glacier is close to the minimum (+0.01) for most years and is generally observed within 10 days.Moreover, the visual inspection of the albedo time series (Fig. 8) confirms a well-defined minimum on Chhota Shigri Glacier, whereas the seasonal cycle of the mean albedo on Mera Glacier is poorly resolved and less marked.
Third, the relation between mass balance and AMAAG is poorly constrained for summer accumulation-type glaciers, where accumulation and ablation occur simultaneously.For Mera Glacier, the range of albedo values used for calibration is narrower than that on Chhota Shigri Glacier (0.47 to 0.55 vs. 0.21 to 0.47).It is also noteworthy that the albedo is calculated from fewer pixels on Mera Glacier than on Chhota Shigri Glacier (16 and 89, respectively), and so the mean albedo obtained for Mera Glacier has a higher uncertainty (Fig. 9b).The AMAAG-B a relation is also calibrated over a shorter period for Mera Glacier than for Chhota Shigri Glacier (6 years vs. 11).Therefore the regression parameters are less constrained on Mera Glacier than on Chhota Shigri Glacier (Fig. 11), which leads to higher uncertainty.
A future application of the albedo method could be to calculate mass balance of neighboring glaciers.This would require more studies to examine regression parameter from the relation between B a and AMAAG in the context of glacier topography and climate.These studies should preferentially be conducted on glaciers that exhibit a strong albedo seasonal cycle, like Chhota Shigri Glacier.Ming et al. (2012) presented the seasonal albedo cycle of 11 Himalayan glaciers, including 7 of them with a well-defined albedo cycle that could be suitable candidates for such a study.Two of their studied glaciers (Naimona'nyi and East Rongbuk glaciers) exhibit an albedo cycle with a phase opposition compared to the others (i.e., the albedo reaches a maximum during summer and a minimum during winter), and further research into the application of the albedo method is required for summeraccumulation type glaciers.

Conclusions
In this study we presented and assessed a method to reconstruct the annual mass balance of Himalayan glaciers from MODIS images.This method could be applied to other recently surveyed glaciers and would allow B a to be estimated from the hydrological year 2000 to present.The relation between B a and the AMAAG must be calibrated over several years with mass balance values obtained from other methods, such as the glaciological method.
Despite the intrinsic limitations of moderate resolution imagery (subpixel variability, exclusion of pixels at the edge of the glacier, presence of ice and rocks mixed pixels, applicability only on large glaciers), our study has successfully applied the MODImLab albedo retrieval algorithm to two Himalayan sites.MODImLab is therefore believed to be a suitable tool to monitor glacier surface albedo in the Himalayas.
Our albedo method for mass balance reconstruction was successful for Chhota Shigri Glacier (a winter-accumulationtype glacier), but questionable in the case of Mera Glacier (a summer-accumulation type glacier).While our albedoderived estimates of glacier mass balance are validated with geodetic mass balance observations, time series of ground-based mass balance data with sufficient length are still required to construct robust regression models.Spatial extrapolation of these results to climatologically similar regions and validation with geodetic mass balance estimates would provide a better indication of the transferability of our results, and provide important information on the inter-annual variability glacier mass balance in unmonitored regions of the Hindu Kush Himalayas.

Figure 1 .
Figure 1.Map of the study area.The altitude ranges given in the legend are in meters.The glacierized areas (from the Randolph Glacier Inventory version 3.2, Pfeffer et al. (2014)), are in blue and the two studied glaciers are represented by red dots.

Figure 2 .
Figure 2. Example of an albedo map retrieved by MODImLab over Chhota Shigri Glacier on 7 September 2001 05:50:00 (UTC).The bold line delimits the pixels used to average the albedo (n = 89, corresponding to 4.875 km 2 ).The northern lower part of the glacier is covered with debris and the tributaries exhibit surfaces that include many rock outcrops.The blue dot shows the location of the AWS.

Figure 3 .
Figure 3. Example of an albedo map retrieved by MODImLab over Mera Glacier on 11 August 2006 04:15:00 (UTC).The bold red line delimits the mask 1 (n = 44 pixels, corresponding to 2.75 km 2 ).The bold black line delimits the mask 2 (n = 16 pixels, corresponding to 1 km 2 ) with pixels located below 5550 m a.s.l.The blue dot shows the location of the AWS.

Figure 5 .Figure 6 .
Figure 5. Measured albedo vs. MODIS-retrieved albedo for Chhota Shigri Glacier (a) and Mera Glacier (b).For Mera Glacier RMSE are 0.10 and 0.08 for the MODImLab and the MOD10, respectively.Please note that the number of MOD10 and MODImLab images is different because the two data set are filtered with different cloud detection algorithms.The inset represents the distribution of the measured albedo (grey = satellite, black = AWS).

Figure 7 .
Figure7.MOD10 daily albedo (light blue dots) and 10-day running average (solid blue line) for the highest pixel of Chhota Shigri Glacier (upper panel; pixel at 5191 m a.s.l.) and Mera Glacier (lower panel; pixel at 6008 m a.s.l.).The solid black line represents the 13 year linear trend in albedo (+4.3 × 10 −4 ± 7 × 10 −5 yr −1 for Chhota Shigri Glacier and −3.6 × 10 −3 ± 3 × 10 −5 yr −1 for Mera Glacier).The red dots are the MCD43 8 day albedo constructed from multiangle acquisitions over 16 days.Only 22 % of MCD43 data points correspond to good quality products on Chhota Shigri Glacier and only 4 % on Mera Glacier.Note that many data points acquired during monsoon time above Mera Glacier, or during winter time above Chhota Shigri Glacier, probably correspond to undetected cloudy pixels.

Figure 8 .
Figure 8. Seasonal variations of the mean albedo for Chhota Shigri (top) and Mera with mask 2 (bottom) glaciers.The solid grey line represents the average of the interpolated daily mean albedo with a 1-σ envelope (dashed grey lines).Note that the scaling of albedo is different on both panels.The size of dots is proportional to the number of good quality images available for the Julian day.The vertical red bars indicate the dates of annual minima for each year under consideration.

Figure 9 .
Figure9.AMAAG as a function of B a for Chhota Shigri Glacier (a) and Mera Glacier with mask 1 (b) and mask 2 (c).The error bars for the mass balance are ±0.40 m.w.e. and ±0.28 m.w.e. for Chhota Shigri and Mera glaciers, respectively.The uncertainty on the mean albedo is calculated as the accuracy of the MODImLab albedo established on Mera Glacier (RMSE = 0.10) divided by the root mean square of the number of pixels in the average (n = 89, n = 44 and n = 16 for Chhota Shigri and Mera glaciers masks 1 and 2, respectively).At Mera Glacier, the 2012-2013 mass balance has been calculated excluding the heavy snow fall that occurred 13-15 October 2013 to be consistent with the end of the hydrological year ending with the ablation season.It is noteworthy that this mass balance also has a larger uncertainty (not shown on the figure) because of ablation measurements were partly compromised by stakes buried in the snow.Note that the scales are different on the two figures.

Figure 10 .
Figure 10.Reconstructed mass balance for Chhota Shigri (top) and Mera (bottom) glaciers.The green and yellow lines represent the annual mean of the geodetic mass balance from Gardelle et al. (2013) and the annual mean of the reconstructed mass balance, respectively.The purple line corresponds to the annual mass balance values calculated with a degree-day model(Azam et al., 2014a).The error bars for the reconstructed mass balance correspond to the 90 % confidence interval.For Mera Glacier some reconstructed mass balance values (for the years 2001 and 2007) correspond to albedo values outside the range of calibration and are therefore believed to be questionable.Note that the scaling of mass balance is different on both panels.

Figure 11 .
Figure 11.Regression parameters obtained from the Monte Carlo simulation (N = 10 000).Every data point represents a slopeintercept couple obtained as the best fit of B a as a function of AMAAG.

Table 2 .
Acquisition dates of the ASTER images used to assess the quality of the different cloud masks.

Table 3 .
Values of the observed minima of the mean albedo and dates at which they are observed.At Mera Glacier, minimum albedo dates can occur between 1 August and 30 June of the following year.The data used for the calibration are in bold.