Variations of Lys Glacier (Monte Rosa Massif, Italy) from the Little Ice Age to the Present from Historical and Remote Sensing Datasets

Alpine glaciers respond to climate imbalance by adjusting their mass and length. In turn, these changes modify the glacial and periglacial environment, leading to increased supraglacial debris cover, the development of glacial lakes and glacier fragmentation. In this research, we investigated the evolution of Lys Glacier (Monte Rosa Group), by studying length, area and volume changes, and evolution of its supraglacial debris cover and proglacial lakes by means of historical sources and high-resolution aerial and satellite orthophotos. Lys Glacier retreated almost continuously, by nearly 2 km, from its maximum Little Ice Age position. More recently, the glacier lost 11.91% of its area between 1975 and 2014 and underwent fragmentation in 2009. Over the same period, glacier fragmentation and tongue stagnation affected the formation and rapid growth of a series of ice-contact lakes and led to a non-linear debris cover evolution. The glacier was also subjected to strong volume losses, with more than 135 m thinning on the ablation tongue from 1991 to 2014. Analysis of the meteorological records (1927–present) from the closest weather station reveals a considerable increase in average annual temperatures by more than 1°C from the mean of 1971–1989 to the mean of 1990–2017.


Introduction
Observations show that the strongest influence of global climate change is recorded in alpine environments and glaciers are the first to be affected In the current context of temperature warming and glacier regression, long-term data reporting changes in glacier properties represent an important asset, as terminus, area and volume fluctuations can provide important information concerning the response of a glacier to climate imbalance. These data are a fundamental input for models that enable reconstructing past climate (from terminus fluctuations [18,19]) and projecting future glacier changes [20], the availability of meltwater for domestic use and the production of hydroelectric energy [21,22], as well as the formation of future lakes and potential hazards [23] and the impact of glacier change on tourism [24].
In this study, to describe the evolution of Lys Glacier, a multiple approach was followed. Terminus fluctuations since the early nineteenth century were analyzed using a variety of sources including detailed bulletins with reports of glaciological campaigns; over a more recent period , the area changes of the glacier were estimated by using remote sensing datasets, i.e. satellite and aerial orthophotos, while volume changes were evaluated by comparing a pair of digital elevation models (DEMs) obtained from cartography and satellite images. In addition, these sources permitted us to evaluate the evolution (i.e. surface cover and patterns) of the supraglacial debris cover and of the ice-contact lakes over the same period of observation and gain insights into the geomorphological evolution of the glacier.

