1. INTRODUCTION
Glacier mass balance (MB) is highly sensitive to atmospheric conditions and constitutes a direct and perceptible indicator of climate change (Haeberli and Beniston, Reference Haeberli and Beniston1998). On a global scale, the mass loss of glaciers (apart from the Greenland and Antarctica ice sheets) constitutes the major contribution to present sea-level rise (IPCC, 2013). For example, over the well-observed 2003–2009 period, glacier mass loss was 259 ± 28 Gt a−1, representing about 30 ± 13% of total sea-level rise (Gardner and others, Reference Gardner2013). Between 1902 and 2005, glacier mass loss, reconstructed by Marzeion and others (Reference Marzeion, Leclercq, Cogley and Jarosch2015), corresponds to 63.2 ± 7.9 mm sea-level equivalent (SLE). Glacier mass losses are projected to increase in the future and are expected to range between 155 ± 41 mm SLE (for the emission scenario RCP4.5) and 216 ± 44 mm SLE (for RCP8.5) between 2006 and 2100 (Radić and others, Reference Radić2014). Moreover, on continental surfaces, glacier behaviour strongly determines watershed hydrology: glacier meltwater impacts on the runoff regime of mountainous drainage basins and consequently on the availability of the water resource for populations living downstream (Huss, Reference Huss2011). Many scientific, economic and societal stakes are thus related to glaciers (Arnell, Reference Arnell2004; Kaser and others, Reference Kaser, Grosshauser and Marzeion2010), hence justifying regular monitoring.
Traditionally, the annual and seasonal glacier MBs are measured using glaciological in situ measurements from repeated snow accumulation surveys and ice ablation stake readings on individual glaciers, that are subsequently extrapolated over the entire glacier area (Østrem and Brugman, Reference Østrem and Brugman1991; Kaser and others, Reference Kaser, Fountain and Jansson2003). Nevertheless, this laborious technique requires heavy logistics, is time consuming and moreover is subject to potential systematic errors (e.g. Thibert and others, Reference Thibert, Blanc, Vincent and Eckert2008). Glaciological MB measurements are thus limited to a few easily accessible glaciers, unequally distributed around the globe: only about 150 glaciers among the world's ~200 000 glaciers are regularly monitored in the field and present at least ten consecutive years of MB measurements (Pfeffer and others, Reference Pfeffer2014; Zemp and others, Reference Zemp2015).
Several remote sensing techniques are used to increase data coverage of glacier mass change at a global scale. One of them is the geodetic method, which consists of monitoring glacier elevation changes by differencing multitemporal DEMs. These DEMs can be derived from aerial photos and/or airborne Light Detection and Ranging (LiDAR) (Soruco and others, Reference Soruco2009; Abermann and others, Reference Abermann, Fischer, Lambrecht and Geist2010; Jóhannesson and others, Reference Jóhannesson2013; Zemp and others, Reference Zemp2013). However, the main limitation of airborne sensors is their geographic coverage, restricted to areas accessible with airplanes; some large remote areas such as High Mountain Asia can hardly be monitored. LiDAR data from the Ice Cloud and Elevation Satellite (ICESat) altimeter can also provide information on glacier elevation changes, but these data are too sparse to allow reliable monitoring of the mass budget of individual glaciers (Kropáček and others, Reference Kropáček, Neckel and Bauder2014; Kääb and others, Reference Kääb, Treichler, Nuth and Berthier2015). DEMs can also be derived from medium-to-high resolution optical satellite images such as Satellite pour l'Observation de la Terre-High Resolution Stereoscopic (SPOT5/HRS) images or radar satellite images as done during the Shuttle Radar Topographic Mission (SRTM). However, the vertical precision of these DEMs (~5 m) is such that it remains challenging to accurately estimate the MB of individual small to medium-sized glaciers (glaciers < 10 km2) (Paul and others, Reference Paul, Frey and Le Bris2011; Gardelle and others, Reference Gardelle, Berthier, Arnaud and Kääb2013). More recently, Worldview and Pléiades sub-meter stereo images have been used to provide accurate glacier topography with a vertical precision of ±1 m and even ±0.5 m on gently sloping glacier tongues (Berthier and others, Reference Berthier2014; Willis and others, Reference Willis, Herried, Bevis and Bell2015). Nevertheless, repeatedly covering all glaciers on Earth with Worldview and Pléiades stereo images would be very expensive given their limited footprint for a single scene (e.g. 20 km × 20 km for Pléiades). Due to DEM errors and uncertainties in the density value for the conversion from volume change to mass change (Huss, Reference Huss2013), this method currently provides MB at a temporal resolution typically, of 5–10 a. This does not provide an understanding of glacier response to climatic variations at seasonal and annual timescales (Ohmura, Reference Ohmura2011).
An alternative method to estimate the MB change of a glacier is based on the fluctuations of the glacier equilibrium line altitude (ELA) (Braithwaite, Reference Braithwaite1984; Kuhn, Reference Kuhn and Oerlemans1989; Leonard and Fountain, Reference Leonard and Fountain2003; Rabatel and others, Reference Rabatel, Dedieu and Vincent2005). The ELA is the altitude on a glacier that separates its accumulation zone, where the annual MB is positive, from its ablation zone, where the annual MB is negative (Ohmura, Reference Ohmura2011; Rabatel and others, Reference Rabatel2012). The ELA and annual glacier-wide MB are strongly correlated. Early studies have shown that the ELA can be efficiently approximated by the snowline altitude at the end of the ablation season (i.e. at the end of the hydrological year) on mid-latitude glaciers (LaChapelle, Reference LaChapelle1962; Lliboutry, Reference Lliboutry1965). The snowline altitude measured at the end of the ablation season thus constitutes a proxy for estimating the annual MB. This method has been applied to glaciers located in the Himalayas, Alaska, Western North America, Patagonia and the European Alps amongst others, using aerial photographs or optical remote sensing images from Landsat and MODIS (Ostrem, Reference Ostrem1973; Mernild and others, Reference Mernild2013; Rabatel and others, Reference Rabatel, Letréguilly, Dedieu and Eckert2013; Shea and others, Reference Shea, Menounos, Moore and Tennant2013). However, the Landsat temporal resolution of 16 d can be a strong limitation for accurately studying the dynamics of snow depletion, especially in mountainous regions with widespread cloud coverage (Shea and others, Reference Shea, Menounos, Moore and Tennant2013); and the methodologies used with MODIS (available since 2000) to detect ELA often require complex algorithms implying heavy classification (Pelto, Reference Pelto2011; Shea and others, Reference Shea, Menounos, Dan Moore and Tennant2012; Rabatel and others, Reference Rabatel, Letréguilly, Dedieu and Eckert2013).
Unlike the above-mentioned techniques, only a few recent studies have focused on retrieving seasonal glacier MBs (Hulth and others, Reference Hulth, Rolstad Denby and Hock2013; Huss and others, Reference Huss2013). Seasons are in fact, more relevant timescales to better understand the glacier response to climate fluctuations, allowing the partitioning of recent glacier mass loss between changes in accumulation (occurring mainly in winter) and changes in ablation (occurring mainly in summer), at least in the mid- and polar latitudes (Ohmura, Reference Ohmura2011). Hulth and others (Reference Hulth, Rolstad Denby and Hock2013) combined melt modelling using meteorological data and snowline tracking measured in the field with GPS to determine winter glacier snow accumulation along snowlines. Huss and others (Reference Huss2013) calculated the evolution of glacier-wide MB throughout the ablation period also by combining simple accumulation and melt modelling with the fraction of snow-covered surface mapped from repeated oblique photography. Both methods have been applied for only one or two glaciers, and a few years as they are time-consuming and require various datasets (meteorological data, remotely sensed data) that are not widely available.
In this context, the consistent archive of optical SPOT/VEGETATION (SPOT/VGT) images, spanning 1998–2014 with a daily and nearly global coverage at 1 km resolution (Sylvander and others, Reference Sylvander, Henry, Bastien-Thiry, Meunier and Fuster2000; Maisongrande and others, Reference Maisongrande, Duchemin and Dedieu2004; Deronde and others, Reference Deronde2014), remains an unexploited dataset to estimate the year-to-year seasonal MB variations of glaciers. Instead of detecting the exact physical snowline altitude of a glacier as done in earlier studies, our new method focuses on monitoring a ‘mean regional’ seasonal altitude of snow (Z), in a region surrounding a glacier, during summer (from early May to the end of September) and winter (from early October to the end of April). The interannual dynamics of this altitude Z, estimated from SPOT/VGT images, is then compared with the interannual variation of observed seasonal MBs (derived from direct glaciological measurements), available continuously every year over 1998–2008, for 55 glaciers in the Alps (WGMS, 2008, 2012, 2013; Huss and others, Reference Huss, Hock, Bauder and Funk2010a, Reference Huss, Usselmann, Farinotti and Bauderb, Reference Huss, Dhulst and Bauder2015). The performance of the 55 individual linear regressions between Z and MB are analysed in terms of the linear determination coefficient R 2 and the RMSE. For each of the 55 regressions, we then perform a cross-validation of the temporal robustness and skill of the regression coefficients. For all 55 glaciers, seasonal MBs are then calculated outside of the calibration period, i.e. each year of the 2009–2014 interval. A validation is realized for some glaciers with observed seasonal MB over 2009–2014. Finally, we discuss the possible causes of differences in the results between seasons, between glaciers and present the limits and perspectives of our methodology.
2. DATA
2.1. Satellite data
This study uses optical SPOT/VGT satellite images provided by two instruments: VEGETATION 1 from April 1998 to the end of January 2003 (VGT 1 aboard the SPOT 4 satellite launched in March 1998), and VEGETATION 2 from February 2003 to the end of May 2014 (VGT 2 aboard the SPOT 5 satellite launched in May 2002). The instruments provide a long-term dataset (16 a) with accurate calibration and positioning, continuity and temporal consistency. The SPOT/VGT sensors have four spectral bands: the blue B0 (0.43–0.47 µm), the red B2 (0.61–0.68 µm), the Near Infra-Red B3 (0.78–089 µm) and the Short Wave Infra-Red SWIR (1.58–1.75 µm). Both sensors provide daily images of almost the entire global land surface with a spatial resolution of 1 km, in a ‘plate carree’ projection and a WGS84 datum (Deronde and others, Reference Deronde2014). In this study, we used the SPOT/VGT-S10 products (10-daily syntheses) freely available between April 1998 and May 2014 (http://www.vito-eodata.be/PDF/portal/Application.html#Home). The SPOT/VGT-S10 products, processed by Vlaamse Instelling voor Technologisch Onderzoek (VITO), result from the merging of data strips from ten consecutive days (three 10-d syntheses are made per month) through a classical Maximum Value Composite (MVC) criterion (Tarpley and others, Reference Tarpley, Schneider and Money1984; Holben, Reference Holben1986). The MVC technique consists of a pixel-by-pixel comparison between the 10-daily Normalized Difference Vegetation Index (NDVI) images. For each pixel, the maximum NDVI value at top-of-atmosphere is picked. This selection allows minimizing cloud cover and aerosol contamination. Due to atmospheric effects, an atmospheric correction is applied to the VGT-S10 products (Duchemin and Maisongrande, Reference Duchemin and Maisongrande2002; Maisongrande and others, Reference Maisongrande, Duchemin and Dedieu2004) based on a modified version of the SMAC (Simple Method for the Atmospheric effects) code (Rahman and Dedieu, Reference Rahman and Dedieu1994). The SPOT/VGT-S10 Blue, Red and SWIR spectral bands have been gathered from 1 April 1998 to 31 May 2014, over a region including the Alps and stretching from 43 to 48.5°N and from 4 to 17°E (Fig. 1).
Fig. 1. Map of the mean NDSI for the period 1998–2014, over the European Alps, and location of the studied glaciers with MB balance measurements (sources 1–3; Table S1). (a) Entire study area. (b) Enlargement of the studied area (red rectangle in (a)).
The terrain elevation over the studied region has been extracted from the SRTM30 DEM, derived from a combination of data from the SRTM DEM (acquired in February 2000) and the U.S. Geological Survey's GTOPO30 dataset (Farr and Kobrick, Reference Farr and Kobrick2000; Werner, Reference Werner2001; Rabus and others, Reference Rabus, Eineder, Roth and Bamler2003). The SRTM30 spatial resolution is 30 arcsec (~ 900 m) and the geodetic reference is the WGS84 EGM96 geoid. The DEM has been resampled using the nearest-neighbour method to a 1 km spatial resolution in order to match the resolution of the S10 SPOT/VGT images.
2.2. Glacier MB data
Glacier-wide winter and summer MB data of 55 different glaciers are used in this study (Table S1 in the Supplementary Material; Fig. 1), covering a total area of ~400 km2, corresponding to 19% of the 2063 km2 glacier area in the Alps (Pfeffer and others, Reference Pfeffer2014). Winter MB (B w) is measured through to the end of April. Summer MB (B s) is the difference between annual MB determined at the end of the ablation period (end of September) and winter MB (Thibert and others, Reference Thibert, Eckert and Vincent2013). Winter MB corresponds to the snow/ice mass accumulated during the winter season (from Octoberyear−1 to Aprilyear of measurement); thereafter, for more clarity, B w_obs refers to the winter MB measured in April of the year of measurement and B s_obs to the summer MB measured in September of the year of measurement.
The different MB datasets all start before 1998. However, as we later compare these MB data with SPOT/VGT snow maps available since April 1998, we used seasonal MB only after 1998.
The first MB dataset is composed of direct glaciological measurements provided by the World Glacier Monitoring Service (WGMS, 2013). Seasonal MB of 44 Austrian, French, Italian and Swiss glaciers of the European Alps are available between 1950 and 2012 (for the longest series) but with discontinuous data. At the time of access to the database, continuous seasonal MB time series were only available for seven glaciers since 1998 (labelled ‘source 1’ in Table S1).
The second MB dataset provided by Huss and others (Reference Huss, Dhulst and Bauder2015) is derived from recently re-evaluated glaciological measurements. For seven Swiss glaciers (labelled ‘source 2’ in Table S1), long-term continuous seasonal MB series from 1966 (at least) to 2014 are inferred from point measurements extrapolated to the entire glacier using distributed modelling, with year-to-year variability directly given by the in situ measurements.
The third MB dataset, composed of 41 Swiss glaciers, is a comprehensive set of field data homogenized using distributed modelling, with year-to-year variation constrained by meteorological data (Huss and others, Reference Huss, Hock, Bauder and Funk2010a, Reference Huss, Usselmann, Farinotti and Bauderb). The seasonal MB time series of 21 glaciers covering 30% of the total glacier area of Switzerland (labelled ‘source 3a’ in Table S1), are available for each year of the 1908–2008 period (Huss and others, Reference Huss, Hock, Bauder and Funk2010a). For 20 other glaciers located in the southeastern Swiss Alps (labelled ‘source 3b’ in Table S1), seasonal MB time series are provided each year from 1900 to 2008 (Huss and others, Reference Huss, Usselmann, Farinotti and Bauder2010b). The continuous seasonal MB series of these 41 glaciers are derived from three different data sources:
-
(1) Geodetic data of volume change for periods of four to ~ 40 a; three to nine DEMs have been used for each glacier since 1930, mostly originating from aerial photogrammetry or terrestrial topographic surveys for the first DEM.
-
(2) (Partly isolated) in situ measurements of accumulation and ablation for about half of the glaciers; they were used to constrain MB gradients and winter accumulation.
-
(3) Distributed accumulation and temperature-index melt modelling (Hock, Reference Hock1999; Huss and others, Reference Huss, Bauder, Funk and Hock2008) forced with meteorological data of daily mean air temperature and precipitation, recorded at various weather stations close to each of the glaciers. The model is constrained by ice volume changes obtained from DEM differencing and calibrated with in situ measurements when available. Huss and others (Reference Huss, Bauder, Funk and Hock2008) provide more details on the model and calibration procedure.
In order to have the most complete set of continuous interannual MB data to calibrate our method, we chose to combine MB time series from all sources (1, 2, 3a, b) over the same period. Since 41 glaciers (annotated 3a, b in Table S1) provide continuous seasonal MB until 2008, the calibration period extends until 2008. As we later compare MB data with SPOT/VGT images provided since April 1998, 1999–2008 constitutes the winter calibration period and 1998–2008 the summer calibration period. In total, 55 glaciers with 10 a of winter MB (B w_obs) over 1999–2008, and 11 a of summer MB (B s_obs) over 1998–2008 have thus been used to calibrate our approach. Coordinates, median elevation and glacier area are available for all 55 glaciers (Table S1).
Furthermore, seasonal MBs are also available for some of these 55 glaciers outside the calibration period (2009–2014, i.e. the end of the period covered by SPOT/VGT data). In winter, 66 additional B w _ obs data are available in total for 19 glaciers over 2009–2014 and in summer, 49 B s_obs measurements are available for 13 glaciers over 2009–2013 (summer 2014 is not covered by SPOT/VGT) (WGMS, 2013; Huss and others, Reference Huss, Dhulst and Bauder2015).
Errors associated with B w _ obs and B s_obs are not provided individually per glacier. Zemp and others (Reference Zemp2013) find an average annual random error of 340 mm w.e. a−1 for 12 glaciers with long-term MB measurements. Huss and others (Reference Huss, Hock, Bauder and Funk2010a) have quantified the uncertainty in glacier-wide winter balance (see Supplementary Information in their paper) as ±250 mm w.e. a−1. This number is not directly applicable to all glaciers (depending on the sampling), but to most of them. Systematic errors associated with in situ glaciological MBs are expected to be within 90–280 mm w.e. a−1 if cautious measurements are realized with sufficient stake density (Dyurgerov and Meier, Reference Dyurgerov and Meier2002), but they are generally assumed to range between 200 and 400 mm w.e. a−1 for one glacier (Braithwaite and others, Reference Braithwaite, Konzelmann, Marty and Olesen1998; Cogley and Adams, Reference Cogley and Adams1998; Cox and March, Reference Cox and March2004). In this study, we consider an error (noted E obs) of ±200–400 mm w.e. a−1 for observed MB.
3. METHODS
3.1. Cloud filtering and temporal interpolation of the snow cover maps
Snow discrimination is based on its spectral signature, which is characterized by a high reflectance in the visible and a high absorption in the SWIR wavelength. Consequently, the Normalized Difference Snow Index (NDSI) first introduced by Crane and Anderson (Reference Crane and Anderson1984) and Dozier (Reference Dozier1989) for the Landsat sensor, defined as a band ratio combining the visible (green) and SWIR Landsat bands, constitutes a good and efficient proxy to map the snow cover with optical remote sensing. The NDSI has since been widely used with different sensors (e.g. Fortin and others, Reference Fortin, Bernier, El Battay, Gauthier and Turcotte2001; Hall and others, Reference Hall, Riggs, Salomonson, DiGirolamo and Bayr2002; Salomonson and Appel, Reference Salomonson and Appel2004). The NDSI value is proportional to the snow cover rate of the pixel and allows a monitoring of the spatial and temporal snow cover variations (Chaponnière and others, Reference Chaponnière2005). In this study, we build an NDSI adapted to the SPOT/VGT sensor (without green channel), inspired by Chaponnière and others (Reference Chaponnière2005). This modified NDSI is computed from the mean of the red B0 and blue B2 channels (to recreate an artificial green band) and from the SWIR:
$$NDSI = \displaystyle{{(B0 + B2)/2 - SWIR} \over {(B0 + B2)/2 + SWIR}}.$$
A different (and more standard) formulation of the NDSI is used for the MODIS sensor and calculated only with the red B0 channel and the SWIR (Xiao and others, Reference Xiao, Shen and Qin2001). Our entire methodology was also tested with this standard NDSI but led to results of inferior quality (not shown).
Clouds have a spectral signature similar to snow and they can be misclassified as snow, especially at the edge of the snow pack. We apply a cloud mask proposed by Dubertret, (Reference Dubertret2012) (and refer to as D-12 cloud mask below) in order to flag cloudy pixels and to avoid overestimating snow coverage. This algorithm, based on various reflectance bands threshold tests, has been adapted to SPOT/VGT images from different cloud masks initially created for higher resolution imagery. It corresponds to the crossing of Sirguey's cloud mask (developed for MODIS sensor) (Sirguey and others, Reference Sirguey, Mathieu and Arnaud2009) and of Irish and Zhu's cloud masks, both developed for the Landsat sensor (Irish, Reference Irish2000; Zhu and Woodcock, Reference Zhu and Woodcock2012). When compared with other classical cloud detection algorithms (e.g. Berthelot, Reference Berthelot2004), the D-12 cloud mask performs the best cloud identification, along with Lissens and others (Reference Lissens, Kempencers, Fierens and Van Rensbergen2000) cloud mask. However, Dubertret (Reference Dubertret2012) concluded that the D-12 cloud mask is more conservative than that of Lissens and others (Reference Lissens, Kempencers, Fierens and Van Rensbergen2000) and allows detection of more-snow covered pixels (by flagging less clouds). We therefore choose to apply the cross-cloud mask on the S10 syntheses, which are composites of the S1 daily images pixels over 10 d.
A temporal interpolation is then computed for the ‘cloudy’ pixels when possible. If a pixel is detected as cloudy in a S10 synthesis, its value is replaced by the mean of the same pixel value in the previous synthesis (S10 t−1), and in the next one (S10 t+1) if these pixels are not cloudy. If they are, we compute the mean of the S10 t−2 and the S10 t+2 synthesis pixel values. In order to produce maps of winter/summer NDSI, we average all 10-d NDSI syntheses included between 1 October and 30 April for each winter of the 1999–2014 period, and between 1 of May and 30 of September for each summer of the 1998–2013 period.
3.2. Calibration
We first superimpose each interannual mean seasonal NDSI map derived from SPOT/VGT for the period 1998–2014 on the SRTM30 DEM. Then, considering different sized of square windows varying from 5 to 401 km side lengths (with steps of 2 km) of P × P pixels and centred on each glacier, we derive the altitudinal distribution of the mean seasonal NDSI. For each glacier, we thus obtain the interannual dynamics of the NDSI altitudinal distribution within different Windows of Snow Monitoring (WOSM) sizes surrounding the glacier (16 curves; Fig. 2). The WOSM side lengths tested are always odd numbers for the glacier to be situated in the central pixel of the square window. Then for each WOSM size, from the intersection between a NDSI value (varying from 0.2 to 0.65, with a step of 0.01) and each seasonal curve of the NDSI altitudinal distribution, a ‘mean regional’ altitude of snow (Z) can be deduced each year, from 1998 to 2014. We do not focus on detecting the exact physical snowline elevation of the glacier at the end of the melt season to estimate annual MB, as done in previous studies (e.g. Rabatel and others, Reference Rabatel, Dedieu and Vincent2005). In fact, the SPOT/VGT 1-km resolution is not adapted to monitor the snowline elevation and its high spatial variability. Moreover, the snowline approach does not allow us to retrieve seasonal MB, because at the end of the accumulation period (end of April/beginning of May in the Alps), the entire glacier is generally covered by snow. For these reasons, we aim at estimating for winter and for summer a statistical mean regional altitude of snow in a region surrounding a glacier. Finally, for each WOSM size and each NDSI value that are tested, we compute a linear regression between the mean regional snow altitudes Z and the observed MBs (Fig. 3), over the calibration period (1999–2008 for winter, 1998–2008 for summer). This linear regression allows an interannual estimation of the glacier's seasonal MB, as a function of Z inferred from SPOT/VGT, as written in the linear equation (generalized for both seasons):
$$B_{{\rm w/s\_VGT}} = \alpha _{{\rm w/s}} \times Z_{{\rm w/s}}\; + \; \beta _{{\rm w/s}}.$$
B w/s_VGT is the winter/summer MB estimated for the year y. Z w/s is the winter/summer ‘mean regional’ altitude of snow; the slope coefficient α w/s expressed in mm w.e. a−1 m−1 represents the sensitivity of a glacier winter/summer MB towards Z w/s and β w/s is the winter/summer intercept term expressed in mm w.e. a−1.
Fig. 2. Altitudinal distribution of NDSI for each year since 1998. The NDSI has been averaged in a square window centred on Griesgletscher, central Swiss Alps. The red horizontal line represents the NDSI value from which the mean regional snow altitude Z (represented by the red vertical line) is inferred for each year. (a) Winter NDSI over 1999–2014 (WOSM size: 367 × 367 km²; NDSI value: 0.43). (b) Summer NDSI over 1998–2014 (WOSM size: 117 × 117 km²; NDSI value: 0.54).
Fig. 3. Observed (a) winter and (b) summer MB of Griesgletscher, central Swiss Alps, as a function of the mean regional snow altitude Z for each year of the calibration period represented by coloured dots. Dashed thin lines represent the 95% confidence intervals for linear regression (solid line).
Coefficients of determination R 2 w/s and RMSE w/s_cal are computed to assess the quality of the regression over the calibration period for each season.
To summarize, for each glacier, a linear regression is computed for each plausible NDSI value and all WOSM sizes. Following some initial tests on a subset of the 55 glaciers, the interval of plausible NDSI values is set to [0.2–0.65]. This NDSI range has also been chosen to include the reference NDSI value of 0.4 commonly accepted to classify a pixel as snow-covered (and associated with a pixel snow cover rate of 50%) (Hall and others, Reference Hall, Riggs and Salomonson1995, Reference Hall, Foster, Verbyla, Klein and Bensont1998; Salomonson and Appel, Reference Salomonson and Appel2004; Hall and Riggs, Reference Hall and Riggs2007; Sirguey and others, Reference Sirguey, Mathieu and Arnaud2009). Then, for each glacier, and each size of WOSM tested, the seasonal NDSI value optimizing the RMSE w/s_cal is selected. This NDSI value adjustment for each individual site (and each season) allows a better adaptation to the glacier-specific environment (e.g. local land cover type, local topography, etc.).
After adjusting the NDSI value, the size of the P × P pixels WOSM surrounding each glacier is also adjusted for each glacier and each season. Here again, the WOSM size minimizing the mean RMSE w/s_cal is selected (Table S1 in Supplementary Material). The WOSM side lengths tested are always odd numbers for the glacier to be situated in the central pixel of the window. A mandatory condition for the WOSM size selection is the continuity of the NDSI altitudinal curves, given the 100 m altitudinal step. The cost function f optimizing the RMSE w/s_cal for each glacier is thus:
$$\eqalign{& f\,(NDSI{^\ast},WOSM{^\ast},\alpha _{{\rm w/s}},\beta _{{\rm w/s}}) = \cr & \quad \min \left( {\sqrt {\displaystyle{{\mathop \sum \nolimits_{y = 1}^N {(B_{{\rm w/s\_VGT}\;{\rm y}} - B_{{\rm w/s\_ref\;} \;{\rm y}})}^2} \over N}}} \right).} $$
NDSI* and WOSM* are respectively the NDSI value and the WOSM side length minimizing the RMSE w/s_cal. Allowing the adjustment of both the NDSI value and the WOSM size is a means to select an optimized quantity of snow-covered pixels (where snow dynamics occurs) that are less affected by artefacts such as residual clouds, aerosols and/or directional effects. Changing NDSI value is a means to scan the terrain in altitude, while changing the windows size is a means to scan the terrain in planimetry.
Consequently, for each glacier, a unique ‘optimized’ linear regression (regarding the parameters α and β) allows us to estimate the seasonal MB from the mean regional snow altitude deduced with SPOT/VGT images.
3.3. Validation
To validate our results, we performed two types of evaluation: (1) cross-validation over the period 1998/1999–2008 and (2) evaluation against recent glaciological MB measurements (not used in the calibration) over the period 2009–2013/2014.
3.3.1. Period 1998/1999–2008
In order to validate the 55 individual optimized linear regressions, we use a classical leave-one-out cross-validation method based on the reconstruction of MB time series where each MB value estimated for the year y is independent of the observed MB for the same year y (Michaelsen, Reference Michaelsen1987; Hofer and others, Reference Hofer, Mölg, Marzeion and Kaser2010; Marzeion and others, Reference Marzeion, Hofer, Jarosch, Kaser and Mölg2012). The cross-validation constitutes an efficient validation mechanism for short time series. For each glacier, we first determine the decorrelation time lag t lag (a), after which the autocorrelation function of the observed MB drops below the 90% significance interval, i.e. for which the serial correlation in the observed MB data is close to zero (for each glacier, t lag = 1). After that, for each glacier with N available observed MBs (N = 10 in winter and N = 11 in summer) we perform N linear regressions between the regional mean snow altitude Z w/s and the observed MB B w/s_obs, leaving each time a moving window of 1 a ± t lag (i.e. 3 a) out of the data used for the regression. The removed value (Z w/s,y and B w/s_obs,y ) for the year y has to be at the centre of the moving window such that the remaining values used for the regression are independent of the removed value. We then obtain N values for the regression coefficients α and β (termed α w/s_cross and β w/s_cross) and N values of reconstructed MB (B w/s_VGT_cross). Standard deviation of the N regression coefficients σ (α w/s_cross) and σ (β w/s_cross) are computed to assess the temporal stability and the robustness of the parameters α w/s and β w/s. The mean regression coefficients α w/s_cross_best and β w/s_cross_best, representing the average of the N α w/s_cross,y and β w/s_cross,y are also calculated. For each glacier, we compute R 2 w/s_cross (as the mean of the N R 2 w/s_ obtained for the N regressions) and an estimate of the error RMSE w/s_cross defined as:
$$RMSE_{{\rm w/s\_cross}} = \left( {\sqrt {\displaystyle{{\mathop \sum \nolimits_{y = 1}^N {(B_{{\rm w/s\_VGT\_cross,}\;y} - B_{{\rm w/s\_obs,}\;y})}^2} \over N}}} \right).$$
We then estimate the skill score of the linear regression as
$$SS_{{\rm w/s}} = 1 - \displaystyle{{RMSE_{{\rm w/s\_cross}}^2} \over {RMSE_{{\rm w/s\_ref\_cross}}^2}}, $$
where RMSE w/s_obs_cross is the mean square error of a reference model. As the reference model, we determine B w/s_ref_cross,y for each year y by averaging the observed MB values leaving out B w/s_obs,y , (allowing that B w/s_ref_cross,y is independent of the observed B w/s_obs,y ). The skill score can be interpreted as a parameter that measures the correlation between reconstructed and observed values, with penalties for bias and under (over) estimation of the variance (Wilks, Reference Wilks2011; Marzeion and others, Reference Marzeion, Hofer, Jarosch, Kaser and Mölg2012). A negative skill score means that the relationship computed has no skill over the reference model to estimate the seasonal MB.
3.3.2. Period 2009–2013/2014
As the period covered by SPOT/VGT data stretches until June 2014, it is possible to calculate the seasonal MB after 2008 outside our calibration period for each glacier. From the intersection between the optimized NDSI value fixed for each glacier and the altitudinal distribution of the NDSI, we deduce a value of seasonal mean regional snow altitude for each year during 2009–2014 in winter and during 2009–2013 in summer. Therefore, with the individual ‘optimized’ relations computed over the calibration period, seasonal MB inferred from SPOT/VGT can be estimated over 2009–2013/2014. In winter, 66 additional B w _ obs data are available for 19 glaciers over 2009–2014 and in summer, 49 B s_obs measurements are available for 13 glaciers over 2009–2013. Therefore, the annual and global RMSE (RMSE w/s_eval) and Mass Balance Error (MBEw/s_eval) can be calculated over each evaluation period.
4. RESULTS
4.1. Performance over the calibration period 1998/1999–2008
We analysed the performance of the linear regression model for the entire glacier dataset during winter and summer. The analysis was first performed individually for each glacier before considering the average results of the 55 glaciers.
4.1.1. Individual glaciers
Figure 4 presents individual performances for the 55 glaciers ranked by increasing RMSE over the calibration period. In winter, correlations between Z w and B w _ obs are high (Fig. 4a): the mean R 2 w for the 55 glaciers is 0.84. R 2 w ranges between 0.31 (for Aletsch, #54 in Table S1) and 0.97 (for Seewjinen #4), with high first quartile (0.79) and median values (0.88). Except for Aletsch, all R 2 w are superior to 0.5. Glaciers from source 3 result in higher mean R 2 w (0.88) than glaciers from source 1 (0.76) and 2 (0.72). Furthermore, the mean RMSE w_cal of the estimated B w_VGT for the 55 glaciers is smaller (158 mm w.e. a−1) than the range of annual errors acknowledged for the glaciological MB measurements (E obs = ±200–400 mm w.e. a−1). RMSE w_cal ranges between 50 mm w.e. a−1 (for Gorner, #1) and 260 mm w.e. a−1 (for Ciardoney, #55), with third quartile and median values of 189 and 159 mm w.e. a−1, respectively. Even the highest RMSE w_cal are in the same range as the lower limit for E obs. Glaciers less well-ranked in terms of RMSE w_cal (Ciardoney, Aletsch, Cantun, respectively #55, #54, #53) are not necessarily the least well-ranked in terms of R 2 w. Cantun for example presents a relatively high RMSE w_cal (237 mm w.e. a−1) but also a high R 2 w of 0.89 because average winter MBs are higher for this glacier. Thus, the two metrics (R 2 and RMSE) are not redundant for characterizing the linear relationships between Z and observed MB for all the glaciers.
Fig. 4. Regression results for the 55 glaciers in terms of R 2, RMSE and seasonal MB. (a) RMSE (black curve) and R 2 (red curve) over the winter calibration period 1999–2008, for the 55 studied glaciers ranked by increasing RMSE. Time series of observed winter MB (b) and winter MB estimated with VGT (c) over 1999–2008. For (b) and (c), each horizontal row of coloured rectangles represents the MB time series of a glacier (the rectangle colour indicating the MB value). Lower panels (d)–(f): respective summer equivalents of graphs (a)–(c), over 1998–2008, for the 55 glaciers ranked by increasing RMSE. The glaciers ranking is not the same for the two seasons (winter and summer glacier ranking in Table S1).
As shown in Figure 4d, the correlation between Z s and B s_obs is also high in summer but not as much as in winter: the mean R 2 s for the 55 glaciers is 0.74, and the difference is also significant for the first quartile and median values (respectively 0.70 and 0.77 in summer instead of 0.79 and 0.89 in winter). R 2 s is between 0.52 (Clariden, #55) and 0.88 (Sarennes, #23). The difference between the minimum and maximum of R 2 s (0.36) is inferior to the same value in winter (0.66), indicating a lower spread of the correlation for summer. The mean RMSE s_cal of the estimated B s_VGT for the 55 glaciers (314 mm w.e. a−1) is twice as large as in winter (158 mm w.e. a−1) but still acceptable compared with the error range of glaciological MB measurements. RMSE s_cal ranges between 134 and 528 mm w.e. a−1 (with third quartile and median values of respectively 355 and 299 mm w.e. a−1). Glaciers from source 2 (Table 1) present higher mean RMSE s_cal and lower R 2 s (391 mm w.e. a−1 and 0.70 respectively) than glaciers from source 1 (325 mm w.e. a−1 and 0.77 respectively) and source 3 (300 mm w.e. a−1 and 0.74 respectively). The summer RMSE s_cal range (394 mm w.e. a−1) is larger than in winter (232 mm w.e. a−1) indicating a wider spread in RMSE s_cal during this season. Moreover, we observe that glaciers with a poor performance in summer do not necessarily present the same poor performance in winter (Table S1).
Table 1. Summary of the linear regression results over the calibration period 1998/1999–2008 according to individual data sources and averaged for the 55 glaciers
Figure 4b, c show B w_obs and B w_VGT time series estimated with SPOT/VGT, over the calibration period 1999–2008, for winter and for the 55 glaciers. By comparing Figure 4b with c, it is seen that for all glaciers the B w _ obs and B w_VGT time series are similar, illustrating that the ‘optimized’ relationships computed from SPOT/VGT allow a good estimation of the interannual variations in winter MB over 1999–2008. We note that Aletsch (the biggest glacier in the European Alps), ranked 54th in winter (and with the lowest R 2), presents a ‘flat’ winter observed MB time series, with low interannual variability.
According to Figures 4e, f, B s_VGT time series are in agreement with B s_obs time series over 1998–2008. As in winter, the mean differences between B s_obs and B s_VGT are the highest for 2 a with extreme summer MBs: 2003 (with strongly negative MB values) and 2007 (a year with above average MB). The very negative and atypical summer MB time series of Sarennes (#23) is well estimated with SPOT/VGT (R 2 s of 0.88 and RMSE s_cal of 282 mm w.e. a−1).
By comparing Figures 6c, d (blue dots), we observe that the regression coefficient α values of all glaciers are much more scattered in summer than in winter. In winter, α w ranges between −2.7 and −0.6 mm w.e. a−1 m−1, whereas α s ranges between −15 and −1 mm w.e. a−1 m−1 in summer. For both seasons, we observe that the lower the absolute α value, the lower the RMSE over the calibration period.
4.1.2. Performance for the average of all observed glaciers
The mean MB time series averaged for all glaciers is also interesting to study as it illustrates glacier behaviour at the scale of the entire European Alps.
In winter (Fig. 5a), 〈B w_VGT〉 globally fits well with 〈B w_obs〉. Absolute mean MB errors |MBE w_cal | (i.e. the mean difference between B w_VGT and B w_obs for all the glaciers) are maximal in 2005 (175 mm w.e. a−1), 2001 (110 mm w.e. a−1), and to a lesser extent, in 2007 (95 mm w.e. a−1) (Table 2). Nevertheless, these errors are small compared with the standard deviation of the observed winter MB of all glaciers σ w_obs computed for 2001, 2005 and 2007, respectively equal to 615, 373 and 272 mm w.e. a−1.We also note that the two contrasting winter balance years in 2000 and 2001 are well captured by the model. In summer (Fig. 5b), 〈B s_VGT 〉 fits less well with 〈B s_obs〉 than in winter although the largest interannual variations are captured. In 2003, we observe the lowest MB values reflecting the exceptional summer heat wave of 2003. Absolute mean MB errors |MBE s_cal | are the highest for 2007 (448 mm w.e. a−1) and 2003 (356 mm w.e. a−1). These summer errors are higher than winter errors but still inferior to σ s_obs of 2007 and 2003 (respectively 594 and 559 mm w.e. a−1) (Table 2). Larger errors in summer can be partly explained by higher interannual variability in summer observed MBs. If we consider the annual MB (sum of winter and summer MB; Fig. 5c), the largest errors occur for 2007 and 2003 (respectively 582 and 411 mm w.e. a−1), as in summer.
Fig. 5. Time series of mean observed MB (red) over the calibration period and of mean VGT MB (blue) over the period covered by SPOT/VGT, averaged for the 55 glaciers. The dashed red curves represent the time series of observed MB±the standard deviation for all glaciers. (a) Winter MB for the calibration period, 1999–2008 and the period covered by SPOT/VGT, 1999–2014. (b) Summer MB for the calibration period, 1998–2008 and the period covered by SPOT/VGT, 1998–2013. (c) Annual MB (sum of winter and summer MBs) over 1999–2013.
Table 2. Mean mass balance error MBE averaged for all glaciers and standard deviation σ of the observed MBs, per year and for the overall calibration period 1998/1999–2008
The agreement between VGT MB estimations and observed MB over the calibration period presented above is satisfying. In order to test the robustness of our approach, we first perform a cross-validation of the 55 regressions to assess the temporal robustness and skills of the regression coefficients. We then compare seasonal VGT MB time series and independent observed MB data (not used for calibration) over the evaluation period 2009–2014.
4.2. Cross-validation over the calibration period 1998/1999–2008
Cross-validation results for the 55 optimized relations initially calibrated over 1998/1999–2008 indicate no negative skill score SS, except for Seewjinen in summer (SS s = −0.08; Fig. 6b). This is consistent with the low R 2 s and high RMSE s_cal calculated for this glacier over the calibration period (Table S1). For most glaciers, the optimized individual relationships thus have skills to estimate the seasonal MB over a simple average of the observed MB. In winter, SS values are higher (〈SS w〉 = 0.76) than in summer (〈SS s〉 = 0.55) (Table 3, Figs 6a, b). In winter, glaciers from Source 3 perform better in terms of 〈SS w〉 and 〈R 2 w_cross〉 than others. This is consistent with their best performance for the calibration (〈R 2 w_cal〉 = 0.88, Table 1). RMSE w_cross are slightly higher than RMSE w_cal (Fig. 6a), but they remain satisfactory with regard to E ref. In summer, RMSE s_cross are superior to RMSE s_cal (Fig. 6b) and also on average slightly superior to E ref (〈RMSE s_cross〉 = 440 mm w.e. a−1; Table 3). Moreover, for both seasons, the higher the SS, the lower the RMSE and the lower the difference between the two RMSEs (derived from both calibration and test-cross) (Figs 6a, b). Only Aletsch Glacier in winter and Cengal Glacier in summer present satisfactory RMSEs (as regards to E ref) and low SS. Therefore, for these two glaciers, despite their RMSE, their calibrated relationships present no particular skill over a simple average of their observed MBs.
Fig. 6. Results of the cross-validation for all glaciers. Skill score as a function of RMSE computed for the calibration (blue) and for the test cross (red), in (a) winter and (b) summer. Alpha derived from the calibration (blue) and from the cross-validation (red) as a function of RMSE computed for the calibration, in (c) winter and (d) summer. Standard deviation of the alpha derived from the test-cross as a function of RMSE computed for the calibration, in (e) winter and (f) summer.
Table 3. Summary of the cross-validation results over the calibration period 1998/1999–2008 according to individual data sources and for the 55 glaciers
In winter, we observe high consistency between α w_cross_best (mean of the N α w_cross) and the calibrated α w (Fig. 6c). Moreover, the standard deviations σ(α w_cross) in winter are low and close to an order of magnitude smaller than for α w_cross_best (Fig. 6e; Table 3); the same applies for σ(β w_cross). These results allow us to conclude that the calibrated relationships are robust and temporally stable in winter (except for Aletsch), despite the shortness of the time series used for the linear regression. The outcomes show the potential skill of the individual relations to accurately estimate the winter MB outside the calibration period. In summer, the consistency between α s_cross_best and α s is also quite high, as in winter (Fig. 6d). Nevertheless, unlike in winter, the standard deviations σ(α s_cross) in summer are high, especially for glaciers with high RMSE s_cal (superior to 350–400 mm w.e. a−1) (Fig. 6f; Table 3). Therefore, in summer, we can conclude on the robustness and the temporal stability of about 60% of the 55 calibrated relationships (presenting an RMSE s_cal < 400 mm w.e. a−1 and a skill score >0.65). In order to test the robustness of our approach outside the calibration period 1998–2008, we then compare seasonal VGT MB time series and independent observed MB data (not used for calibration) over the evaluation period 2009–2014.
4.3. Performance over the evaluation period 2009–2013/2014
The winter and summer MB averaged for all 55 glaciers estimated with SPOT/VGT (in blue) after 2008 are shown in Figure 5. In winter (Fig. 5a), 〈B w_VGT〉 values are generally above average over 2009–2014, except for 2012, where we note a strong decrease in 〈B w_VGT〉. In summer (Fig. 5b), 〈B s_VGT〉 is close to the average of 1998–2008, except for a notably negative MB during summer 2012. Table S2 (Supplementary Material) presents for each glacier, estimates of both winter and summer MB over 2009–2014 and 2009–2013 respectively. In order to evaluate these VGT MB calculations (over 2009–2014 for winter and 2009–2013 for summer), we now compare them with observed MBs available for a subset of the 55 glaciers used for the calibration (19 glaciers in winter and 13 in summer).
For winter, 70% of the computed MB present an error smaller than the estimated uncertainty in the direct measurements E obs max (Fig. 7a). We find the highest errors MBE w_eval for Sarennes and Ciardoney (−1015 mm w.e. in 2009 and −850 mm w.e. in 2011). This result is not surprising as these glaciers are subject to a high RMSE w_cal over the calibration period (Table S1). 2011 and 2013 are the most poorly estimated years (RMSE w_eval of 587 and 501 mm w.e.), whereas 2010 and 2014 are best represented (RMSE w_eval of 252 and 245 mm w.e.; Table 4). RMSE w_eval calculated from all data in the evaluation period (n = 66) is 411 mm w.e. a−1 and thus nearly three times larger than the average RMSE w_cal (150 mm w.e. a−1) computed for the same subset of 19 glaciers over the calibration period. However, RMSE w_eval is still comparable with E obs max and the mean MB error MBE w_eval is close to zero (16 mm w.e. a−1), suggesting that B w_VGT estimations are unbiased. To sum it up, our approach is able to estimate the winter MB out of the calibration period for glaciers in the Alps with an acceptable mean overall error and without bias.
Fig. 7. Seasonal mass balances estimated with SPOT/VGT as a function of observed (a) winter and (b) summer MBs for glaciers over the evaluation period 2009–2014. The 1:1 agreement is plotted (bold line). The uncertainty in each MB measurement E obs max (±400 mm w.e. a−1) is not represented for the sake of clarity.
Table 4. MB error MBE and RMSE per year and for the overall evaluation period 2009–2013/2014
In summer, 60% of the computed MB present an error inferior or equal to E obs max (Fig. 7b). Largest errors MBE s_eval are about two (to three) times superior to E obs max (e.g. Gietro with an error of 1695 mm w.e. a−1 in 2011) (Table 4). 2009 and 2011 are poorly estimated (resp. RMSE s_eval of 756 and 660 mm w.e.), but summer 2013 is correctly reproduced (low RMSE s_eval of 230 mm w.e.). RMSE s_eval calculated for the 49 evaluation points (561 mm w.e. a−1) is ~1.5 times larger than the RMSE calculated over the calibration period for the subset of 13 glaciers used for evaluation (RMSE s_cal = 368 mm w.e. a−1). RMSE s_eval is also higher than winter RMSE w_eval and than E obs max. The mean MBE s_eval for the 49 points is not negligible (162 mm w.e. a−1), which indicates that B s_VGT estimations are slightly positively biased during 2009–2013. This positive bias is mainly due to a strong bias (664 mm w.e. a−1) during summer 2009. No obvious explanation was found for this anomalous year. Excluding summer 2009, the MBE s_eval is reduced to 52 mm w.e. a−1. However, validation is made with fewer points in summer than in winter. Furthermore, glaciers used for evaluation in summer are not the best performers over the calibration period: the mean RMSE s_cal of the 13 evaluation glaciers (368 mm w.e. a−1) is higher than the mean RMSE s_cal of the entire (55 glaciers) dataset (318 mm w.e. a−1).
Observed MB measurements for more glaciers and more years will be welcome to improve the robustness of the linear regressions between MB and Z in summer and to provide a more representative error RMSE s_eval on the MB estimated with SPOT/VGT.
5. DISCUSSION AND PERSPECTIVES
5.1. Method performance
Our method performs better in winter than in summer because it is based on the interannual dynamics of the altitudinal snow cover distribution to retrieve the interannual variation of seasonal MBs. In summer, the interannual dynamics are more difficult to capture as there is less snow. The lower number of snow-covered pixels in summer than in winter constitutes in fact a limit to retrieve smooth curves of altitudinal NDSI distribution. To justify this hypothesis, we estimate for each season, the mean number of pixels N p with a NDSI higher than 0.2 in each optimized WOSM sizes centred on Clariden Glacier. This glacier was chosen for its much lower performance in summer (R 2 s = 0.52; RMSE s_cal = 528 mm w.e. a−1) than in winter (R 2 w = 0.92; RMSE w_cal = 130 mm w.e. a−1). At Clariden, the optimized WOSM for winter (29 km × 29 km), gives a value of N p that is ~6 times higher in winter (840) than in summer (144). With the optimized summer WOSM (225 km × 225 km), N p is ~12 times higher for winter (19 971) than for summer (1541). Thus, N p increases with WOSM size for both seasons.
This justifies the enlargement of the optimized WOSM* in summer in order to integrate more snow cover variability (Fig. 8). The distribution of WOSM* sizes is more spread towards larger windows: the median WOSM* side length is 119 km in summer against 43 km in winter. The need to increase the window size in summer suggests relatively homogeneous snow variations in summer across the Alps. This result is in agreement with the homogeneous summer ablation observed on glaciers across the entire Alpine Arc (Vincent and others, Reference Vincent2004). Pelto and Brown (Reference Pelto and Brown2012) also noted a similarity of summer ablation across glaciers in the North Cascades.
Fig. 8. Distribution of glaciers (%) as a function of their WOSM* side length (a) in winter and (b) in summer.
Another reason for the difference in results between the seasons might be the cloud interpolation. For each season, we calculate the percentage of interpolated pixels for glaciers with low RMSE (inferior to the first quartile value) and high RMSE (superior to the third quartile), averaged over the calibration period. In winter, the mean percentage of interpolated pixels is large (4%) and the difference between high and low RMSE glaciers is small: the mean difference over the 1999–2014 period is 0.13%. In summer, a striking difference is observed between glaciers with high and low RMSE: the mean difference of the interpolated pixels percentages (0.64%) is more than four times larger than in winter. Therefore, the temporal interpolation seems to have more impact in summer than in winter. Our hypothesis to explain this observation is that the cloud mask performs less well for some glaciers in summer, but this needs to be further explored.
The effect of pixels contaminated by undetected clouds can also be another reason explaining the reduced performance for summer. The spectral signatures in the available bands of clouds and snow are close, indicating that cloudy and snow-covered pixels have similar NDSI values. The error in a pixel NDSI value caused by a cloud is smaller for a snow-covered pixel than for a snow-free pixel. As there is less snow in summer than in winter, undetected clouds may impact more NDSI and snow detection in summer.
Despite the difference in performances between the seasons, we can highlight the ability of the approach to capture the extreme summer balance in 2003 and the two contrasting winter balance years in 2000 and 2001.
The performance of our method during the calibration (in terms of R 2 and RMSE) also varies with the source of the MB data. The best results are obtained with glaciers from source 3 (particularly for winter). The year-to-year variability in this dataset is based on meteorological time series (temperature, precipitation) and distributed modelling. As NDSI and interannual snow cover dynamics are also driven by meteorological parameters, we expected to find a good agreement but we are unable to conclude on the robustness of the method to estimate MB. However, a satisfying mean error for the evaluation period, assessed only with MB data from sources 1 and 2, strengthens the robustness of the method. In fact, these MBs are closer to ‘reality’ as they are composed of in situ glaciological measurements. Nevertheless, MB data from sources 1 and 2 are also subject to uncertainties. In particular, for Aletsch Glacier (from source 2), the largest glacier of the Alps (83 km2), the lowest R 2 w (0.31) for winter was found. A possible explanation for the poor performance of our approach for Aletsch is that its ablation area reaches to relatively low elevation. It is thus snow-free for a considerable part of the period used for winter balance estimate (October to April). Winter MB data for Aletsch Glacier thus include both accumulation and melting, whereas the NDSI approach for winter is optimized to represent snow accumulation, as for the other glaciers.
5.2. Perspectives
A limitation of the proposed methodology is the impossibility (in some rare cases) to estimate the seasonal MB after the calibration period. This happens only in summer, specifically in 2009 and 2011 for Sarennes, and for Vernagt in 2012. In those particular cases, the NDSI value optimizing the linear relationships between Z and MB over the calibration period 1998–2008 does not intersect the NDSI altitudinal distribution curves: in both cases the curves are below the NDSI value. Thus, no MB can be retrieved. We note that for these two glaciers the optimized WOSM side lengths are among the smallest (5 km for Sarennes and 21 km for Vernagt). If we increase the WOSM sizes, we can recover summer MB estimations for these years and these two glaciers but at the cost of degrading the regression quality (R 2 s and RMSE s_cal) over the calibration period, mainly for Vernagt (R 2 s decreases from 0.75 to 0.60 and RMSE s_cal increases from 262 to 332 mm w.e. a−1). An objective criterion is being analysed to determine a WOSM side length threshold >5 km in order to get NDSI altitudinal curves representative enough of the altitudinal snow cover distribution.
An important outlook of this study is to achieve a better cloud detection as an under/over estimation of cloudy pixels impacts snow-cover monitoring, especially in summer.
Another perspective to improve our summer results could be to couple our methodology with a melt model that is able to provide information on the short-term dynamics of melting.
One of this study's aim is to perform real-time seasonal MB monitoring with daily low-resolution optical imagery for a large sample (here 55) of Alpine glaciers. Low resolution optical satellite images are available 2–3 d after the acquisition dates, implying that seasonal snow cover and seasonal MB could be estimated a few weeks after the end of each season. Moreover, the method could also allow for glaciers with recently initialized MB series to extend them backwards in timer, for example by using the Advanced Very High Resolution Radiometer satellite data available since 1978. Interrupted MB time series could be reconstructed as well.
The ultimate and long-term goal of this work is to apply the method at large scales (i.e. thousands of glaciers in a mountain range for which no or very few direct MB data are available). Therefore, the next step is to assess how accurately the approach can estimate the MB of an unmeasured glacier. For that purpose, it is necessary to build and apply a generic transposable relation between Z and MB with fixed parameters. The sensitivity of the coefficients of the linear regressions towards explicative factors (in particular topographic factors such as glacier size, aspect, etc.) or the variations of the optimized NDSI value and WOSM size need to be investigated in order to determine explicit and objective criteria to build a generic relation. Currently, the variability of alpha and beta remains a little too large (even in winter) to easily and in a satisfactory way obtain a simple generic relationship.
6. CONCLUSION
In this study, we have described an empirical method to estimate seasonal MBs of Alpine glaciers from kilometric resolution optical SPOT/VGT images. From seasonal snow cover maps of the Alps derived for each year over 1998–2014, a regional mean snow altitude of a region surrounding each glacier is derived for 55 glaciers with MB data. Promising linear relationships between this regional mean snow altitude and observed seasonal MBs have been found over 1998–2008. The explained variance in winter is high for all glaciers (R 2 = 0.84 on average) and the mean RMSE is low (161 mm w.e. a−1). Results are not as good for summer but the explained variance is acceptable (R 2 = 0.73) and the mean RMSE (318 mm w.e. a−1) is still in the range of errors associated with glaciological MB measurements (typically from ±200 to 400 mm w.e. a−1). Cross-validation of the 55 individual linear regressions allows assessment of the temporal stability and the robustness of all the derived relationships in winter and of ~60% in summer. Estimations of seasonal MB over 2009–2014 is also performed for these 55 glaciers, and a mean global error is calculated for some estimates, based on a more limited dataset of observed MBs. We are able to estimate winter MB with an acceptable mean global error (405 mm w.e. a−1) and without bias. In summer, the mean error of MB estimation over 2009–2014 is higher (561 mm w.e. a−1) but less observed MB measurements are available for validating the regressions during this season. Moreover, the subset of glaciers used for evaluation in summer tends to perform more poorly during calibration. Still, these results are promising and estimates of the seasonal MB could be performed as soon as the SPOT/VGT data (or data from similar satellites) are available and processed. A real-time seasonal MB ‘monitoring’ is thus conceivable, a few weeks after the end of each season for a large sample (here 55) of Alpine glaciers.
Our method performs better in winter than in summer. This is mainly explained by the fact that interannual dynamics of altitudinal snow cover distribution is more difficult to capture in summer as there is less snow. We also highlight the fact that interpolation of cloudy pixels has more impact in summer than in winter. However, an in-depth analysis needs to be carried out. The greater accuracy in winter emphasizes the value of this method over the classical snowline approach that does not provide any estimate in winter when accumulation area ratio is 100%.
The SPOT/VGT mission ended in May 2014 but the PROBA-V satellite, launched in May 2013, ensures continuity and has been providing data since October 2013, with images at 1 km, 300 m and 100 m resolution. A comparison of snow cover and MB estimates derived from SPOT/VGT and from PROBA-V 1 km for the period of overlap (winter 2013–2014) is underway in order to extend MB data after 2015.
Supplementary material
The supplementary material for this article can be found at http://dx.doi.org/10.1017/jog.2016.78.
ACKNOWLEDGEMENTS
We acknowledge the WGMS for the mass-balance data. We also thank all the observers who contributed to collect the seasonal mass balances and shared them with the community through the WGMS database. We are grateful to VITO and Belspo for the SPOT/VGT satellite images distribution, and the financial and scientific support to carry out this study. We acknowledge support from the CNES/TOSCA programme (in particular Juliette Lambin) as well as funding of a Ph.D. fellowship by VITO/CLS (and especially Eric Gontier, VITO and Estelle Obligis, CLS). We thank the Scientific Editor, David Rippin and two anonymous reviewers for their comments and suggestions which significantly improved the manuscript. This article is also a tribute to Gilbert Saint and his vision during the early 90s, of what an operational satellite mission should be.
Open access
