Changes in glacier equilibrium-line altitude in the western Alps from 1984 to 2010 : evaluation by remote sensing and modeling of the morpho-topographic and climate controls

We present time series of equilibrium-line altitude (ELA) measured from the end-of-summer snow line altitude computed using satellite images, for 43 glaciers in the western Alps over the 1984–2010 period. More than 120 satellite images acquired by Landsat, SPOT and ASTER were used. In parallel, changes in climate variables, summer cumulative positive degree days (CPDD) and winter precipitation, were analyzed over the same time period using 22 weather stations located inside and around the study area. Assuming a continuous linear trend over the study period: (1) the average ELA of the 43 glaciers increased by about 170 m; (2) summer CPDD increased by about 150 PDD at 3000 m a.s.l.; and (3) winter precipitation remained rather stationary. Summer CPDD showed homogeneous spatial and temporal variability; winter precipitation showed homogeneous temporal variability, but some stations showed a slightly different spatial pattern. Regarding ELAs, temporal variability between the 43 glaciers was also homogeneous, but spatially, glaciers in the southern part of the study area differed from glaciers in the northern part, mainly due to a different precipitation pattern. A sensitivity analysis of the ELAs to climate and morpho-topographic variables (elevation, aspect, latitude) highlighted the following: (1) the average ELA over the study period of each glacier is strongly controlled by morphotopographic variables; and (2) the interannual variability of the ELA is strongly controlled by climate variables, with the observed increasing trend mainly driven by increasing temperatures, even if significant nonlinear, low-frequency fluctuations appear to be driven by winter precipitation anomalies. Finally, we used an expansion of Lliboutry’s approach to reconstruct fluctuations in the ELA of any glacier of the study area with respect to morpho-topographic and climate variables, by quantifying their respective weight and the related uncertainties in a consistent manner within a hierarchical Bayesian framework. This method was tested and validated using the ELA measured on the satellite images.


Introduction
Historically, the glacier's annual surface mass balance and equilibrium-line altitude (ELA) have been computed from direct field measurements of snow accumulation and snow/ice ablation through point measurements using a network of ablation stakes and snow pits on individual glaciers.The longest mass balance data series started in the late 1940s for Storglaciären (Sweden), Taku Glacier (USA), Limmern and Plattalva glaciers (Switzerland), Sarennes Glacier (France) and Storbreen (Norway).The World Glacier Monitoring Service (WGMS) has mass balance data on around 260 glaciers worldwide (note that the World Glacier Inventory contains more than 130 000 glaciers), but uninterrupted time series spanning more than 40 yr are available for only 37 glaciers (WGMS, 2012).The small number of available data series is due to the cost, in terms of money, manpower, and time, of this laborious method.The small number of data series is an obstacle that needs to be overcome if we are to improve our knowledge of the relationship between climate and A. Rabatel et al.: Changes in glacier equilibrium-line altitude in the western Alps glacier changes from mountain-range up to regional scale, as well as our knowledge of the contribution of glaciers to water resources in the functioning of high-altitude watersheds.Remote sensing provides a unique opportunity to address the question of glacier changes at regional scale.Since the availability of satellite data in the 1970s, several methods have been developed to compute changes in glacier volume at multiannual to decadal scale using digital elevation models (DEMs) (e.g.Echelmeyer et al., 1996;Baltsavias et al., 1999;Berthier et al., 2004;Gardelle et al., 2012), or the annual ELA and mass balance (e.g.Østrem, 1975;Rabatel et al., 2005Rabatel et al., , 2008Rabatel et al., , 2012a)).However, satellite derived DEMs are not accurate enough to compute variations in the annual volume of mountain glaciers.Another alternative emerged from recent studies based on ELA and distributed mass balance modeling using meteorological input fields derived from local weather stations, data reanalysis, or regional climate models (e.g.Zemp et al., 2007;Machguth et al., 2009;Marzeion et al., 2012).However, such modeling approaches still have a number of limitations, e.g.accumulation is underestimated in mountainous regions.On the other hand, for mid-latitude mountain glaciers, the end-of-summer snow line altitude (SLA) is a good indicator of the ELA and thus of the annual mass balance (Lliboutry, 1965;Braithwaite, 1984;Rabatel et al., 2005).This enables ELA changes to be reconstructed for long time periods from remote-sensing data (Demuth and Pietroniro, 1999;Rabatel et al., 2002Rabatel et al., , 2005Rabatel et al., , 2008; Barcaza et al., 2009;Mathieu et al., 2009), because the snow line is generally easy to identify using aerial photographs and satellite images (Meier, 1980;Rees, 2005).Consequently, it is possible to study the climate-glacier relationship at a massif-or regional scale (Clare et al., 2002;Chinn et al., 2005), which is particularly useful in remote areas where no direct measurements are available.
In the current study, we rely on previous studies conducted in the French Alps (Dedieu and Reynaud, 1990;Rabatel et al., 2002Rabatel et al., , 2005Rabatel et al., , 2008)), to reconstruct ELA time series for more than 40 glaciers over the 1984-2010 period, using the end-of-summer snow line detected on satellite images.Our aim are (1) to quantify at a regional scale the temporal and spatial changes of the ELA; (2) to characterize the relationships between ELA and both morphotopographic and climate variables; and (3) to reconstruct the spatio-temporal variability of annual ELA time series by incorporating the above-mentioned relationships in an expansion of Lliboutry's variance decomposition model (1974).
2 Study area, data and method