Study site
Lys Glacier drains the southern flank of the Mont Rosa Group (45°54' N, 07°50′ E) (Figure 1). The most recent Italian glacier inventory ( [14], data from 2005) reports a glacier area of 9.58 km 2 , a south westerly aspect, an elevation range between 2392 and 4323 m a.s.l. and a length of 5.71 km. While Lys Glacier is presently the fourth largest Italian glacier, comparatively few studies have been conducted to investigate its evolution over the past century and recent decades. Strada [25] described terminus fluctuations until the early 1980s; Pelfini et al. [26] estimated the glacier response time using dendrochronology; Rota et al. [27] used DEMs from cartographic sources to calculate volume changes between 1925 and 1994. None of these studies or iconographic sources show evidence of continuous supraglacial debris cover until the late 1980s. Since then, debris cover has been present on the lower sector of the ablation tongue, initially as the continuation of medial moraines at an elevation below 2550 m a.s.l. The medial moraines formed below the glacier icefalls (Figure 2), thanks to the debris supplied by the surrounding deglaciated rock walls (owing to macrogelivation processes, permafrost degradation and structural rock falls). In more recent years, debris has come to cover the entire glacier tongue below 2550 m a.s.l., until its detachment from the upper sector of the glacier, which occurred in 2009 [28]. The occurrence of debris cover is thus a recent phenomenon and probably a consequence of the present deglaciation, which has affected Lys Glacier as well as the other alpine glaciers [1,15]. Although Lys Glacier cannot be considered a debris-covered glacier even in recent times (as only a small sector of the glacier surface is continuously debris-covered, see [29]), the supraglacial morphologies and the processes affecting the debris-covered glacier tongue are comparable to those of larger debris-covered glaciers in Asia and in the Alps, including thermokarst processes and the formation of cavities and ice-contact lakes [29].
Glaciers and the Polar Environment 4

Terminus fluctuations
In this study, the terminus fluctuations of Lys Glacier were retrieved from different sources; early data concerning the nineteenth century are the same used in Strada [25] and references therein, i.e. early cartographic sources (map of the states of Sardinia from 1818 and cadastral maps), photographs and descriptions, particularly those in Monterin [30]; from the early twentieth century, length changes were obtained by extracting the field measurements available from the journals published by the Italian Glaciological Committee (CGI) [31,32]. Surveys of the glacier terminus positions consist of tape measurements from fixed  reference points in the glacier forefield, made at or near the end of the balance year, i.e. in the autumn season, with an estimated accuracy about 0.5 m, which can become worse in case of bad environmental conditions, poorly documented switches to new reference points, measurements taken over very long distances or at points on the snout outside the main flow line, or residual snow patches [12]. While the actual uncertainty from these issues is difficult to quantify, glacier terminus fluctuations remain an important asset to assess global climate and environmental change [19].
We followed recommendations by Citterio et al. [12] to perform qualitychecking of the measurements extracted from CGI journals, to minimize the issues related to changing reference points or surveyors over the years. When multiple measurements were available from the same survey, an averaged variation was calculated. In case of missing years in the record, periods up to 5 years where the reference point had remained the same were filled by uniformly distributing the total terminus variation over the missing years.
In case of a gap associated with an undocumented change in a reference point, the gap was not filled, and the first measurement following the gap was considered as the starting point for comparison with the measurement resulting for the following year, as in Citterio et al. [12]. The data record for Lys Glacier from CGI journals is rather uninterrupted and permits analyzing the glacier history in detail over the last century.

Data sources and preprocessing
Maps, aerial orthophotos and satellite images were analyzed to calculate the area and volume changes of Lys Glacier and changes in the area of the proglacial lakes. The available cartographic sources include large-scale maps from the Val d'Aosta region, produced in 1975 and 1991 at a nominal scale of 1:10000 ( system based on the ED50 datum; contour lines and elevation points were digitized as shapefiles, and a DEM was produced from the 1991 map with a pixel spacing of 10 m using ArcMap Topo to Raster utility. Aerial orthophotos from 1988 to 2012 (Table 1) were obtained from Geoportale Nazionale (http://www.pcn.minambiente.it/mattm/); these have a pixel size of 0.5 m and were downloaded in the UTM32 coordinate system based on the WGS84 datum. Most images were acquired in late summer months (August and September; Table 1), with minimum snow cover, while the orthophotos from 1988 and 2012 were acquired in early summer (June and July) and show evidence of residual snow, which however does not affect the glacier tongue. In addition to the aerial orthophotos, we obtained a stereo pair acquired from the Pleiades satellite constellation from the European Space Agency. The stereo pair was imaged on 1 September 2014 and provided at level 1B processing stage; we further processed it to obtain a DEM and an orthorectified image using NASA's AMES stereo pipeline (ASP). Eleven ground control points were selected from those available from Val d'Aosta regional authority (http://geonavsct.partout.it/pub/GeoNavITG/monografie.asp) to improve the geolocation accuracy of the DEM through bundle adjustment in the software; the DEM was then produced using the semi-global matching algorithm available in ASP [33], and the raw multispectral images were projected onto the DEM for orthorectification. Both the Pleaides DEM and multispectral image have a final resolution of 2 m.

Calculation of area and volume changes
Based on the available planar data ( Table 1), the glacier and proglacial lake outlines were estimated over multiple years: in the 1975 and 1991 technical maps, they were already drawn on the map, while in the other datasets, they were manually digitized as part of this study. The spatial resolution of the orthophotos allowed us to clearly identify the outlines and to distinguish the debris-covered parts of the tongue from the proglacial areas. Before digitization, all available data were reprojected to the UTM32N coordinate system based on the WGS84 datum for consistency. We further evaluated the accuracy of the manual delineation using the buffer method proposed by Paul et al. [34], by allowing the glacier and lake outlines to grow and shrink with a buffer of 2 pixels.
The analysis of elevation and volume changes was based on the comparison of the DEMs from 1991 and 2014. Before performing this comparison, the DEMs were co-registered using the approach developed by Berthier et al. [35] and applied by Fugazza et al. [36]. In this approach, one DEM ('slave') is iteratively shifted with respect to a reference DEM ('master'), to minimize the standard deviation of elevation differences over stable areas located outside the glacier (σ Δh ). We selected the oldest DEM as the reference and resampled both DEMs to a common resolution of 10 m. By applying the co-registration, we obtained a residual σ Δh of 4.18 m; the average elevation difference over stable areas was 7 m, which was subtracted from the reference DEM. Based on the elevation difference on the glacier surface, we then calculated the volume change as: where k = 1 … K are the pixels of the reference surface, A k is the area of each pixel and d k is the pixelwise difference between the 1991 and 2014 DEMs.
To express the uncertainty of volume changes, we used the approach described by Fischer et al. [37]: the uncertainty in elevation change ( σ ∆h ) is not considered as entirely correlated but scaled to account for the effective correlated area A cor , based on the correlation length (i.e. the distance at which two pixels are correlated). We calculated the correlation length using an empirical semivariogram in R software as 200 m and the effective correlated area as 0.12 km 2 . The uncertainty in volume change σ ∆h_cor is then expressed as: where A 1991 is the glacier area in 1991. The uncertainty calculation can in our case only be applied to the glacier tongue, as at altitudes above the glacier equilibrium line the presence of snow with little contrast deteriorates DEM reconstruction and leads to larger errors [38], particularly with the oldest DEM from 1991.

Mapping debris cover and its changes
To characterize the evolution of debris cover on Lys Glacier, we adopted the methodology proposed by Azzoni et al. [9] for glaciers in the Ortles-Cevedale Group. Supraglacial debris was mapped by employing a maximum likelihood (ML) supervised classification approach. Four classes were included in the classification: debris, ice, snow and shadow. Although we did not aim at mapping the other surface types, a higher number of classes was chosen to permit an improved mapping of debris cover; the shadow class was also included as some images are affected by topographic shadow; this occurred mostly in the accumulation basins except for the orthophoto from 2003, where shadows occur over the whole glacier; thus, this orthophoto was excluded from the analysis ( Table 1). We also excluded the maps from 1975 and 1991, as the information on the presence of debris was not available. For each of the satellite/aerial orthophotos, we independently selected 10 to 15 training areas. To estimate the accuracy of the classification approach, we performed a validation test on the 2006 and 2012 orthophotos, by selecting 100 random points for each test and comparing the results of ML classification against those from a manual classification. The 2 years were chosen because they represent best (2006) and worst (2012) conditions for debris cover mapping, in terms of image quality, presence of snow and shadows. For both years, we calculated overall accuracy (the ratio of the number of correctly classified points to the number of total points), as well as producer's accuracy (PA) and user's accuracy (UA), for the different classes. PA denotes the ratio of the number of correctly classified points to the total (manually classified) point for a class and is the complement of omission errors; UA denotes the ratio of the number of correctly classified points to the number of predicted points for a class and is the complement of commission errors.

Meteorological data
In order to interpret the glaciological data of the study area, we analyzed meteorological variables recorded by the weather station installed in 1927 at Gressoney d'Ejola, 1850 m a.s.l., at a distance of 3.7 km from Lys Glacier (Figure 1). This long dataset allows investigating the recent climate behaviour, which influences the evolution of the glacier. At this weather station, manual observations of air temperature, liquid precipitation, total snow depth and thickness of fresh snow, atmospheric pressure, relative humidity, wind speed and velocity and cloud type and cover were collected three times a day (at 8 am, 2 pm and 7 pm) from 1927 to 2012. The data record is almost uninterrupted except for a short period from 1962 to 1970, when the station was temporarily moved downvalley, at 1730 m a.s.l. in the Orsia village. In 2002, an automatic weather station was installed at the same site, and it replaced the manual weather station after its dismissal in 2012. The automatic weather station measures air temperature, total precipitation (by a heated gauge) and snow height (by an automatic webcam) every hour. The contemporary presence of manual and automatic instruments has allowed ensuring the homogeneity and continuity of the series. Although a weather station also exists at a higher elevation and closer proximity to Lys Glacier (Alpe Courtlys, 1992 m a.s.l., 2.5 km distance), this station was installed in 2001 and does not permit a long-term analysis of climate variables. In this study, we calculated monthly, seasonal and annual averages of air temperatures and monthly, seasonal and annual cumulated liquid/ solid precipitation for the analysis of climatological conditions.

Terminus fluctuations
The history of the advances and retreats of Lys Glacier is summarized in Figure 3. We here report the complete record of data but also focus on the most recent period

Area and volume changes
Over the period of observation, Lys Glacier underwent large changes in area. From 1975 to 1999, little change is evident on the western section of the glacier tongue (Figure 4), which reached the lowest elevations (ca. 2350 m a.s.l.). However, the easternmost part of the glacier tongue had already started retreating by 1999. Far greater changes than in the earlier period occurred from 1999 to 2014, with marked regression both in the western and eastern sections of the glacier tongue (Figure 4a). On the western section, the glacier tongue detached from the upper portion of the glacier at the base of an icefall, whereas on the eastern section the retreat exposed steep rocks underneath the ice. In 2014, the two main branches of the glacier still maintained a small connection at the base of a large rock outcrop (Figure 4a) and in the accumulation basins. Over the years, changes in Lys Glacier area were almost linear in rate, except for the first year, 1975, when the glacier had its largest area (11.3 km 2 ) of the period of observation (Figure 5). Since 1988, the area decreased at a rate of −0.045 km 2 y −1 . In view of the high spatial resolution of the imagery, the uncertainty in the glacier outlines is below 2%, which supports the evidence of the general trend. Beside changes in the glacier area, recent years have also seen a transformation of the geomorphological setting of the glacier outwash plain. While the exact date marking the origin of glacier lakes is difficult to establish, by 1999 five small ice-contact lakes had formed on the glacier tongue. Progressive stagnation in the ablation area allowed the lakes to grow and new ones to develop. Thus, from five scattered ponds in 1999 (totalling 1527 m 2 ) a larger overall area is seen together with the coalescence in two larger lakes in 2006 (8730 m 2 in total). In 2014, three further lakes formed while the existing ones grew in size (Figure 4b, c). The area of glacier lakes increased almost exponentially (summing the area of all lakes) from 1999, reaching 58,560 m 2 in 2014 ( Figure 5). The uncertainty in the lake outlines is larger than that of the glacier outlines, up to 30% in 1999 because of their relatively small size.    On the eastern part of the ablation tongue, the decrease in ice thickness reached instead 41 m. Aside from the glacier tongue, the other areas of the glacier show a more complex pattern, with areas of both positive and negative changes occurring. A loss in glacier thickness is seen particularly in the lower sectors of the western glacier tongue and in a narrow band in the uppermost reaches of the western accumulation basin, as well as on the eastern tongue in its lower easternmost area (Figure 6). A positive difference from 1991 to 2014 can instead be seen particularly on the western accumulation basin, where an apparent increase in thickness is observed between 10 and 20 m, and a maximum of above +75 m was recorded. Considering the entire glacier, the ice loss signal is still predominant, with an average of −4.34 m and a volume change of −47.06 × 10 6 m 3 .

Debris cover evolution
On Lys Glacier, debris initially increased from 1988 to 2000, when it reached almost 1.3 km 2 on the glacier tongue (Figures 5 and 7a). In 1988, debris was present in narrow medial moraines and more abundantly at the terminus (Figure 2). In 2000, a more homogeneous coverage of the glacier tongue can be seen (Figure 7b), as a mantle of debris appears covering the glacier at the higher elevations at the margins of the western tongue, also suggesting a possible input from the lateral valley walls; conversely, on the eastern tongue and parts of the eastern accumulation basin, coverage is sparse and patchy. Since 2000, the expansion of supraglacial debris appears to have halted as the total area decreased ( Figure 5): this effect was probably caused by the shrinking of the stagnating glacier tongue, which was however entirely debris-covered by 2005 (Figure 7c and d). The spread of supraglacial debris appears to have slowed down also at the higher elevations, although limited evidence for increasing coverage is seen for the eastern tongue of the glacier. The relatively high slope and the consequent presence of seracs might limit debris accumulation in those areas.
The accuracy of debris cover maps was evaluated separately for 2006 and 2012, as seen in Tables 2 and 3. In both years, PA and UA are very high for debris: both are 100% in 2006, while PA is 89% and UA 97% for 2012. Overall accuracy was 83% in 2006 and 71% in 2012. The limited overall accuracy in both tests is mostly caused by the difficulty in distinguishing between snow and ice (Figure 7b, where no ice was identified), which however is not crucial for the analysis of debris cover. We estimate the accuracy of debris cover mapping in other years to lie between the two values reported for 2006 and 2012; however, it is also possible that the debris cover

Climatological analysis
Due to its high altitude and position, Lys valley experiences cold winters and temperate summers. Heavy rainfall can occur when south humid Mediterranean winds blow and collide against the orographic barrier of the southern slopes of Monte Rosa. Perturbations from the west and northwest are more frequent, but they discharge their rain/snow content mainly on the Mont Blanc and Valais areas, leaving the northeastern extremity of the Aosta Valley almost dry.
From 1928 to 2018 (excluding the period 1962-1970 when the station was temporarily moved downvalley), the mean annual temperature at Gressoney d'Ejola station was +4.4°C ranging from +2.6°C in 1984 to +6.1°C in 2015. Since we cannot focus on the commonly used 30-year reference period (1961-1990), we considered 1952-1961 and 1971-1990. Over this period of observation, the mean annual temperature was +4.0°C, slightly lower than the average of the past 30 years (1989-2018), which was +4.6°C. Climate warming is more evident when summer temperatures are compared: from +11. 6°C during 1952-1961/1971-1990 to +12.7°C during 1989-2018, with a mean of +12.2°C over the whole period. The hottest month is July with a mean temperature of +13.2°C followed by August (+12.7°C) and June (+10.7°C). Generally, every 3 years the maximum daily of +25°C is recorded, even if the absolute maximum (+28.2°C) was observed on 4 August 2017 and 11 August 2003; 2003 was the hottest summer on record (average + 15.4°C; Figure 8) with relatively low amounts of liquid precipitation (216.6 mm corresponding to 74% of the mean summer total, 293.6 mm). Conversely, the coldest summer was 1977 (+9.9°C) with very heavy rainfall (428.6 mm corresponding to 146% of the mean summer total).  Table 3.

Accuracy estimation for debris cover classification based on the aerial orthophoto from 2012.
During the winter seasons, the monthly mean temperature is −2.6°C in December, −3.6°C in January and − 2.8°C in February with an absolute minimum up to −25.0°C recorded on 10 February 1986. The mean winter temperature is −3.0°C, and the coldest season was 2009-2010 with an average temperature of −5.6°C (and relatively low precipitation: 174.7 mm of total precipitation, 67.5 cm of snow depth, 181 cm of fresh snow; Figure 8), while the warmest one was 1948-1949 (0.0°C), also characterized by the lowest total precipitation (48.4 mm corresponding to 24% of the total mean winter amount, 201.4 mm; Figure 8), the lowest mean snow depth and cumulative fresh snow (7.7 cm and 52 cm corresponding to 12% of the mean winter amount-63.9 cm-and 28% of the mean winter total fresh snow, 188.7 cm, respectively) and the lowest number of days with snow cover (58 days corresponding to 67% of the mean total winter days, 86.8 days).
Generally, the coldest day is on 5 January (−4.5°C on average). Frost days (T max < 0°C) generally occur from October to April, even if days with T min < 0°C can occur even in July. Thaw (T ave > 0°C) begins at the end of June.
Comparing the two 30-year periods, the mean annual cumulated precipitation  Considering the whole year (average of 1927-2018), 1134.6 mm of rain and melted snow are distributed in 111 rainy and/or snowy days. Total precipitation is lower than the amount falling in the areas downvalley, which are even more exposed to the southerly humid winds, but much higher than the amount received in the dry Aosta Valley where 500 mm per year are generally recorded owing to its intra-alpine position. From November to April, there is an average 371.7 cm of fresh snowfall each year, almost equally divided between the various months: 61.9 cm per month ranging from 54.9 cm in November to 65.9 in December. Heavy snowfalls occur with Mediterranean temperate humid winds, with the daily maximum of 120 cm recorded on 1 January 1986 and the monthly maximum of 250 cm in April 1989.

Discussion
Since the LIA, and particularly over recent times, Lys Glacier has undergone remarkable changes in all the investigated parameters of length, area, volume and debris cover. The terminus fluctuation curve for Lys Glacier is strikingly similar to that of other published curves for Mer de Glace (Mont Blanc region, France) and Unterer Grindelwald (Bernese Oberland, Switzerland), which share some of the longest records of terminus fluctuations for alpine glaciers [39,40], suggesting that in spite of the distance between these glaciers, the climatic setting with a predominant influence from westerly winds is similar. All three glaciers share two distinct advance phases during the LIA, interrupted by a period of retreat, which according to Vincent et al. [41] was caused by a decrease in winter precipitation; small differences however exist in the timing and magnitude of such advances: the maximum length of Lys Glacier in the past two centuries was reached in 1821, similarly to Mer de Glace (and Rosenlauigletscher, [40]), while the extent of Unterer Grindelwald and most other alpine glaciers peaked around the 1850s [40]. Well documented for all three glaciers is also the rapid retreat following the end of the LIA around 1860, although Mer de Glace also shows a small readvance around 1867, for which no evidence exists for Lys Glacier. Both Lys and Mer de Glace then enter a period of relative stability until the 1930s (although marked by several small advances and retreats); a period of marked retreat follows, attributed to enhanced solar radiation [42] lasting until the short-lived advance phase of the late 1960s/early 1970s, which lasted longer, up to 1995, for Mer de Glace compared to Lys and which is generally observed for mostly alpine glaciers [12], although smaller glaciers tend to have larger readvance periods in the twentieth century, indicating a shorter reaction time [43].
At Gressoney d'Ejola, the late 1970s appear particularly favorable years for glacialism, with cool summers and high amounts of winter snowfall. Since the 1980s, a clear warming trend has emerged for summer temperatures, which are 1.3°C higher from the mean of 1971-1989 to the mean of 1990-2017 (Figure 8) and 1.1°C higher when including 1952-1961/1971-1989 (+12.7°C in 1990-2017 compared to +11.6°C), while no clear signal can be seen in winter temperatures, total and solid precipitation. This is in line with trends previously observed for Italy [44] and high elevation regions [45] and explains the large retreat rates seen since 1985 (−13.85 my −1 in 1985-2017). Considering a longer set of temperature data analyzed for the Alps (1856-1998), the whole twentieth century was characterized by rising temperatures, at a rate of 0.50°C per century (considering summer temperatures [46]), even before the record-breaking decades of the 2000s and 2010s.
As concerns glacier area, the values reported here are always larger than those of the recent Italian glacier inventory (9.58 km 2 in 2005 [14]), even in 2014, because their outlines did not take into account a small area in the eastern part of the glacier, which was here included for consistency with the 1975 and 1991 outlines drawn on the maps. The rate of change of Lys Glacier found in this study (−0.44 kmy −1 or − 0.4% y −1 ) is however comparable to that of alpine glaciers, both in Italy and in the other alpine countries [14,46,47]. The increase in debris cover is lower than that reported by Azzoni et al. [9] for the Ortles-Cevedale region, where 38 glaciers were reported to have on average a 13.3% higher proportion of their area covered in debris in 2012 (reaching 30%) than 2003, while Lys Glacier went from 7.9% in 1988 to 12.4% in 2000, decreasing again to 8.5% in 2014. These values are also lower than those reported by Shukla et al. [5] for Samudratapu Glacier in Indian Himalayas, where debris cover nearly doubled over less than 3 years and in line with observations of glaciers in Caucasus by Stokes et al. [11], who describe a 3-6% increase in debris cover and a 57% increase in supraglacial and proglacial lake area. The formation of lakes on stagnating debris-covered glacier tongues was also observed by Kirkbride and Warren [48] and is attributed to the presence of ice-cored moraine which prevents meltwater runoff and favors the accumulation of water in depressions left by melting ice. Similarly to Stokes et al. [11], we also found that debris cover has not halted glacier retreat, counter to the evidence that a thick debris cover is known to reduce ablation (see, e.g. [49]). A field campaign conducted in 2006 revealed that debris thickness is generally above 10 cm and up to 60 cm, well above the critical threshold for which the insulating effect prevails on the albedo effect [49]. Thus, while it is possible that mass wasting would have been even higher without debris, it is more likely that glacier retreat occurred owing to the presence of ice-contact lakes and cavities, which enhance melt through backwasting [50].
Unlike retreat rates, mean thickness and volume changes for Lys Glacier are noticeably lower than in similar studies conducted on alpine glaciers: D'Agata et al. [21] report a decrease in ice thickness of −14.91 m for glaciers in Sondrio Province, Central Italian Alps, from 1981 to 2007, corresponding to −0.57 my −1 , while we estimated thinning of Lys Glacier to be 0.19 my −1 . In their study of all glaciers in the Swiss Alps, Fischer et al. [37] report an area-weighted mass balance of −0.62 m w.e. y −1 from 1980 to 2010, while the geodetic mass balance of Lys Glacier (using a conversion factor of 0.85 accounting for the average density of ice and firn as done by [37]) would be −0.16 m w.e. y −1 , which is at the low end of the scale for glaciers analyzed in the study by Fischer et al. [37], although a few Swiss glaciers do share a similar geodetic mass balance. The geodetic mass balance of Lys Glacier would also be lower than that reported by Berthier et al. [51], i.e. -1.05 ± 0.37 m w.e. y −1 for glaciers in the Mont Blanc massif using ASTER, SPOT and Pleiades DEMs between 2000 and 2014. If we exclude the possibility that debris cover has slowed down thinning for reasons stated above, the large differences and the apparent low average elevation changes of Lys Glacier can be explained by (1) the relatively large size of the glacier accumulation basins and the presence of seasonal snow on the Pleiades image (Figure 2a) and (2) interpolation errors, especially in the oldest cartographic DEM, also seen particularly affecting the accumulation basins of the glacier where little change is expected to occur over the years [38], while we observe areas with apparent thickening of 20 m and up to 75 m. Considering the glacier tongue, the glacier thinning is however evident: to compare our findings against those of Rota et al. [27], we selected an area below 2660 m a.s.l., corresponding to 0.7 km 2 . The average ice loss was computed as 62.92 ± 0.81 m from 1991 to 2014, equal to a thinning rate of 2.74 ± 0.03 my −1 , and higher than the values reported by Rota et al. [27], i.e.

Conclusion
In this study, we analyzed the evolution of Lys Glacier, one of the largest glaciers in the Italian Alps, by looking at a variety of parameters: terminus fluctuations were studied from historical sources and glaciological bulletins from 1812 to 2017; changes in surface, debris cover and area of supraglacial/proglacial lakes together with volume changes were examined from cartographic and remote sensing datasets from 1975 to 2014. The glacier length variations were found to be similar to those of large glaciers in the Alps such as Mer de Glace and Unterer Grindelwald, indicating a similar climatic setting in spite of the distance of these glaciers; the worst conditions for the glacier development occurred after the end of the Little Ice Age and since 1985 (−443 m from 1985 to 2017) reflecting increasing temperatures as seen from the closest weather station located at Gressoney d'Ejola. Overall, Lys Glacier has retreated by almost 1.6 km since the LIA.
All the other glaciological findings point to a strong glacier reduction, which is interpreted as an evident impact of climate change: the rate of area change was −0.04 km 2 y −1 since 1988, while glacier volume decreased by −47 × 10 6 m 3 from 1991 to 2014. The glacier debris cover increased from 1988 to 2000, when it covered 12.4% of the glacier area, and then started decreasing again, as a result of glacier shrinking, while the area of the proglacial lakes grew exponentially over the same period. We consider the changes in area and debris cover as highly reliable in view of our accuracy assessment (max 2% error in the glacier outlines and accuracy between 90% and 100% when mapping debris cover), while the uncertainty in volume variations is larger because of the lower quality of the input DEM from 1991.
In view of the present conditions of the glacier, which prevent reaching the glacier tongue, remote sensing remains as the only viable option to investigate the glacier variations in the future, while the detachment of the glacier tongue has further complicated studying terminus fluctuations. To further our understanding of the glacier past conditions, other historical sources should be considered, including pictorial documents to lengthen the record of glacier terminus position and aerial photography from the past century to provide more accurate estimates of volume changes.