Changes in glacier Equilibrium-Line Altitude (ELA) in the western Alps over the 1984-2010: evaluation by remote sensing and modeling of the morpho-topographic and climate controls

The paper employs an analytical perspective that will probably be unfamiliar to most glaciologists. A reviewer prepared to evaluate the statistical analysis could probably suggest how the authors could make the paper more accessible to a broader range of readers. The paper uses factor analysis, which occurs only seldom in the glaciological literature, for instance, but devotes not a word to describing it.


Conclusions References
Tables Figures

Back Close
Full

Abstract
We present time series of equilibrium-line altitude (ELA) measured from the end-ofsummer snowline 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 parameters (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 parameters (elevation, aspect, latitude) highlighted the following: (1) the average ELA over the study period of each glacier is strongly controlled by morpho-topographic parameters; and (2) the interannual variability of the ELA is strongly controlled by climate parameters, 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 parameters, 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 130 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, 2011).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 glacier changes at the scale a mountain range and 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.Introduction

Conclusions References
Tables Figures

Back Close
Full Here, 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 snowline.Indeed, for mid-latitude mountain glaciers, the end-of-summer snowline 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 snowline 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 Sect. 2 we describe the study area, and the data and methodologies used in the current work.In Sect.3, we detail temporal and spatial changes in each one of the parameters: ELA, summer cumulative positive degree days and winter precipitation.
In Sect.4, we focus first on the relationship between the ELA and morpho-topograpic parameters, and second, on the relationship between the ELA and climate parameters.Finally, relying on the hierarchical framework proposed by Eckert et al. (2011) for mass balance point measurments, we incorporate these relationships in an expansion of Lliboutry's variance decomposition model (1974) to reconstruct the spatio-temporal variability of annual ELA time series.We compare and validate the modeled data with the ELA measured from satellite images.

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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 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 snowline every year during the study period; (2) glaciers with all aspects had to be represented; and (3) glaciers located in all the glaciated mountain ranges in the French Alps from the southern-most ( 44• 50 N) to the northern-most ( 46• 00 N) had to be included.Among the selected glaciers, three belong to the GLACIOCLIM observatory (http: //www-lgge.ujf-grenoble.fr/ServiceObs/index.htm) 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 three glaciers were used to assess the representativity 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 snowline 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 completly 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 Supplement: Table 1 for a detailed list of the images used).Introduction

Conclusions References
Tables Figures

Back Close
Full For the snowline delineation on multispectral images, a test of band combinations and band ratios to facilitate the identification of the snowline 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 short-wave infrared bands (see Table 2 for details about the wavelengths of each satellite) was the most appropriate to identify the limit between snow and ice.The snowline was delineated manually on each glacier.

DEM, computation of the altitude of the snowline and of the glacier morpho-topographic parameters
To avoid border effects on the glacier banks (Fig. 2), which could mean that the position of the snowline depended on local conditions (shadows from surrounding slopes, additional snow input due to avalanches or wind drift), the SLA was only measured on the central part of each glacier.The average altitude of each snowline was calculated by overlaying the shapefiles containing the digitized snowlines on a DEM.For each glacier, the DEM was also used to compute the following morpho-topographic parameters used in this study: surface area, mean elevation, aspect, latitude and longitude of the glacier centroid.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.Consequently, either one of the DEMs can be used to compute the altitude of the snowline without influencing the results.
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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 SLA 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 i is the glacier, t is the year, β t is a term depending on the year only which is common to all the glaciers analyzed, α i is the interannual mean for each glacier, 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 Introduction

Conclusions References
Tables Figures

Back Close
Full terms, but to estimate each missing value SLA i t 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 "hat" 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 dataset 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 ε i t , representing the model errors was 73 m.This is less than half the RMSE of the measured SLAs (174 m), showing that the use of the bi-linear model makes sense on this dataset.

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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 glaciated areas (except Lyon-Bron weather station which was selected as a regional reference).Among the 40 weather stations, only 20 which had continuous series of winter precipitation and summer temperatures over the period 1984-2010 were finally used.Additional weather stations in Switzerland were also used.The climate parameters 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 snowline.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).

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