Selection of the study sites
A recent update of the glacier inventory of the French Alps lists 593 glaciers (Rabatel et al., 2012b).The total glacial coverage in the French Alps was about 340 km 2 in the mid- 1980s and had decreased to about 275 km 2 in the late 2000s.In the current study, 43 glaciers located in the French Alps or just next to the border with Switzerland and Italy were selected (Table 1, Fig. 1).The selection was based on the following criteria: (1) glaciers had to have a high enough maximum elevation to allow observation of the snow line every year during the study period; (2) glaciers with all aspects had to be represented; and (3) glaciers located in all the glacierized mountain ranges in the French Alps from the southernmost (44 • 50 N) to the northern-most (46 • 00 N) had to be included.
Among the selected glaciers, three belong to the GLACIO-CLIM observatory, which runs a permanent mass balance monitoring program.These three glaciers are Saint Sorlin Glacier, monitored since 1957 (# 31 in Table 1 and Fig. 1), Argentière Glacier monitored since 1975 (# 4 in Table 1 and Fig. 1), and Gébroulaz Glacier monitored since 1983 (# 28 in Table 1 and Fig. 1).Mass balance data available on these Table 1.List of the glaciers studied and the weather stations used in this study.Numbers (for the glaciers) and codes (for the weather stations) refer to Fig. 1.All glaciers are located in the French Alps except when specified in brackets: CH = Switzerland and IT = Italy.For the weather stations, letters in brackets identify the data producer: MF = Météo-France, EDF = Electricité de France, MS = Météo-Suisse.The coordinates are given in degrees, minutes, for latitude and longitude and in m a.s.l. for elevation.

Glaciers
Weather three glaciers were used to assess the representativeness of the end-of-summer SLA computed from the satellite images as an indicator of the ELA and of the annual mass balance computed from field measurements (see Sect. 3.1 below).

Satellite images and snow line delineation
A total of 122 images of the 43 glaciers were used to cover the 27 yr study period.Unfortunately, in some years usable images were not available for all the glaciers because of (1) cloud cover hiding the underlying terrain; and (2) snowfalls that can occur in late summer and completely cover the glaciers.This was mainly the case in the 1980s and 1990s when fewer satellites were in orbit than in the 2000s.The images we used were acquired by the following satellites: Landsat 4TM, 5TM, 7ETM+, SPOT 1 to 5 and ASTER, with spatial resolutions ranging from 2.5 to 30 m (see supplementary material: Table 1 for a detailed list of the images used).
For the snow line delineation on multispectral images, a test of band combinations and band ratios to facilitate the identification of the snow line on the images is described in Rabatel et al. (2012a, see Fig. 4).These authors concluded that the combination of the green, near-infrared, and shortwave infrared bands (see Table 2 for details about the wave- lengths of each satellite) was the most appropriate to identify the limit between snow and ice.The snow line was delineated on the central part of the glaciers to avoid border effects on the glacier banks: shadows from surrounding slopes, additional snow input by avalanches, overaccumulation due to wind (Fig. 2), which could generate equilibrium-line position dependence on local conditions (Rabatel et al., 2005).
The delineation was performed manually because automatic methods hardly succeed in distinguishing the snow line from the firn-line when both to be presented on the glacier.This distinction results to be also visually difficult when the pixel size of the satellite images is too coarse (see the discussion section).

Digital elevation model, computation of the altitude of the snow line and of the glacier morpho-topographic variables
The average altitude of each snow line was calculated by overlaying the shapefiles containing the digitized snow lines on a DEM.For each glacier, the DEM was also used to compute the following morpho-topographic variables used in this study: surface area, mean elevation, aspect, latitude and longitude of the glacier centroid.The mean elevation of each glacier has been computed as the arithmetic mean of the elevation of each pixel of the DEM included within the glaciers outline.This mean elevation is rather close to the median elevation of the glaciers.Indeed, the difference between the two variables is 10 m on average for the 43 studied glaciers.This shows that the studied glaciers have on average an almost symmetrical area-altitude distribution (Braithwaite and Raper, 2010).Several DEMs are available for the study area.Among these, one comes from the French National Geographical Institute (IGN) with a resolution of 25 m, and dates from the early 1980s.Another one is the ASTERGDEM with a resolution of 30 m, and dates from the mid-2000s.The average difference between the two DEMs at the mean elevation of the SLA of the 43 glaciers studied over the whole study period, i.e. approx.3050 m a.s.l., is about 20 m.This matches the vertical accuracy of the DEMs, and is small in comparison to the interannual variability of the SLA (the standard deviation of the measured SLA over the whole study period ranges between 75 and 255 m depending on the glacier).However, because the glacier surface lowering can reach several tens of meters at lower elevation (up to 80 m at 1500 m a.s.l.) due to the important glacier shrinkage over the last decades (Thibert et al., 2005;Paul and Haeberli, 2008;Vincent et al., 2009;Berthier and Vincent, 2012), when the SLA is at lower elevation as was the case in 2001 (2839 ± 129 m a.s.l.) it is better to use a DEM as close as possible to the date of the used images.Accordingly, the DEM from the French IGN was used for the first half of the period (till the late 1990s) and the ASTERGDEM was used for the second part of the period.Changing from one DEM to the other in the late 1990s has no impact on the results because at that time, the SLA was located between 3000 and 3200 m a.s.l., an altitudinal range for which the difference in elevation between the DEMs is lower than their vertical accuracies.
An uncertainty was estimated for each SLA.Uncertainty results from different sources of error (Rabatel et al., 2002(Rabatel et al., , 2012a)): (1) the pixel size of the images, which ranged between 2.5 and 30 m depending on the sensor; (2) the slope of the glacier in the vicinity of the SLA, which ranged from 7 • to 32 • depending on the glacier and the zone where the SLA was located in any given year; (3) the vertical accuracy of the DEM, which was about 20 m; and (4) the standard deviation of the calculation of the average SLA along its delineation, which ranged from 10 m to 170 m depending on the glacier and on the year concerned.The last source of error (4) is the most important.The resulting total uncertainty of the SLAs (root of the quadratic sum of the different independent errors) ranged from ±15 to ±170 m, depending on the year and glacier concerned.The uncertainty was greater when the SLA was located on a part of a glacier where the slope is steepest because in this case, the standard deviation of the computed SLA is high.

Method used to fill data gaps in the snow line altitude time series
In spite of the large number of images used in this study, there were still some gaps in the SLA time series.In the 27 yr study period, for the 43 glaciers studied, the SLA was measured in 85 % of the 1081 cases.Thus, data were missing in 15 % of the cases.
To fill these gaps, a bi-linear interpolation was applied.Originally, this approach was used by Lliboutry (1965Lliboutry ( , 1974) ) to extract the spatial and temporal terms from a very incomplete table of mass balance measurements (70 % gaps) in the ablation zone of Saint Sorlin Glacier.This approach is based on the fact that the glacier mass balance depends on both the site and the time of measurement.This approach was subsequently extended to a set of glaciers in the Alps by Reynaud (1980), and later on, to glaciers in other mountain ranges by Letréguilly and Reynaud (1990).Here, we applied it to our SLA table that had far fewer gaps, writing simply where SLA it is the snow line altitude of the glacier i for the the year t, α i is long-term mean for each glacier over The Cryosphere, 7, 1455-1471, 2013 www.the-cryosphere.net/7/1455/2013/ the period of record, β t is a term (depending on the year only) which is common to all the glaciers analyzed, and , a table of centered residuals assumed to be independent and Gaussian (N denotes the normal distribution).The purpose is not simply to obtain the α i and β t terms, but to estimate each missing value SLA it on glacier i and year t as: where T is the number of years of the study period, M is the number of glaciers studied, and the circumflex is a classical statistical estimate.
To maximize the consistency of the reconstruction, αi was calculated using only the years for which no data was missing.However, only six years (out of the 27) had the SLA for each of the 43 glaciers, which is not quite enough to obtain an average SLA that is not too dependant on a specific year.Consequently, the initial data set was divided into three geographical subsets: the Ecrins Massif, the Vanoise area, and the Mont Blanc Massif, and the computation was applied independently to each one of these subsets.In this way, more years with SLAs for all the glaciers of each subset were available (see Table 3, where the statistics of each subset are presented).The year 1995 was the only exception.For this specific year, because very few values were available, the reconstruction was made considering the three subsets together.
The RMSE of ε it , representing the model errors, was 73 m.This is less than half the standard deviation of the measured SLAs (174 m), showing that the use of the bi-linear model makes sense on this data set.

Meteorological data
Table 1 also lists the weather stations we used, which are located in and around the study area (Fig. 1).Among the 105 weather stations available in the Météo-France database (including stations operated by Météo-France and Electricité de France) in the French departments (Hautes-Alpes, Isère, Savoie and Haute-Savoie) in the study area, 40 were selected on the basis of data availability and their proximity to the glacierized areas (except the Lyon-Bron weather station, which was selected as a regional reference).Among the 40 weather stations, only 20 of which had a continuous series of winter precipitation and summer temperatures over the period 1984-2010, and were finally used.Two weather stations in Switzerland were also used.The climate variables used in the analysis were (1) cumulative positive degree days (CPDD) from 15 May to 15 September for each year t, extrapolated to the altitude of 3000 m a.s.l.using a standard gradient of 6 • C km −1 (3000 m a.s.l.being the approximate mean elevation of the SLA of the 43 glaciers studied over the whole study period); and (2) cumulated winter precipitation from 15 September of the year t-1 to 15 May of the year t.During this period, liquid precipitation is negligible at 3000 m a.s.l.
For the CPDD, the date of 15 September was chosen because it matches the average date of the satellite images used to delineate the snow line.Furthermore, it should be noted that at 3000 m a.s.l. in our study area, after 15 September, CPDD are negligible (Thibert et al., 2013).

Modeling of the control of the morpho-topographic and climate variables
From our data set of SLA series, glacier morpho-topographic variables, and climate variables, we propose a simple expansion of Lliboutry's approach (1974) to model and reconstruct the SLA fluctuations of any glacier in our study area over the study period.This approach also enables us to quantify the respective control of each morpho-topographic and climate variables on the SLA spatio-temporal variability.Following the ideas of Eckert et al. (2011), the model is implemented in a hierarchical Bayesian framework (e.g.Wikle, 2003).The main advantage over more empirical approaches is treating the different sources of uncertainty in a consistent manner.For instance, the available information (missing values and repeated observations among correlated glacier series) is "respected" when inferring the spatio-temporal patterns of interest, more specifically: the regression variables relating ELA fluctuations, their morphotopographic and climate drivers, and their respective weight.
In detail, starting from the available SLA table with hollow (i.e.without reconstructed values) and Eq.(1), we further decomposed the spatial term as: where a o is the interannual regional average, X ki are the spatial explanatory variables considered for each glacier, i.e. the morpho-topographic variables: latitude, average altitude, surface area and aspect, and V i ∼ N 0, σ 2 v is the spatial residual, i.e. the local effect for the glacier i which is not explained by the spatial explanatory variables considered.
Similarly, the temporal term is decomposed as: where the X kt are the temporal explanatory variables considered, i.e. the climate variables: summer CPDD and winter precipitation, and Z t ∼ N 0, σ 2 z is the temporal residual, i.e. the annual effect common to all glaciers for any year t, which is not explained by the temporal explanatory variables considered.
For the sake of simplicity, the model is fed with reduced standardized variables (by dividing by the standard deviation the difference between each value and the average of the series).This avoids identifiability problems by granting that t β t = 0, and makes the a k and b k coefficients more directly interpretable, at least for non-heavily correlated explanatory variables (see below).
Without the V i and Z t terms, the model would be a classical spatio-temporal linear model with easy inference, e.g. with likelihood maximization providing analytical solutions for point estimates and asymptotic standard errors.However, this would not allow us to distinguish between explained and unexplained variance in the spatial and temporal terms.Hence, the V i and Z t terms make the model hierarchical, i.e. richer but also more tricky to estimate.This challenge was solved using Bayesian Markov chain Monte Carlo methods (MCMC, Brooks, 1998).As detailed in Eckert et al. (2007Eckert et al. ( , 2010b) ) for similar problems, the robustness of the inference was checked using tests based on starting different simulation runs at different points of the parameter space (Brooks and Gelman, 1998).Poorly informative prior probability distributions were used for all parameters, however, this standard choice allows us to obtain posterior estimates only driven by information conveyed by the data.From the joint posterior distribution of all parameters, latent variables, and missing values, we retained point estimates (the posterior mean), posterior standard deviations, and 95 % credibility intervals, the Bayesian counterpart of the classical confidence interval (Table 4).The construction of the model ensures that α i , β t and ε it are independent.Furthermore V i and k X ki a k are independent, so are Z t and k X kt b k .Hence the total variance of the data table is: By analogy to the classical R 2 , variance ratios can then be computed to evaluate the respective weight of the different Similarly, ratios can be computed for the spatial and temporal terms solely: and, getting rid of the random spatio temporal fluctuations, the overall time/space ratio and determination coefficient in the decomposition can be computed as: All these statistics quantify the weight of mean spatial and temporal effects in the total variability, and the capacity of the chosen explanatory variables to model these mean effects.At the observation level of each year/glacier of the model, local and annual adjustment statistics can also be computed as: respectively representing the annual and local bias, and the annual and local determination coefficient, with M the number of glaciers considered and T the number of years studied.They quantify the goodness of fit of each of the annual/local series.
Variance ratios and adjustment statistics were evaluated by computing their values at each point of the MCMC iterative simulation run, giving access to point estimates and related uncertainties taking into account the number of missing values for each year/glacier (logically, when nearly all observations are missing for a given year/glacier, the corresponding determination coefficients δ 2 t and δ 2 i are unknown).

Results
Here we give the results of both the SLA and the climate variables.In Sect.3.1, we present the temporal and spatial variability of the SLA measured on the 43 glaciers for the 1984-2010 period.In Sect.3.2., we present the variability of climate variables (i.e.summer CPDD and winter precipitation).

Changes in the equilibrium-line altitude
As mentioned above, it has been demonstrated that on midlatitude glaciers where superimposed ice is negligible, the end-of-summer SLA is an accurate indicator of the ELA (see Fig. 6 in Rabatel et al., 2005, andFig. 3B in Rabatel et al., 2008).Because of this strong correlation and for the sake of simplicity, we hereafter only use the term ELA.

Temporal variability of the equilibrium-line altitude
Figure 3a illustrates changes in the ELA over the 1984-2010 period for the 43 glaciers studied (see also see Supplement: Table II where all the data are presented).Considering all the glaciers, the average ELA for the whole period was located at 3035 ± 120 m a.s.l., i.e. the interannual variability of the average ELA of all glaciers was high (120 m on average).The difference between extreme years was as high as 460 m, with the lowest average ELA measured in 1984 (2790 ± 180 m a.s.l.), and the highest average ELA measured in 2003 (3250 ± 135 m a.s.l.).In addition, over the study period, the ELA time series showed an average increasing trend of 6.4 m yr −1 , assuming a linear trend which results to be statistically significant considering a risk of error of 5 %.This is the equivalent of an average increase of 170 m over the 1984-2010 period, i.e. higher than the interannual variability of the average ELA.
The average annual ELA during the first five years of the study period was lower than the average ELA for the whole period (except in 1986, when it was slightly higher than the average).These years match the end of a 15 yr period (from the mid-1970s to the late 1980s) during which alpine glaciers increased in size due to higher winter accumulation and reduced summer ablation (Vincent, 2002;Thibert et al., 2013).Since 2003, the average annual ELA has consistently been above the average ELA for the whole period.Considering the two subsets, i.e. before and after 2003, no significant trend appeared in either of the subsets.The second subset is shorter (8 yr) than the first (19 yr), and trends computed on short periods should be interpreted with caution.However, the average ELA between the two subsets increased by about 140 m, which is 82 % of the total change.Hence, the year 2003 can also be considered as a breakpoint in the time series.It should be noted that the 2003 extreme heat wave resulted in a reduction, or even complete loss on some glaciers, of the firn area, hence introducing a positive feedback on glacier mass balances after 2003, consistent with higher ELAs.
Considering each glacier individually, Fig. 4 shows the rate of increase in the ELA over the whole period for each www.the-cryosphere.net/7/1455/2013/The Cryosphere, 7, 1455-1471, 2013 glacier plotted against its aspect.In this study, the aspect of the glacier matched the average aspect of the area of the glacier, where the ELA fluctuated over the whole study period.It is notable that (1) the increasing trend of the ELA varied between less than 1 m yr −1 and more than 13 m yr −1 depending on the glacier; and (2) glaciers facing east displayed a more pronounced increasing trend of the ELA.However, this result should be interpreted with caution because the distribution of the glaciers studied were not homogeneous with respect to the aspect of the majority of glaciers in the French Alps, which are north facing.Furthermore, one of our selection criteria should also be recalled: to measure the ELA for each year, the glaciers had to reach at least 3250 m a.s.l.
(highest average ELA measured in 2003), and south-facing glaciers reaching this elevation are rare and can only be found in the Mont Blanc Massif.

Spatial variability of the equilibrium-line altitude
Figure 5a shows the results of a factor analysis of the ELA time series, including all the glaciers and all the years in the study period.The factor analysis is a statistical method allowing to describe the variability, differences or similitudes among observed and correlated variables according to a lower number of factors, i.e. unobserved variables, for instance the space and time.For the ELA time series, the first two factors explained 80 % of the common variance with 71 % on the first factor.The high percentage of common variance on factor 1 represents the common temporal signal in the interannual variability of the ELAs that is related to climatic forcing.Factor 2 explained 9 % of the common variance of the glaciers studied and is related to different spatial patterns.Hence, with respect to the second factor, two groups can be distinguished corresponding to the northern Alps (Mont Blanc, Vanoise) and the southern Alps (Ecrins), respectively.One glacier in the Vanoise Massif, Gébroulaz Glacier, is closer to the glacier group of the southern Alps.Its geographic location, i.e. nearest to the southern sector among the Vanoise glaciers, could partly explain the position of this glacier in the graph.The distinction between the two groups could be related to a higher elevation of the ELA for the glaciers of the southern sector (Fig. 6a) which could be associated with differences in winter accumulation between the southern and northern Alps (see below) and warmer summer temperature in the southern sector, which would increase the amount of CPDD and thus the ablation (note that the effect of higher summer temperature on a shift of precipitation from snow to rain is not considered in this study because we do not use summer precipitation in our analysis).Indeed, Fig. 6a shows the average ELA over the 1984-2010 period for each glacier plotted as a function of its latitude: a clear decreasing trend is apparent from south to north, with an average ELA located on average 150 m higher for the glaciers in the southern sector.

Changes in climate conditions
Figure 3b shows the changes in summer CPDD over the 1984-2010 period for the 22 weather stations used.Summer CPDD present an increasing trend, averaging 5.3 ± 1.9 CPDD yr −1 at 3000 m a.s.l.assuming a linear trend over the period 1984-2010.This trend is statistically significant considering a risk of error of 5 %.Such an increase (150 ± 1.9) in the CPDD over the study period is the equivalent of an additional energy supply of 14 ± 6 W m −2 , which can result in an additional ablation of about 0.5 m w.e.(assuming that all of this energy is used to melt the snow, and so that sublimation is negligible).The year 2003 had the maximum summer CPDD in the 1984-2010 period, and 2003 is an outlier in the time series.In this particular year, the ELA was about 300 meters higher than the average for the whole period (in fact above many glacierized summits), and is clearly associated with an unusually warm summer (Beniston and Diaz, 2004) that caused intense ablation, because winter precipitation in 2003 was average for the study period.
Temporal variations in the summer CPDD at 3000 m a.s.l. were highly homogeneous across the study area.Indeed, the factor analysis of these data (Fig. 5b) showed a common variance of 85 % on the first axis indicating similar temporal variations between the stations located up to 150 km apart.
Regarding winter precipitation, no significant trend was observed over the study period (Fig. 3c).However, one can note that, on average, interannual variability was lower after 2001 (standard deviation divided by 2 after 2001, falling from 140 mm to 62 mm), but it should be kept in mind that the time series after that date is short, and this point thus needs to be confirmed with longer data series.Over the study period, maximum winter precipitation occurred in 2001, and the same value was almost reached in 1995.These two years correspond to low ELAs, but not to the lowest of the whole time series, which occurred in 1984 and 1985, and appears to be associated more with cooler summer temperatures.
Variations in winter precipitation at the different stations over time were not as homogeneous as the summer CPDD, even if the first factor of the factor analysis explained 71 % of total variance (Fig. 5c).As was the case for the ELAs, two groups can be distinguished on the second axis, explaining 9 % of the variance.This factor analysis showed that most of the variance was common to all the weather stations and was due to large-scale precipitation patterns, but that local effects also exist.This is particularly true for the weather stations located in the southern part of the study area, which is influenced by a Mediterranean climate in which some winter precipitation is not only due to typical western disturbances, but also to eastern events coming from the Gulf of Genoa such as the one in December 2008 (Eckert et al., 2010a).

Discussion
In this section, we first discuss the bivariate relationships between each of the morpho-topographic variables considered and the average ELA of each glacier for the study period.Then, in the second part, we present the bivariate relationships between the ELA interannual variability and each climate variable.Finally, we exploit the expansion of Lliboutry's approach (1974) that explicitly incorporates both climate and morpho-topographic variables to reconstruct annual ELAs.

The role of morpho-topographic variables
Figure 6 shows the relationship between the average ELA of each glacier computed over the whole study period and the latitude of the glacier (Fig. 6a), the mean altitude of the www.the-cryosphere.net/7/1455/2013/The Cryosphere, 7, 1455-1471, 2013  glacier (Fig. 6b), the glacier surface area (Fig. 6c), and the glacier aspect (Fig. 6d).
We already mentioned the relationship between the ELA and the latitude (Fig. 6a) in Sect.3.1.2,where we noted that the average ELA of the glaciers located in the southern part of the study area was about 150 m higher in altitude than the ELA of the glaciers located in the northern part of the study area.This meridional effect is consistent with the drier and warmer conditions associated with the Mediterranean climate that prevails in the southern part of the study area.
Figure 6b and c show the relationship between the average ELA and the average altitude of the glacier and the glacier surface area, respectively.The correlation between the average ELA and the glacier average altitude was positive, mean-ing that the lower the average altitude of the glacier, the lower the average ELA.Conversely, the correlation between the average ELA and the glacier surface area was negative, meaning that the smaller the glacier, the higher the ELA.To understand these correlations, one has to keep in mind the relationship between the ELA and the accumulation area ratio (AAR).Indeed, the ELA constitutes the lower limit of the accumulation zone, which represents ∼ 2/3 of the total glacier surface area in a steady-state glacier.Accordingly, the wider the accumulation zone, the bigger the glacier, the lower its snout, the lower the mean altitude and consequently the lower the ELA of the glacier.Current ELAs certainly can not represent steady-state conditions, nevertheless, over a long period of time (almost 30 yr in our case), assuming pseudo-stationary conditions over this period, the average surface area of a glacier as well as its average altitude can be deemed to be representative of its average ELA.This is in good agreement with the results presented by Braithwaite and Raper (2010) showing, on the basis of 94 glaciers, that the glacier area, median and mid-range altitudes are accurate proxies of the ELA for a period of time, e.g.several decades.
Figure 6d shows the average ELA of each glacier plotted as a function of the glacier aspect.For this morphotopographic variable, we distinguished between the different massifs in the study area, because when treated all together, the effects of aspect and latitude tended to cancel each other out.Indeed, the average ELA of a glacier facing north in the Ecrins Massif is almost the same as the average ELA of a glacier facing south in the Mont Blanc Massif.As a result, the meridional effect related to the latitude is stronger than the aspect effect.However, considering each massif separately, and even if our samples are not homogeneous with respect to aspect, glaciers facing north and east have a lower average ELA than glaciers facing south and west, as would be expected for mid-latitude glaciers in the Northern Hemisphere.

Climatic control of equilibrium-line altitude interannual variability
Analysis of the variances and correlations between the ELA of each glacier and the climate variables (summer CPDD and winter precipitation) was conducted for the study period.Considering our 27 yr data set, the correlation coefficient r is statistically significant with a risk of error lower than 1 % when r > 0.49, or r < −0.49, and lower than 5 % when r > 0.38 or r < −0.38.
Figure 7 shows the correlation coefficients for the two climate variables obtained by selecting the weather station with the highest value for each glacier.The correlation coefficients were plotted against the latitude of the glaciers in order to evaluate spatial variations in these coefficients.For summer CPDD, Fig. 7a shows: (1) a generally high and always sig-nificant correlation coefficient regardless of the glacier and its location in the study region; (2) a slight but not statistically significant decrease in the correlation coefficient with increasing latitude.This means that the proportion of variance of the ELAs explained by the summer CPDD is, first, important, and secondly, of the same order regardless of the glacier concerned.
This also appears clearly in Fig. 3 which shows strong similarities between annual variability of ELAs (Fig. 3a) and CPDD (Fig. 3b), mostly an increasing trend over the whole study period, although stronger before 1990, which apparently stopped after the 2003 maximum.Hence, the annual means of both variables are strongly correlated (r = 0.73).To further highlight their similarities, a spline regression was performed, which fit the data variability well, apart from the exceptional CPDD in 2003.The correlation between these low frequency signals was enhanced with regards to annual means, reaching r = 0.84, confirming the very similar longterm variations in ELA and CPDD over the study area.
Regarding winter precipitation, Fig. 7b shows the particular behavior of the Ecrins Massif (blue triangles in Fig. 7b), which could be related to a winter precipitation regime influenced by the eastern disturbances coming from the Gulf of Genoa and presenting higher temporal variability than the typical regime in which disturbances come from the west.However, except for three glaciers in the Ecrins Massif, the correlation between the ELA and cumulative winter precipitation was statistically significant considering a risk of error lower than 5 %, and even 1 % in most cases.These three glaciers are Glacier des Violettes (#40), Glacier de la Girose (#33) and Glacier du Mont de Lans (#35); the low correlation with precipitation for these glaciers could be due to topographic effects influencing redistribution of the snow by the wind which can affect accumulation processes in one way or the opposite, i.e. overaccumulation or erosion (Girose and Mont de Lans are located on a dome and Violettes in a narrow and cached mountainside).
Hence, the annual average ELAs (Fig. 3a) and winter precipitation (Fig. 3c) were rather strongly negatively correlated (r = −0.63).Like for CPDDs, a well-supported spline regression enhanced their similarities (r = −0.80).This showed that the irregularities in the increasing trend of the average ELA at the scale of the study region were linked to winter precipitation anomalies: excesses in 1995 and 2001, and deficits in 1990 and 1998, which respectively decreased/increased the annual average ELA.
Finally concerning Fig. 7, one has to keep in mind that (1) the highest correlation coefficient was found for summer CPDD; and (2) except for few glaciers in the Ecrins Massif, significant correlations were also found with winter precipitation but to a lesser extent than with summer CPDD.These points are in good agreement with results found at a glacier scale based on direct field measurements showing that the mass balance variability (and hence the ELA) is primarily controlled by summer ablation variability, which is in turn closely linked with summer CPDD (e.g.Martin, 1974, Vallon et al., 1998;Vincent, 2002).However, if Fig. 3 confirms that the observed increasing trend in the average regional ELA is in fact mainly driven by increasingly warm summer temperatures, it also shows that nonlinear low frequency fluctuations in the average regional ELA are significantly related to winter precipitation anomalies.
In terms of sensitivity of the ELA to climate variables, Fig. 8 shows the scatter plots of ELA anomalies vs. summer CPDD anomalies (Fig. 8a) and winter precipitation anomalies (Fig. 8b).For this figure, both the ELA anomalies and the summer CPDD and winter precipitation anomalies were computed from the annual average values of the 43 glaciers and the six weather stations with the highest correlation coefficients.These two graphs show that the sensitivity of the ELA to summer temperature was 115 m • C −1 (note that the length of the period used to compute the summer CPDD was 120 days in our case, and that the temperature gradient used to compute the CPDD at 3000 m a.s.l. was 6 • C km −1 ), and the sensitivity of the ELA to winter precipitation was 48 m/100 mm.This means that an additional amount of winter precipitation of about 240 mm is needed to compensate for a daily average increase in summer temperature of 1 • C. Note that in our estimate, we do not consider that the increase in summer temperature would lead to a lengthening of the ablation period.
Our estimate of ELA sensitivity to summer temperature, 115 m • C −1 , is in the middle range of values reported in the literature which range from 60-70 m • C −1 (Vincent, 2002) to 160 m • C −1 (Gerbaux et al., 2005).Other values from the European Alps and other mid-latitude temperate regions that can be found in the literature are 90-115 m • C −1 (Braithwaite and Zhang, 2000), 100 m • C −1 (Zemp et al., 2007), 125 m • C −1 (Kuhn, 1981) or 130 m • C −1 (Oerlemans and Hoogendoorn ,1989).Such a difference in values mainly results from the approaches used that vary from empirical regressions from ablation stake measurements on a monitored glacier and temperature data from one weather station located close to the glacier, to modeling approaches based on physical processes and providing information on the sensi-tivity of the ELA to all meteorological variables.Here, our approach is empirical, but relies on large number of glaciers and weather stations and can be considered as relevant at the scale of French Alpine glaciers.4, one can first see that the temporal/spatial decomposition is very good, with a low R 2 res equal to 0.25.This suggests a very clear distinction between the effects of space and time in the total data variability.Furthermore, the average temporal effect for each year and the interannual spatial effect for each glacier are well captured by the symmetric temporal and spatial regression models, with high determination coefficients in time (R 2 time = 0.73) and in space (R 2 space = 0.72), leading to a total determination coefficient (R 2 tot = 0.73).This confirms the capacity of the chosen explanatory variables to model the average spatial and temporal effects.This is particularly true for the temporal explanatory variables: with only two climate variables, the regression explains more than 30 % of total variance (main term), leading to the predominance of interannual fluctuations with regards to interglacier variability (R 2 time/space = 0.56).As a consequence, the average annual ELAs modeled from temporal explanatory variables follow annual fluctuations of the mean ELA (β t series) very well, and roughly match the spline adjustment presented in Fig. 3a.A Spearman test showed that CPDD and winter precipitations are independent, so that, since all the explanatory variables are standardized anomalies, the regression parameters given in Table 4 are directly the marginal correlation coefficient with the β t series, 0.42 and −0.31, respectively.
For the spatial term, the four morpho-topographic variables explain about R 2 topo = 25 % of total variance, and R 2 space = 72 % of the variability of the interannual average ELA for each glacier (α i series).For the spatial explanatory variables, a slight correlation exists between the surface area and the average elevation of the glaciers.However, the values of a k (Eq.5) for each explanatory variable remain close to the correlation values of the bivariate regressions (Sect.4.1.),suggesting that the model is not overparameterized so that the full inferred multivariate relationship can be interpreted with reasonable confidence.
Figures 9c and 10c (Figs.9d and 10d) present the annual and local determination coefficients (the annual and local bias).Since the bias is generally close to zero, and the determination coefficient acceptable, this confirms the overall capacity of the model to reproduce the data set variability, in accordance with the inferred low R 2 res = 0.25 value for each year/each glacier.In other words, the assumption of complete separability of the effects of space and time appears to be correct for most of the years/glaciers whose ELA time series can be reasonably approximated by the regional annual mean series/interannual local mean series, weighted by the interannual mean effect for each glacier/the regional annual effect for each year, respectively.However, for some years or glaciers, the high bias/low determination coefficients show that this is not true, i.e. that the spatio-temporal decomposi-tion is insufficient to fully represent the variability of ELA values measured using the remote sensing images.For these peculiar years or glaciers, certain points are worth emphasizing to explain the non-negligible spatio-temporal residuals: -In Fig. 9c, the years 1990 and 2008 showed low determination coefficients.This means that for these two years, the ELA variability between the glaciers concerned differed from that of the whole period.Indeed, we mentioned in Sect.4.1.that, over the whole study period, the average ELA of the glaciers of the southern sector of the study region was about 150 m higher in altitude than that of the glaciers located in the northern sector.In the years 1990 and 2008, the situation was not the same, the ELA was homogeneous in 2008, and, in 1990, the ELA of the glaciers located in the southern sector was lower than the ELA of the ones in the northern sector.This situation is not clearly explained by the  -In Figure 9d, the only year with an annual bias higher than ±5 m was the year 2002.It should be noted that the year 2002 followed a year with a very high winter accumulation leading to a low ELA.During the field campaigns at the end of summer 2002, a firn line (firn from the year 2001) and the snow line were observed.However, the difference between the firn line and the snow line is difficult to distinguish on images with 30 m resolution (most of the ELA in 2002 were measured on Landsat images).As a consequence, it is quite possible that for some of the glaciers studied, the firn line was delineated instead of the snow line.This used to be the limitation of the use of satellite images, however, with the new sensors launched in the recent years that have a resolution of few meters or even tens of centimeters, this limitation is now reduced.
-In Fig. 10c, some glaciers (e.g.Gébroulaz or Violettes glaciers) show low determination coefficient.This means that for these glaciers, the ELA interannual variability differed from the one computed for all the glaciers studied.For Gébroulaz Glacier, Rabatel et al. (2005) had already reported specific behavior.These authors assumed that the peculiar accumulation conditions on this glacier could lead to a low massbalance gradient.This could be one of the causes of the different interannual variability of the ELA on this glacier.For Violettes Glaciers (the only glacier with a local bias above 15 m, Fig. 10d), as already mentioned in Sect.4.2, the peculiar morpho-topographic conditions could be responsible for the different interannual variability of its ELA.

Conclusion
This study presents time series of ELA for 43 glaciers located in the western Alps over the 1984-2010 time period.ELAs were derived from the end-of-summer snow line altitude, an accurate proxy of the ELA for mid-latitude mountain glaciers.More than 120 satellite images were used to obtain an annual coverage of the 43 glaciers during the 27 yr period.
In parallel, this study analyzed changes in climate variables in terms of summer CPDD and winter precipitation.These two variables were obtained for the same time period from 22 weather stations located in and around the study area.From these time series several points are worth emphasizing: -Assuming a linear trend over the study period: the average ELA increased by about 170 m; the summer CPDD increased by about 150 PDD at 3000 m a.s.l., while winter precipitation remained rather stationary.
-Temporal variability of the ELAs, summer CPDD and winter precipitation were very homogeneous over the study region.Regarding spatial variability, homogeneity was high for summer CPDD, but glaciers in the southern part of the study area differed from the others.
Bivariate and relationships between the ELAs and morphotopographic and climate variables were also established.These relationships allow us to highlight the fact that: -The average ELA over the study period of each glacier is strongly controlled by its morpho-topographic variables, namely its average altitude, surface area, latitude and aspect.
-The interannual variability of the average ELA of the glaciers studied is strongly controlled by climate variables, with the observed increasing trend in the average ELA truly mainly driven by increasingly warm temperatures, even if nonlinear low frequency fluctuations in the average ELA time series appear to be significantly related to winter precipitation anomalies.
From these relationships, we propose an expansion of Lliboutry's approach, implemented within a hierarchical Bayesian framework, to reconstruct the ELA fluctuations of any glacier in the study area.Considering two climate variables and four morpho-topographics variables, we show that our approach is well suited to: (1) isolate an average temporal effect for each year and an interannual mean effect for each glacier; (2) model and explain their fluctuations from a reduced number of climatic and topographic explanatory variables whose respective weight and related uncertainties are quantified in a consistent manner; and (3) well approximate the full variability of our ELA data set for each year and for each glacier and identify specific spatio-temporal behaviors that do not match the assumption of complete separability of space and time effects.
Finally, both methods presented in this study, based on remote sensing and modeling, are powerful tools to provide annual series of ELA that will be useful in hydrological modeling of high mountain watersheds where there is a need to integrate the glacier component in such an approach and when no field data are available.

Fig. 1 .
Fig. 1.The 43 glaciers studied (red) among the 593 glaciers in the French Alpine glacier inventory (in blue) plus additional glaciers on the Italian and Swiss sides of the Franco-Italian and Franco-Swiss borders.The three boxes represent the main glacierized ranges quoted in the text and figures, i.e.Mont Blanc, Vanoise and Ecrins.

Fig. 2 .
Fig. 2. Example of Landsat image used to delineate the snow line (red line) with a combination of spectral bands: green, near-infrared and shortwave infrared.Argentière Glacier (Mont Blanc Massif), Landsat TM5 image acquired on 9 September 1987.

Fig. 3 .
Fig. 3. Changes over the 1984-2010 period in (A) the ELA for the 43 glaciers studied (see Fig. 1 andTable 1); (B) summer CPDD; and (C) winter precipitation recorded by the weather stations used (see Fig.1and Table1).In each graph, the horizontal black bar represents the annual average of the sample, the gray box represents the median interval (Q3 − Q1), and the vertical black lines show the interval between the first and the last decile (D1 and D9).On each graph (described in Sects.3.1.1.and 3.2.), the dashed line represents the smooth underlying trend captured by a cubic smoothing spline regression(Lavigne et al., 2012).The dotted line is the corresponding 95 % confidence interval.The yellow and purple boxes highlight the years in which the average ELA at the scale of the study region is linked with positive (yellow) or negative (purple) winter precipitation anomalies (see Sect. 4.2).

Fig. 4 .
Fig. 4. Average increasing rate of the ELA (m yr −1 ) over the 1984-2010 period of each of the glaciers studied with respect to the aspect of the glacier.The color of the diamond identifies the massif in which the glacier is located: blue = Ecrins, green = Vanoise, red = Mont Blanc.

Fig. 5 .
Fig. 5. Factor analysis of (A) the ELAs measured on satellite images; (B) the summer CPDD computed from weather stations; and (C) winter precipitation computed from weather stations.Detailed information about each graph is provided in the text.

Fig. 6 .
Fig. 6.Average ELA of each one of the 43 glaciers studied over the 1984-2010 period as a function of (A) glacier latitude, (B) glacier mean altitude, (C) glacier surface area, and (D) glacier aspect.The color of the diamonds in (A) and (D) identifies the massif in which the glacier is located: blue = Ecrins, green = Vanoise, red = Mont Blanc.

Fig. 7 .
Fig. 7. (A) Highest correlations among all the possible "glacier vs. station" pairs between each glacier (each symbol represents one glacier) and CPDD, plotted as a function of latitude.(B) The same but with the winter precipitation.

4. 3
Figures 9 and 10 summarize the results of the modeling approach from the symmetric temporal and spatial point of views.Hence, in each figure, A and B subplots illustrate the capacity of our method to model the common pattern of time (respectively the interannual local effect), getting rid of the interannual effect of each glacier (respectively the annual common effect).C and D subplots quantify the goodness of fit at the lower level of each annual/local series in terms of bias and determination coefficient.From Figs. 9 and 10 and Table4, one can first see that the temporal/spatial decomposition is very good, with a low R 2

Fig. 9 .
Fig. 9. Results of temporal modeling.(A) and (B) show the comparison between the ELA measured on the satellite images and the annual average ELA modeled using the temporal explanatory variables only, thus getting rid of the interannual mean specificity of each glacier.a o + β t (Eq.6) is the same plus "white noise" not explained by the temporal explanatory variables.(C) and (D) quantify the goodness of the fit for each annual series in terms of determination coefficient and bias, respectively.

Fig. 10 .
Fig. 10.Results of spatial modeling.(A) and (B) show the comparison between the ELA measured on the satellite images and the average ELA modeled using the spatial explanatory variables only, thus getting rid of the annual effect common to all glaciers.a i is the same plus "white noise" not explained by the spatial explanatory variables.(C) and (D) quantify the goodness of fit for each local series in terms of determination coefficient and bias, respectively.

Table 2 .
Characteristics of the wavelength of the spectral bands used to identify the snow line on the satellite images of the three sensors.

Table 3 .
Statistics of the three subsets used to reconstruct missing values.

Table 4 .
Posterior characteristics of the proposed model applied on the 43 ELA time series for the 1984-2010 period.a k and b k denote the regression parameters for each explanatory variable in brackets.Overall, 2.5 % and 97.5 % denote the lower and upper bound of the 95 % credible interval.