Changes in the ELA
As mentioned above, it has been demonstrated that the end-of-summer SLA is representative of the ELA and is highly correlated with the mass balance in mid-latitude glaciers (Rabatel et al., 2005(Rabatel et al., , 2008)).Figure 3  glaciers in the French Alps for which field mass balance measurements are available for the 1984-2010 period.Because of this strong correlation and for the sake of simplicity, we hereafter only use the term ELA.

Temporal variability of the ELA
Figure 4a 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 linear trend of 6.4 m yr −1 , 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 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 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.Introduction

Conclusions References
Tables Figures

Back Close
Full Considering each glacier individually, Fig. 5 shows the rate of increase in the ELA over the whole period for each 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 in 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 ELA
Figure 6a 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 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. 7a) which could be associated with differences in winter accumulation between the southern and northern Alps (see below) and warmer Introduction

Conclusions References
Tables Figures

Back Close
Full temperature in the southern sector.Indeed, Fig. 7a 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 4b shows the changes in summer CPDD over the 1984-2010 period for the 22 weather stations used.Summer CPDD present an increasing linear trend, averaging 5.3±1.9CPDD yr −1 at 3000 m a.s.l.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).
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 m higher than the average for the whole period, and is clearly associated with an unusually warm summer 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. 6b) 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 linear trend was observed over the study period (Fig. 4c).However, one can note that, on average, interannual variability was lower after 2001 (standard deviation divided by 2 after 2001), 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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. 6c).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 morphotopographic parameters 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 parameter.Finally, we describe and exploit an expansion of Lliboutry's approach (1974) that explicitely incorporates both climate and morpho-topographic parameters to reconstruct annual ELAs.

The role of morpho-topographic parameters
Figure 7 shows the relationship between the average ELA of each glacier computed over the whole study period and the latitude of the glacier (Fig. 7a), the mean altitude of the glacier (Fig. 7b), the glacier surface area (Fig. 7c), and the glacier aspect (Fig. 7d).
We already mentioned the relationship between the ELA and the latitude (Fig. 7a) in Sect.3.1.2,where we noted that the average ELA of the glaciers located in the southern Introduction

Conclusions References
Tables Figures

Back Close
Full part of the study area was about 150 m higher in altitude that 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 7b and c shows 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, meaning 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.As a consequence, the lower the ELA, the wider the accumulation zone, the bigger the glacier, the lower its snout and finally, the lower the mean altitude of the glacier.Current ELAs can certainly 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 can be deemed to be representative of its average ELA.
Figure 7d shows the average ELA of each glacier plotted as a function of the glacier aspect.For this morpho-topographic parameter, 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.Introduction

Conclusions References
Tables Figures

Back Close
Full

Climatic control of ELA interannual variability
Analysis of the variances and correlations between the ELA of each glacier and the climate parameters (summer CPDD and winter precipitation) was conducted for the study period.Considering our 27 yr dataset, 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 8 shows the correlation coefficients for the two climate parameters 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. 8a shows: (1) a generally high and always significant 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. 4 which shows strong similarities between annual variability of ELAs (Fig. 4a) and CPDD (Fig. 4b), 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 long-term variations in ELA and CPDD over the study area.
Regarding winter precipitation, Fig. 8b shows the particular behavior of the Ecrins Massif (blue triangles in Fig. 8b), 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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. over-accumulation or erosion (Girose and Mont de Lans are located on a dome and Violettes in a narrow and cashed mountainside).Hence, the annual average ELAs (Fig. 4a) and winter precipitation (Fig. 4b) 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. 8, 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. 4 confirms that the observed increasing trend in the average regional ELA is in fact mainly driven by increasingly warm 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 parameters, Fig. 9 shows the scatter plots of ELA anomalies vs. summer CPDD anomalies (Fig. 9a) and winter precipitation anomalies (Fig. 9b).For this figure, both the ELA anomalies and the summer Introduction

Conclusions References
Tables Figures

Back Close
Full 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 CPDD 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 additionnal 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 Gerbaux et al., 2005).Such a difference in values mainly results from the approach 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 sensitivity of the ELA to all meteorological parameters.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.

ELA modeling from morpho-topographic and climate parameters
From our dataset of ELA series, glacier morpho-topographic parameters, and climate parameters, we propose a simple expansion of Lliboutry's approach (1974) to model and reconstruct the ELA fluctuations of any glacier in our study area over the study period.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 Introduction

Conclusions References
Tables Figures

Back Close
Full interest, more specifically: the regression parameters relating ELA fluctuations, their morpho-topographic and climate drivers, and their respective weight.
In detail, starting from the available ELA 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 covariates considered for each glacier, i.e. the morpho-topographic parameters: 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 covariates considered.
Similarly, the temporal term is decomposed as: where the X kt are the temporal covariates considered, i.e. the climate parameters 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 covariates considered.
For the sake of simplicity, the model is fed with reduced standardized variables, which 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 covariates (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

TCD Introduction Conclusions References
Tables Figures

Back Close
Full 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 priors were used for all parameters.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 ε i t 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 terms in ELA variability: 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: 10 All these statistics quantify the weight of mean spatial and temporal effects in the total variability, and the capacity of the chosen covariates to model these mean effects.At the "lower" level of the model, local and annual adjutment statistics can also be TCD Introduction

Conclusions References
Tables Figures

Back Close
Full computed as: respectively 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 sequence, 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, see Figs. 10 and 11).This method explains why the sum R2 topo + R2 clim + R2 topo res + R2 clim res + R2 res of posterior estimates (posterior mean) in Table 4 is very close but not fully equal to 1.  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 (R2 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 covariates to model the average spatial and temporal effects.This is particularly true for the temporal covariates: with only two climate parameters, 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 covariates follow annual fluctuations of the mean ELA (β t series) very well, and roughly match the spline adjustment presented in Fig. 4a.A Spearman test showed that CPDD and winter precipitations are independent, so that, since all the covariates 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 parameters 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 covariates, 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 covariate remain close to the correlation values of the bivariate regressions (Sect.4.1.),suggesting that the model is not over-parameterized so that the full inferred multivariate relationship can be interpreted with reasonable confidence.
Figures 10c and 11c  of the effects of space and time appears to be correct for most of the years/glaciers whose ELA time series can be reasonably approximed 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 decomposition 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 spatiotemporal residuals: -In Fig. 10c, 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 climate parameters and further investigations are needed to explain what happened during these two years.
-In Fig. 10d, 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 campains at the end of summer 2002, a firn line (firn from the year 2001) and the snowline were observed.However, the difference between the firn line and the snowline 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 snowline.This used Introduction

Conclusions References
Tables Figures

Back Close
Full 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. 11c, 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 mass-balance 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. 11d), as already mentioned in Sect.4.2., the peculiar morpho-topographic conditions could be responsible for the different interannual variability of its ELA.

Conclusions
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 snowline 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 parameters in terms of summer CPDD and winter precipitation.These two parameters 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.Introduction

Conclusions References
Tables Figures

Back Close
Full -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 morpho-topographic and climate parameters 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 parameters, 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 parameters, 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 parameters and four morphotopographics parameters, 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 covariates whose respective weight and related uncertainties are quantified in a consistent manner; and (3) reproduce the full variability of our ELA dataset 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 Introduction

Conclusions References
Tables Figures

Back Close
Full  Full  Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | illustrates this correlation for three Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Figures 10 and 11 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.Discussion Paper | Discussion Paper | Discussion Paper | From Figs. 10 and 11 and Table (Figs. 10d and 11d) 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 dataset variability, in accordance with the inferred low R Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

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 covariate in brackets.2.5 % and 97.5 % denote the lower and upper bound of the 95 % credible interval.