Michael Zemp

Christophe Delacourt

Simon Gascoin

Delphine Six

Pascal Sirguey

Marie Dumont

Fanny Brun

Marion Réveillet

Maxim Lamare

Jakob Abermann

Katy Medina

Léna Gurriaran

Maxim, Bernard, Thibaud Merci Donc À Gilles

Anne Gag

Kévin, Ugo, Jordi, Maria, Gabbi, Inès, Céline Les Juliens

Lucas Davaze

Antoine Rabatel
email: antoine.rabatel@univ-grenoble-alpes.fr

Yves Arnaud

Romain Hugonnet

Etienne Berthier

Frank Paul

Xiaofeng Li

Prasad S Thenkabail

Article 3 : Region-wide annual glacier surface mass balance for the European Alps from 2000 to 2016

Keywords: Glaciers, Surface mass balance, Optical remote sensing

vi de possibles conditions météorologiques favorisant les fluctuations observées ; (3) calibrer et valider les modèles glacio-hydrologiques utilisés pour estimer les contributions actuelles et futures des glaciers de montagne au fonctionnement hydrologique des bassins versants et à l'élévation du niveau des mers. vii REMERCIEMENTS Au terme de ces trois années de recherche, d'échanges et d'aventures, il est temps de remercier tous ceux qui ont compté.

Je remercie ainsi les membres du jury, pour avoir étudié mes travaux, fait le déplacement pour ma soutenance, mais aussi pour leurs remarques et la discussion qui a suivi la soutenance et que j'ai énormément apprécié. Merci

ABSTRACT

Beyond their iconic role of climate change, mountain glaciers can be considered as Earth' essential component and natural "climate-meter". Despite their small spatial coverage (0.5% of emerged land), mountain glaciers contribute as high as 30% of the observed sea-level rise. In some regions, they are considered as essential issues because of their importance in terms of potable water, agriculture, hydroelectricity or natural hazards. A small share is however monitored in situ (<0.0025%) and their fluctuations at regional scale are poorly known.

Thanks to the development of high spatial resolution satellite sensors (metric to decametric) new methods are today available to significantly increase the number of monitored glaciers. After a state of the art of the existing methods and an identification of the limitations, we focused our attention on the development of two methods.

The first one is based on the automatic detection of the snow/ice interface altitude (i.e. snowline) at the glacier surface from optical satellite images. This altitude, when estimated at the end of summer, is a proxy of the annual glacier-wide mass change at the glacier surface (called surface mass balance, SMB). Using this approach, we estimated the annual SMBs of 239 glaciers in the European Alps and 82 glaciers in the tropical Andes for the period 2000-2016 and 2000-2018, respectively. The mean mass losses are -0.74 and -1.29 m water equivalent per year for the two regions, respectively. This approach allowed deriving the first dataset of annual SMBs for individual glaciers at regional scale from optical remote sensing. We found significant relationships between the computed SMBs and the glacier morphotopographic features (e.g. slope, median altitude …), with steeper and higher glaciers, experiencing less mass losses. Comparison with in situ monitored SMBs revealed an overestimation of mass losses from in situ estimates, due to a low representativeness of steep glaciers (>20°) in the in situ datasets. Our study also revealed heterogeneous inter-annual variability across the European Alps, partially explained by the climatic context of the studied sub-regions, thanks to the analysis of climate reanalysis data.

We developed a second method to derive the annual and summer SMBs from albedo maps, computed from MODIS images. With an application on 30 glaciers in the French Alps, this work opened the way toward a regional application of this method, in order to estimate both annual and summer SMBs. By performing regional applications on different glacierized regions, we developed and validated methods capable of deriving the annual and summer SMBs of individual mountain glaciers at regional scale, from optical remote sensing data. These data could then be used to [START_REF][END_REF] assess the impact of peculiar climatic conditions onto mountain glaciers; (2) investigate possible meteorological conditions driving the documented glacier fluctuations; (3) calibrate and validate glacio-hydrological models used to estimate the current and future contributions of mountain glaciers to the hydrological functioning of mountain catchments and to sea level rise. v RESUME Au-delà de leur rôle d'icône du changement climatique, les glaciers de montagne sont une composante essentielle de notre planète. Ils sont, de plus, de véritables « climat-mètres » naturels. Malgré leur faible superficie (0.5% des terres émergées), les glaciers de montagne contribuent à hauteur de 30% à la hausse du niveau des mers. Dans certaines régions, ils constituent de véritables enjeux quant à l'eau potable, l'agriculture, la production hydroélectrique ou les aléas glaciaires. Peu sont en revanche instrumentés (<0.0025%) et leurs fluctuations à l'échelle de régions entières sont mal connues.

Grâce au développement de capteurs satellitaires à haute résolution spatiale (métrique à décamétrique), le développement de méthodes automatisées permet aujourd'hui d'augmenter considérablement le nombre de glaciers observés. Après avoir dressé un état de l'art des méthodes existantes et identifié les verrous méthodologiques, nous avons développé deux méthodes en particulier.

La première se base sur la détection automatique de l'altitude de la limite glace/neige (i.e. ligne de neige) à la surface du glacier, à partir d'images satellites optiques. Cette altitude, lorsque mesurée à la fin de l'été, est un marqueur du changement de masse à la surface du glacier ayant eu lieu au cours de l'année (appelé bilan de masse de surface). Cette approche nous a permis d'estimer le bilan de masse de surface annuel de 239 glaciers dans les Alpes européennes et de 82 glaciers dans les Andes tropicales pour la période 2000-2016 et 2000-2018, respectivement. La perte moyenne annuelle observée est de -0.74 et de -1.29 m équivalent eau par an pour les deux régions respectivement. A notre connaissance, cette approche a permis d'établir le premier jeu de données de bilans de masse de surface annuels pour des glaciers individuels à échelle régionale à partir d'images satellites optiques. Une dépendance du bilan de masse de surface moyen par glacier à des caractères morpho-topographiques (e.g. pente, altitude médiane …) a été observée, où plus les glaciers sont pentus et hauts en altitude, moins leur perte de masse est importante. Une comparaison avec des mesures in situ dans les Alpes Européennes révèle une surestimation de la perte de masse par ces dernières si on les extrapole spatialement, notamment à cause de la faible représentation de glaciers à forte pente (>20°) dans les mesures in situ. Notre étude sur les Alpes Européennes a de plus permis d'identifier une variabilité interannuelle hétérogène sur cette région, en partie expliquée par des contextes climatiques différents grâce à l'utilisation de données issues de ré-analyses.

Le développement d'une autre méthode a permis, à partir de l'analyse de cartes d'albédo issues du capteur MODIS, de caractériser le bilan de masse de surface annuel et estival de 30 glaciers dans les Alpes françaises. Cette étude ouvre la porte à l'utilisation de cette méthode pour l'analyse du bilan annuel et saisonnier à l'échelle régionale.

Ce travail a permis, à travers des applications dans différentes régions englacées, de développer et valider des méthodes capables, à partir d'images satellites optiques, d'estimer le bilan de masse de surface annuel et saisonnier de glaciers de montagne à l'échelle de régions entières. Ces estimations peuvent ensuite être utilisées pour : [START_REF][END_REF] Au-delà de leurs rôles esthétiques et d'icônes du changement climatique [START_REF] Mackintosh | Reconstructing Climate from Glaciers[END_REF], les glaciers de montagne sont des composants essentiels de notre planète. Intimement liés au climat, ils répondent rapidement aux variations climatiques et sont de véritables "climat-mètres" [START_REF] Vaughan | Observations: Cryosphere[END_REF], faisant d'eux une des variables essentielles d'observation du climat par l'organisation météorologique mondiale (OMM, [START_REF] Bojinski | The Concept of Essential Climate Variables in Support of Climate Research, Applications, and Policy[END_REF]. Le récent dérèglement climatique, en partie dû à l'impact anthropique, a renforcé la nécessité de les observer et de les étudier, notamment dans des régions où les données climatiques directes (i.e. température, précipitations …) sont rares et où les glaciers de montagne constituent un des rares marqueurs.

Les glaciers ont, de plus, une importance hydrologique capitale à l'échelle régionale (e.g. ressource en eau, production d'électricité) et globale (e.g. hausse du niveau des mers). Malgré leur faible couverture spatiale à l'échelle des terres émergées (0.5%, [START_REF] Vaughan | Observations: Cryosphere[END_REF], les glaciers de montagne sont présents dans 26% des grands bassins versants mondiaux, hébergeant près d'un tiers de la population mondiale [START_REF] Beniston | Climatic Change in Mountain Regions: A Review of Possible Impacts[END_REF]. Dans certaines régions des Hautes Montagnes d'Asie (HMA) ou des Andes tropicales, ils constituent des réservoirs d'eau indispensables à l'approvisionnement des populations en eau potable et pour l'irrigation des cultures. Ces régions sont ainsi particulièrement impactées par l'augmentation de la perte de masse de ces glaciers qui ont, pour certains, déjà dépassés leurs contributions maximales aux débits hydrologiques observés (i.e. peak water) lors de la saison de fonte (e.g [START_REF] Huss | Global-scale hydrological response to future glacier mass loss[END_REF] et verront la contribution glaciaire aux débits saisonniers réduite (e.g., [START_REF] Barnett | Potential impacts of a warming climate on water availability in snow-dominated regions[END_REF]Kaser et al., 2010).

A l'échelle globale, une partie de ces eaux de fonte convergent vers les océans et participent ainsi à la hausse du niveau des mers si elles ne sont pas compensées par l'évaporation. En dépit de leur faible étendue comparée aux grandes calottes glaciaires (0.5% versus 11% de la surface totale des terres émergées), les glaciers contribuaient à hauteur de 0.76 ± 0.37 mm par an sur Introduction 3 Historiquement, les glaciers de montagne sont mesurés in situ en estimant leurs variations de masse au cours du temps : le bilan de masse de surface glaciaire (SMB pour l'acronyme en anglais). Cette mesure permet d'estimer avec précision l'impact du climat sur le glacier. Cependant, il existe des données de SMB pour peu de glaciers (environ 0.002%, WGMS, 2017) du fait de l'investissement humain et logistique requis. Dans le but d'affiner notre compréhension des causes, impacts et enjeux du dérèglement climatique sur l'évolution des glaciers à l'échelle régionale, que cela soit pour la ressource en eau, la hausse du niveau des mers, ou une meilleure adaptabilité envers les aléas glaciaires, le nombre d'observations doit augmenter, à la fois spatialement et temporellement. La modélisation est un outil intéressant qui permet de s'affranchir des contraintes spatio-temporelles et des modèles récents permettent de simuler l'évolution passé et future de la totalité des glaciers terrestres (e.g. OGGM, Maussion et al., 2019, GloGEMflow, Zekollari et al., 2019). Ces modèles sont, en revanche, fortement dépendants des mesures de SMB in situ, souvent biaisées négativement à cause de leur faible représentativité spatiale (Gardner et al., 2013;[START_REF] Marzeion | Observation-Based Estimates of Global Glacier Mass Change and Its Contribution to Sea-Level Change[END_REF]. Ce biais peut être atténué par l'apport de méthodes utilisant des images satellites, telle que la géodésie [START_REF] Zemp | Reanalysing glacier mass balance measurement series[END_REF]. Aujourd'hui, il n'existe en revanche pas de méthode fiable permettant, à partir d'images satellites, d'estimer les variations de SMB de plusieurs centaines de glaciers de montagne (échelle régionale) au pas de temps annuel et/ou saisonnier. C'est donc dans ce contexte que mon travail de thèse s'inscrit, afin de développer et de valider des méthodes capables de produire des jeux de données de SMB par glacier à l'échelle régionale et aux pas de temps annuels et saisonniers, en profitant de la récente grande disponibilité d'images satellites optiques à haute résolution.

ORGANISATION DU DOCUMENT

Ce manuscrit suit un fil général allant de la présentation des méthodes aujourd'hui disponibles pour estimer le SMB des glaciers de montagne, à la présentation de deux méthodes développées ainsi que l'analyse des jeux de données associés. Le manuscrit débute par un état de l'art épistémologique (Sect. 2.1) narrant les premières observations glaciologiques et les grandes questions qui ont fait débats au sein de la communauté glaciologique. Une présentation des principales méthodes permettant l'estimation du SMB est ensuite présentée (Sect. 2.2 et 2.3), dans le but de mettre en exergue les capacités et limites de chaque méthode (Sect. 2.4). Cette synthèse permet de rappeler mes objectifs de recherche, et l'angle d'attaque choisi pour lever les verrous scientifiques, méthodologiques et techniques identifiés. Un résumé de mes recherches est ensuite présenté Sect. 3. Il contient un résumé des différents articles que j'ai pu écrire ou co-écrire et vise à replacer le lecteur dans la démarche méthodologique que j'ai suivie pendant ces trois années.

Sect. 4, 5, 6 et scientifique pour les glaciers de montagne prend de l'ampleur avec, notamment, Jean-Jacques Scheuchzer et Johann Heinrich Hottinger [START_REF] Rémy | Mais comment s'écoule donc un glacier ? Aperçu historique[END_REF][START_REF] Zryd | Les glaciers en mouvement: la population des Alpes face au changement climatique, Collection Bibliographie 209 le savoir suisse[END_REF]. Contrairement aux préoccupations scientifiques actuelles telles que l'estimation de la contribution des glaciers à la hausse du niveau des mers, ou l'estimation d'aléas d'origines glaciaires, le Un obstacle épistémologique décrit un concept, une loi ou un théorème très complexe à comprendre en mobilisant uniquement ses connaissances ou observations. Franchir un obstacle épistémologique revient à le comprendre dans son intégralité et donc à intégrer son sens profond.

Box 1 : Obstacle épistémologique premier obstacle épistémologique (Box 1) auquel les scientifiques se sont heurtés, et qui occupera toute l'attention de ces derniers jusqu'à la fin du XIXème siècle, est le mécanisme à l'origine de l'écoulement du glacier. La première théorie, décrite par Scheuchzer en 1705 et qui sera, plus de 100 ans plus tard soutenue par Louis Agassiz, est celle de la dilatation [START_REF] Agassiz | Remarks on glaciers[END_REF].

Scheuchzer explique que l'écoulement glaciaire est dû aux cycles successifs de gel-dégel d'eau liquide dans les anfractuosités de la glace, déclenchant le mouvement observé d'amont en aval. Cette théorie expliquait notamment les variations de vitesse observées à la surface du glacier entre l'été (plus rapide) et l'hiver (plus lent), confirmées de nos jours par des mesures in situ (e.g. [START_REF] Vincent | Sliding velocity fluctuations and subglacial hydrology over the last two decades on Argentière glacier, Mont Blanc area[END_REF]. Altmann, puis Gruner et De Saussure en 1760 élaborent parallèlement la théorie de la gravitation [START_REF] Seligman | Research on Glacier Flow[END_REF]. Cette théorie concilie la rigidité apparente de la glace avec le mouvement observé en proposant un écoulement uniquement dû au glissement d'un seul tenant du glacier sur son socle. Enfin, en 1773, André César Bordier est le premier à [START_REF] Walker | Early Discoverers XXXV: Descent of Glaciers: Some Early Speculations on Glacier Flow and Ice Physics[END_REF]. "Abri construit par Hugi en 1842 sur la moraine médiane de l'Unteraar, sous la confluence du Lauteraar (à gauche) et du Finsteraar (à droite). Cet abri, construit à la base de l'éperon rocheux au centre de l'image a atteint l'endroit représenté en 13 années. Cette moraine a été l'endroit où Agassiz installa l'Hôtel des Neuchâtelois pendant ses recherches en 1840." Méthodes 7 reconnaître une certaine viscosité de la glace, la comparant à de la "cire amollie flexible" et introduit la notion d'écoulement du glacier d'amont en aval (Bordier, 1773). La théorie de Bordier sera reprise en 1840 par le chanoine Rendu (Rendu, 1840) (Hugi, 1842, Fig. 1).

En 1840, Louis Agassiz retourne sur l'Unteraar pour poursuivre les mesures d'écoulement, en construisant une station de mesure sur une moraine médiane, "l'Hôtel des Neuchâtelois". Il fore pour cela plusieurs trous n'excédant pas sept mètres et y insère des jalons pour mesurer le déplacement de la glace et s'affranchir de l'impact des débris supraglaciaires contrairement à Hugi [START_REF] Clarke | A short history of scientific investigations on glaciers[END_REF]. En retournant sur le glacier l'année suivante (1841) et notamment accompagné du physicien écossais James David Forbes, il s'aperçoit que les jalons ne sont, pour la plupart, plus enchâssés dans la glace à cause de la fonte en surface du glacier. Sans le savoir, il pose les bases de la méthodologie encore utilisée aujourd'hui pour mesurer in situ le bilan de masse de surface glaciaire (SMB, Le bilan de masse glaciaire correspond à la somme de tous les flux de masse entrant et sortant du glacier. Le taux de changement de masse peut être calculé sur la surface du glacier comme la somme des flux de masse de surface, interne et basal (Cuffey and Paterson, 2010) sur un intervalle donné. Ce taux est exprimé en équivalent en eau grâce à une hypothèse sur la densité de la glace/neige. Le terme dominant étant, pour la plupart des glaciers, le bilan de masse de surface, le bilan de masse glaciaire est souvent approximé comme le bilan de masse de surface glaciaire (SMB).

Box 2 : Bilan de masse de surface (SMB) Box 2). C'est finalement Forbes qui, en exportant les méthodes d'Agassiz sur la Mer de Glace en 1842, donne les premières clefs du mystère de l'écoulement des glaciers, en proposant une déformation de la glace sous son propre poids, entraînant l'écoulement de celui-ci. Il se base pour cela sur des mesures de vitesse de surface effectuées au théodolite sur la Mer de Glace, validant ce que André César Bordier et le chanoine Louis Rendu avaient déjà compris par le raisonnement dès 1773 et 1840. Malgré son apparente rigidité, la glace a un comportement viscoplastique, lui permettant de se déformer, mais aussi de se fracturer ! De plus, Forbes complète ses observations par des expériences et admet que l'écoulement peut être facilité par la présence d'eau liquide à la base du glacier, agissant comme lubrifiant [START_REF] Forbes | Travels Through the Alps of Savoy and Other Parts of the Pennine Chain: With Observations on the Phenomena of Glaciers[END_REF]. Peu consensuelle à son origine, cette théorie sera néanmoins confortée par Weinberg en 1907 qui, appliquant les lois de Navier-Stokes au glacier, validera la nature visqueuse de la glace. extremity at any time seems to retire, this implies nothing more than that the portion of ice, melted by the heat of summer, is more considerable than that which the glacier brings along with it in its progress" (Agassiz, 1839, p. 383).

Ce postulat est d'une importance capitale car il exprime simplement le fonctionnement d'un glacier, ainsi que les causes qui peuvent l'amener à subir des changements de longueur et/ou de volume. De plus, il met en exergue le lien étroit entre le SMB, et le climat, faisant des glaciers un observatoire du climat, et dont l'estimation du SMB est la mesure la plus pertinente.

Les premières mesures se sont majoritairement concentrées sur les variations de position du front du glacier, dans plusieurs régions simultanément (Europe, Scandinavie et Amérique du Nord). Ces mesures ont été étoffées par l'engouement des peintres de l'époque pour ces géants de glace, dont Samuel Birmann a été l'un des plus actifs, et qui permis, bien plus tard, la reconstruction des positions du front de plusieurs glacier (e.g.

Méthodes 9

Nussbaumer [START_REF] Nussbaumer | The Little Ice Age history of the Glacier des Bossons (Mont Blanc massif, France): a new high-resolution glacier length curve based on historical documents[END_REF]Zumbühl et al., 2008;Zumbühl and Nussbaumer, 2018) à une époque où peu de mesures existaient, comme sur le versant Nord du massif du Mont-Blanc par exemple (Fig. 2). Malgré leur lien avec les changements de volume glaciaire, les variations de longueur n'en sont qu'une mesure indirecte. Ces variations de longueurs sont fortement dépendantes de la géométrie du glacier, de ses caractéristiques morphotopographiques ainsi que de sa dynamique d'écoulement (Cuffey and Paterson, 2010), qui conditionnent son temps de réponse au forçage climatique, propre à chaque glacier.

C'est véritablement François-Alphonse Forel qui, dès les années 1890, met en place des campagnes de mesures glaciologiques régulières (Rapports sur les variations périodiques des glaciers, Forel, 1895) dans le cadre de la création de la Commission Internationale des Glaciers, fondée en 1894 [START_REF] Hall | A Glacial Committee for North America and Europe[END_REF] ancêtre du World Glacier Monitoring Service (WGMS), collectant les données de la plupart des glaciers instrumentés. Les mesures effectuées dès cette période ne se cantonnent pas uniquement à la position du front. Dans son discours préliminaire, Forel introduit la notion de variation de grandeur, pour quantifier les variations de longueur, de largeur, d'épaisseur et donc de volume. Un protocole expérimental intitulé "Méthode du glacier du Rhône" est ainsi élaboré. Développé par des ingénieurs suisses, il est basé sur des mesures périodiques de retraits de la langue glaciaire, de variations d'épaisseur le long de profils transversaux ainsi que de vitesses d'écoulement en surface. Cela mènera aux premières séries de SMB. Il faudra attendre la fin de la seconde guerre mondiale pour voir cette méthode se perfectionner puis être appliquée à un plus grand nombre de glaciers [START_REF] Haeberli | CLIMATE CHANGE AND ITS IMPACTS ON GLACIERS AND PERMAFROST IN THE ALPS[END_REF].

MESURER LE BILAN DE MASSE DE SURFACE A PARTIR DE DONNEES IN SITU

La méthode glaciologique

Au cours du XXème siècle, des méthodes ont été développées afin d'estimer in situ le SMB. Une de ces méthodes, la méthode glaciologique consiste à mesurer in situ, l'accumulation et l'ablation à des sites choisis (appelés balises) en fonction de la morpho-topographie de la surface du glacier, puis d'extrapoler ces valeurs à l'ensemble du glacier. Ces balises sont implantées sur la zone d'ablation et d'accumulation, comme illustré Fig. 3 sur plusieurs glaciers français. On peut noter un nombre de balises généralement plus réduit dans la zone d'accumulation par rapport à la zone d'ablation. Cette différence peut provenir de la difficulté d'accès aux balises ainsi que l'effort plus important pour mesurer l'accumulation. Ce sous-échantillonnage dans la zone d'accumulation peut conduire à une mauvaise caractérisation de l'accumulation, notamment à cause de la grande variabilité spatiale de l'enneigement, due à la topographie locale parfois complexe ainsi qu'au processus éoliens. Les relevés de balises sont effectués le plus près possible de la fin de la saison d'ablation et d'accumulation (respectivement vers le 1er octobre et 1er mai dans les Alpes). Le protocole de mesure diffère en fonction du type de balise. Pour les balises disposées en zone d'ablation, des jalons sont introduits dans la glace (autour du 1er octobre dans les Alpes, en fonction des conditions météorologiques), à l'aide d'une foreuse ou d'une sonde à vapeur. Un an d'après, la partie émergée des jalons est mesurée, permettant ainsi de quantifier la hauteur de glace ayant fondue pendant l'année hydrologique (schéma illustratif Fig. 4a). Dans la zone d'accumulation, le protocole décrit ne peut être utilisé, car le jalon serait perdu d'une année à Méthodes 11 l'autre. Vers le 1er mai, l'épaisseur de neige accumulée pendant la période hivernale est estimée, soit par carottage, soit grâce à des puits de neige (i.e. "snowpit"), comme illustrée Fig. 4b. Afin de déterminer une hauteur indépendante de la densité de la neige, l'épaisseur mesurée est convertie en équivalent en eau (w.e.) grâce à des mesures de densité. La différentiation entre la neige accumulée durant l'année et le névé est effectuée visuellement grâce à l'identification d'un horizon plus sombre ainsi qu'à une variation importante de densité et de la taille des grains. L'hypothèse est ensuite faite que chaque balise (dénommée i) est représentative des processus de surface au sein d'une zone de surface Si (illustration Fig. 5). Pour connaître le SMB à l'échelle du glacier (appelé glacier-wide) à partir de mesures ponctuelles aux balises, l'équation 1 est utilisée : Cette méthode, recommandée par le WGMS et largement utilisée dans toutes les régions englacées, est notamment détaillée dans les articles de Østrem and Brugman (1991), et [START_REF] Kaser | A manual for monitoring the mass balance of mountain glaciers[END_REF]. Ainsi en 2015, 166 glaciers sont instrumentés pour déterminer leurs bilans de masse de surface, en utilisant majoritairement cette méthode, constituant ainsi un jeu de données de référence dont 40 glaciers ont aujourd'hui plus de 40 ans d'observation [START_REF] Wgms | global glacier change bulletin -world glacier monitoring service[END_REF]. La plupart des incertitudes attachées à cette méthode proviennent de :



La distribution et la représentativité des balises à la surface du glacier et la méthode d'extrapolation utilisée [START_REF] Escher-Vetter | Four decades of winter mass balance of Vernagtferner and Hintereisferner, Austria: methodology and results[END_REF][START_REF] Hock | Application of Kriging Interpolation for Glacier Mass Balance Computations[END_REF][START_REF] Thibert | Glaciological and volumetric mass-balance measurements: error analysis over 51 years for Glacier de Sarennes, French Alps[END_REF]  Les changements de superficie et d'altitude de la surface du glacier pendant la période d'observation [START_REF] Zemp | Reanalysing glacier mass balance measurement series[END_REF]. Ces erreurs peuvent être réduites en utilisant des modèles numériques de terrain (MNT) pour contraindre le SMB à la date des acquisitions photogrammétries. Cette méthode, appelée méthode géodésique, peut être utilisée indépendamment des mesures in situ. Nous y reviendrons par la suite.

 L'erreur de mesure à chaque balise, provenant de la mesure de la densité ainsi que de la rugosité de surface [START_REF] Thibert | Glaciological and volumetric mass-balance measurements: error analysis over 51 years for Glacier de Sarennes, French Alps[END_REF]. Une approche statistique alternative a été développée par Louis [START_REF] Lliboutry | Multivariate Statistical Analysis of Glacier Annual Balances[END_REF] pour quantifier le bilan de l'ensemble du glacier à partir des mesures aux balises. Elle peut également être utilisée pour reconstruire des données manquantes aux balises, notamment en lien avec les aléas de mesure (balises manquantes ou détruites). Cette approche se base sur les mesures acquises aux balises, permettant de capturer la variabilité spatiale et temporelle des bilans de masse aux balises. L'hypothèse est ensuite faite que la variabilité spatiale observée sera reproduite chaque année, permettant en cas d'absence de mesures, de reconstruire celle-ci (voir la description de la méthode appliquée aux bilans de masse, [START_REF] Reynaud | Mass-Balance Measurements: Problems and two New Methods of Determining Variations[END_REF] [START_REF] Basantes-Serrano | Slight mass loss revealed by reanalyzing glacier mass-balance observations on Glaciar Antisana 15<span class="italic">α</span> (inner tropics) during the 1995-2012 period[END_REF]. Elle est aussi adaptée à d'autres types de données, pour compléter des données manquantes dans des séries de l'altitude de la ligne d'équilibre glaciaire par exemple (Rabatel et al., 2013a[START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF].

La méthode hydrologique

Une autre méthode pour mesurer le SMB in-situ consiste à mesurer les flux d'eau entrant et sortant du glacier, indépendamment de mesures directes sur le glacier (Cuffey and Paterson, 2010;[START_REF] Tangborn | A comparison of glacier mass balance by glaciological, hydrological and mapping methods, South Cascade Glacier, Washington[END_REF] [START_REF] Tangborn | A comparison of glacier mass balance by glaciological, hydrological and mapping methods, South Cascade Glacier, Washington[END_REF] non pris en comptes par la méthode glaciologique, la perte par infiltration [START_REF] Viani | Glacier-wide summer surface mass-balance calculation: hydrological balance applied to the Argentière and Mer de Glace drainage basins (Mont Blanc)[END_REF] ou liée à des défaillances instrumentales, notamment pour estimer la quantité de précipitations [START_REF] Sicart | Glacier mass balance of tropical Zongo glacier, Bolivia, comparing hydrological and glaciological methods[END_REF]. [START_REF] Braithwaite | Calculation of Glacier Ablation from Air Temperature, West Greenland[END_REF], soit en considérant aussi, par exemple, la radiation solaire pour mieux reproduire la variabilité spatiale de la fonte [START_REF] Hock | A distributed temperature-index ice-and snowmelt model including potential direct solar radiation[END_REF][START_REF] Pellicciotti | An enhanced temperatureindex glacier melt model including the shortwave radiation balance: development and testing for Haut Glacier d'Arolla, Switzerland[END_REF].

Les approches empiriques sont en revanche moins nombreuses pour quantifier le bilan hivernal notamment à cause de la forte variabilité de l'accumulation en fonction de la topographie de surface du glacier et de son contexte orographique. Des approches très simples se basant sur la corrélation entre les précipitations en vallée et l'accumulation de neige à la surface de glaciers sont assez largement utilisées; elles se basent sur des données issues de stations météorologiques, de réanalyses, où des modèles climatiques régionaux [START_REF] Hock | A distributed temperature-index ice-and snowmelt model including potential direct solar radiation[END_REF][START_REF] Machguth | Distributed glacier mass-balance modelling as an important component of modern multi-level glacier monitoring[END_REF][START_REF] Vincent | Snow accumulation and ice flow at Dôme du Goûter (4300 m)[END_REF]. Des approches plus complexes ont récemment été développés, prenant en compte le contexte topographique de la surface [START_REF] Huss | Determination of the seasonal mass balance of four Alpine glaciers since 1865[END_REF][START_REF] Sold | Methodological approaches to infer end-of-winter snow distribution on alpine glaciers[END_REF][START_REF] Winstral | Simulating wind fields and snow redistribution using terrain-based parameters to model snow accumulation and melt over a semi-arid mountain catchment[END_REF]) mais de grandes incertitudes demeurent sur les bilans hivernaux estimés (Réveillet et al., 2018), pouvant avoir un impact non-négligeable sur le bilan annuel modélisé [START_REF] Réveillet | Which empirical model is best suited to simulate glacier mass balances?[END_REF].

D'autres approches assimilent des données de télédétection comme l'altitude de la ligne de neige afin d'améliorer la paramétrisation des modèles degrés-jour, en prenant en compte l'état de la surface du glacier (e.g. Barandun et al., 2018;[START_REF] Tawde | Estimation of Glacier Mass Balance on a Basin Scale:An Approach Based on Satellite-Derived Snowlines and a Temperature Index Model[END_REF]. Certaines de ces approches permettent d'établir, des projections du devenir des glaciers alpins pour le XXIème siècle (e.g. [START_REF] Zemp | Alpine glaciers to disappear within decades?[END_REF], la contribution passée et future des glaciers au niveau des mers (e.g. Marzeion et al., 2012;[START_REF] Radić | Regional and global projections of twenty-first century glacier mass changes in response to climate scenarios from global climate models[END_REF], ou à l'échelle des grands bassins versants englacés (Bolibar et al., en discussion).

2. Les modèles basés sur la quantification des processus de surface.

Certains de ces modèles se basent notamment sur la résolution du bilan d'énergie de surface (BES -Box 3) pour estimer le SMB. Ce bilan est aussi utilisé pour estimer la fonte à partir de données de stations météorologiques in-situ qui permettent le calcul direct (méthode d'eddy covariance, [START_REF] Swinbank | The measurement of vertical transfer of heat and water vapor by eddies in the lower atmosphere[END_REF] ou indirect [START_REF] Fitzpatrick | Surface Energy Balance Closure and Turbulent Flux Parameterization on a Mid-Latitude Mountain Glacier, Purcell Mountains, Canada[END_REF][START_REF] Gravesen | A new programme for monitoring the mass loss of the Greenland ice sheet[END_REF][START_REF] Thornthwaite | The determination of evaporation from land and water surfaces[END_REF][START_REF] Van As | Warming, glacier melt and surface energy budget from weather station observations in the Melville Bay region of northwest Greenland[END_REF][START_REF] Wagnon | Energy balance and runoff seasonality of a Bolivian glacier[END_REF] de toutes les composantes du BES.

TELEDETECTION DU BILAN DE MASSE

"Si un observateur pouvait être transporté à une assez grande hauteur au dessus des Alpes, pour embrasser d'un coup d' oeil celles de la Suisse, de la Savoie et du Dauphiné, il verrait cette chaîne de montagne sillonnée par de nombreuses vallées et composée de plusieurs chaînes parallèles, la plus haute au milieu et les autres décroissant graduellement, à mesure qu'elles s'en éloignent."

Horace Bénédicte de Saussure (1803) invoquant le soutien de la télédétection pour l'étude des glaciers de montagne.

Les avancées technologiques dans le domaine de l'imagerie spatiale et aéroportée ont permis de quantifier les changements de longueur et de volume glaciaires dès le début des années 1970. Les premières observations se sont principalement concentrées sur l'observation de caractéristiques de surface telles que l'altitude de la ligne d'équilibre glaciaire, la taille de la zone d'accumulation, aussi dénommée AAR pour accumulation area ratio (Box 4 and 5;[START_REF] Chinn | Glacier snow line variations in the Southern Alps, New Zealand[END_REF]Dyurgerov et al., 2009;[START_REF] Kaser | Changes of the equilibrium-line altitude in the tropical Cordillera Blanca, Peru, 1930-50, and their spatial variations[END_REF]LaChapelle, 1962;Meier and Post, 1962;[START_REF] Østrem | Erts Data in Glaciology-An Effort to Monitor Glacier Mass Balance from Satellite Imagery[END_REF], ou encore les variations de longueur à partir de photographies aériennes (e.g. [START_REF] Hubley | Glaciers of the Washington Cascade and Olympic Mountains; Their Present Activity and its Relation to Local Climatic Trends[END_REF]LaChapelle, 1962). Certaines études ont, par exemple, permis de mettre en évidence certaines crues glaciaires (e.g., [START_REF] Bindschadler | Geometry and Dynamics of a Surge-type Glacier[END_REF][START_REF] Post | The Recent Surge of Walsh Glacier, Yukon and Alaska[END_REF] et aléas glaciaires [START_REF] Kääb | REMOTE SENSING OF GLACIER-AND PERMAFROSTRELATED HAZARDS IN HIGH MOUNTAINS: AN OVERVIEW[END_REF]Kääb et al., 2018). L'apparition de satellites civils avec notamment le lancement des missions LANDSAT dès 1972 permit d'augmenter drastiquement le nombre d'observations satellitaires et d'être graduellement capable d'estimer le bilan de masse par télédétection (Braithwaite, 1984;[START_REF] Chinn | Glacier snow line variations in the Southern Alps, New Zealand[END_REF][START_REF] Ruyter De Wildt | A method for monitoring glacier mass balance using satellite albedo measurements: application to Vatnajökull, Iceland[END_REF][START_REF] Greuell | Assessment of the surface mass balance along the K-transect (Greenland ice sheet) from satellite-derived albedos[END_REF][START_REF] Østrem | Erts Data in Glaciology-An Effort to Monitor Glacier Mass Balance from Satellite Imagery[END_REF]Rabatel et al., 2005[START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF]Rabatel et al., , 2017)). De plus, l'apparition de satellites capables d'acquisitions La ligne d'équilibre glaciaire (ELA pour equilibrium-line altitude) est l'altitude à laquelle l'accumulation est exactement compensée par l'ablation et donc où le SMB au cours d'une année est nul.

Box 4 : ELA stéréoscopiques a ouvert la porte à l'application de la méthode géodésique à partir de données spatiales, permettant de mesurer directement les variations de volume glaciaire (e.g. Berthier et al., 2016;[START_REF] Brun | A spatially resolved estimate of High Mountain Asia glacier mass balances, 2000-2016[END_REF]Gardelle et al., 2013;[START_REF] Menounos | Heterogeneous Changes in Western North American Glaciers Linked to Decadal Variability in Zonal Wind Strength[END_REF][START_REF] Ragettli | Heterogeneous glacier thinning patterns over the last 40 years in Langtang Himal, Nepal[END_REF]Shean et al., 2016). Plusieurs méthodes sont aujourd'hui "opérationnelles" pour estimer le SMB des glaciers de montagne par télédétection et sont présentées succinctement dans les paragraphes suivants.

Ces méthodes de télédétection peuvent être divisées en deux catégories : la première regroupe les méthodes mesurant directement une variation de masse et/ou de volume, alors que la seconde regroupe les méthodes utilisant des caractéristiques de surface comme indicateur du SMB.

Les méthodes directes

La méthode gravimétrique

La seule méthode en théorie capable, à partir de données satellites, de mesurer le bilan de masse de glaciers et plus généralement un bilan hydrologique est la méthode dite gravimétrique. La gravimétrie consiste à mesurer le champ de pesanteur ainsi que ses variations. En effet, la pesanteur est fortement dépendante de la masse de l'objet considéré. Ainsi, si une variation de masse a lieu à un endroit donné du globe (vidange d'un lac, changement de quantité d'eau dans un aquifère, fonte d'un glacier), ceci se traduit par un changement de pesanteur associé. En mars 2002 ont été lancés deux satellites, GRACE-1 et GRACE-2 (Gravity Recovery And Climate Experiment), dans le but d'observer les anomalies de champ de gravité. En d'autres termes, ils permettent de mesurer les variations de masse en surface et subsurface de la Terre. Pour un usage glaciologique, ces données ont tout d'abord été utilisées pour estimer le bilan de masse du Groenland et de l'Antarctique (e.g. [START_REF] Velicogna | Increasing rates of ice mass loss from the Greenland and Antarctic ice sheets revealed by GRACE[END_REF][START_REF] Velicogna | Measurements of Time-Variable Gravity Show Mass Loss in Antarctica[END_REF][START_REF] Wouters | GRACE observes small-scale mass loss in Greenland[END_REF].

Le "accumulation area ratio" est un ratio permettant de quantifier la part d'un glacier étant en zone d'accumulation par rapport à sa superficie totale. Par exemple, un ratio nul pour un glacier donné et une année donnée correspond à une situation où le glacier n'a plus de zone d'accumulation. De récentes études ont montré que des anomalies de champ de gravité pouvaient être utilisées pour quantifier le bilan de masse de régions englacées non-polaires (e.g. Gardner et al., 2013;Jacob et al., 2012;[START_REF] Luthcke | Recent glacier mass changes in the Gulf of Alaska region from GRACE mascon solutions[END_REF][START_REF] Matsuo | Time-variable ice loss in Asian high mountains from satellite gravimetry[END_REF][START_REF] Yi | Evaluation of glacier changes in high-mountain Asia based on 10 year GRACE RL05 models[END_REF]. En revanche, du fait de la faible résolution spatiale, des incertitudes importantes peuvent être associées aux résultats et ne rendent pas possible l'utilisation de ces données sur des glaciers individuels de moins de 100 km2, comme le montre la carte représentant les variations de masse à l'échelle de toutes les hautes montagnes d'Asie, Fig. 6 (issue de Jacob et al., 2012). De plus, il est parfois délicat de différencier la contribution glaciaire des autres contributions hydrologiques telles que le stockage d'eau souterrain, l'impact de la mousson, le rebond isostatique, le soulèvement tectonique ou érosion de chaînes de montagne [START_REF] Chanard | Modeling deformation induced by seasonal variations of continental water in the Himalaya region: Sensitivity to Earth elastic structure[END_REF][START_REF] Yi | Evaluation of glacier changes in high-mountain Asia based on 10 year GRACE RL05 models[END_REF]. D'autres méthodes sont ainsi plus adaptées à l'estimation du bilan de masse glaciaire par télédétection, en particulier pour les glaciers de montagne [START_REF] Bamber | The land ice contribution to sea level during the satellite era[END_REF]. pour estimer les variations d'épaisseur [START_REF] Schutz | Overview of the ICESat Mission[END_REF], notamment en Antarctique, au Groenland (e.g. [START_REF] Howat | Rates of southeast Greenland ice volume loss from combined ICESat and ASTER observations[END_REF], ainsi que sur d'autres champs de glace et glaciers de tailles suffisantes (e.g. [START_REF] Bolch | Mass loss of Greenland's glaciers and ice caps 2003-2008 revealed from ICESat laser altimetry data[END_REF]Gardner et al., 2013) et al., 2006;[START_REF] Favey | Evaluating the Potential of an Airborne Laser-scanning System for Measuring Volume Changes of Glaciers[END_REF][START_REF] Jóhannesson | Ice-volume changes, bias estimation of mass-balance measurements and changes in subglacial lakes derived by lidar mapping of the surface of Icelandic glaciers[END_REF][START_REF] Rees | Mass balance and dynamics of a valley glacier measured by highresolution LiDAR[END_REF] mais sont souvent onéreuses et logistiquement exigeantes, et, à ma connaissance, seuls certains glaciers d'Alaska bénéficient aujourd'hui d'un suivi annuel systématique aéroporté, même si celui-ci se cantonne parfois à un transect central [START_REF] Arendt | Rapid wastage of Alaska glaciers and their contribution to rising sea level[END_REF][START_REF] Echelmeyer | Airborne surface profiling of Alaskan glaciers[END_REF]. Le LiDAR est davantage utilisé pour la prévention et le suivi d'aléas en lien avec les glaciers (e.g. [START_REF] Vincent | Volume and frequency of ice avalanches from Taconnaz hanging glacier[END_REF], le permafrost (e.g. Guerin et al., 2017;Rabatel et al., 2008), ou encore la quantification de la variabilité spatiale de l'enneigement saisonnier (e.g. [START_REF] Painter | The Airborne Snow Observatory: Fusion of scanning lidar, imaging spectrometer, and physically-based modeling for mapping snow water equivalent and snow albedo[END_REF][START_REF] Revuelto | Mapping the annual evolution of snow depth in a small catchment in the Pyrenees using the longrange terrestrial laser scanning[END_REF].

La méthode géodésique

Cette méthode consiste à mesurer les variations totales du volume du glacier par mesures répétées de sa topographie, non ponctuellement comme l'altimétrie laser, mais avec une couverture, partielle ou complète. La production régulière de MNTs permet d'exprimer un changement de volume correspondant à la somme de toutes les variations de volume observées à l'échelle d'un pixel durant la période étudiée, comme illustré Fig. 8 sur le Mont

Rainier pour la période 2000-2009, 2009-2018et 2000-2018[START_REF] Menounos | Heterogeneous Changes in Western North American Glaciers Linked to Decadal Variability in Zonal Wind Strength[END_REF]. Afin de calculer le SMB, il est nécessaire de faire une hypothèse sur la densité de surface du glacier. Il est pour cela courant d'utiliser une valeur moyenne représentative de la densité de surface, dont l'estimation faite par Huss, (2013) avec surf  =850±60 kg m -3 est la plus utilisée. Les SMB obtenus en utilisant la méthode géodésique peuvent ainsi être exprimés en m w.e.

Historiquement, cette méthode provient de la photogrammétrie terrestre (e.g. [START_REF] Lang | Die Volumenänderungen des Hintereisferners (Ötztaler Alpen) im Vergleich zur Massenänderung im Zeitraum 1953-64[END_REF] ou aéroportée (e.g. [START_REF] Finsterwalder | Zur Höhenänderung von Ostalpengletschern im Zeitraum 1969-1979[END_REF][START_REF] Vincent | Dynamic behaviour analysis of glacier de Saint Sorlin, France, from 40 years of observations, 1957-97[END_REF], et est toujours utilisée actuellement. La méthode géodésique comporte plusieurs sources d'erreurs :

1. Les erreurs liées au capteur, à savoir sa résolution radiométrique et sa possible saturation sur des surfaces très réfléchissantes comme la neige en zone d'accumulation [START_REF] Raup | Generation of data acquisition requests for the ASTER satellite instrument for monitoring a globally distributed target: glaciers[END_REF], l'encodage des données ou encore les perturbations du signal radiométrique dans l'atmosphère.

2. Les erreurs liées à la production des MNTs, telles que le géoréférencement, la projection ou reprojection des images. La coregistration et l'orthorectification des images peut aussi induire des erreurs sur les cartes de variations d'épaisseur obtenues même si des méthodes existent pour réduire leurs impacts (e.g., Nuth and Kääb, 2011). [START_REF] Thibert | Glaciological and volumetric mass-balance measurements: error analysis over 51 years for Glacier de Sarennes, French Alps[END_REF][START_REF] Zemp | Reanalysing glacier mass balance measurement series[END_REF]méthode ELA, Rabatel et al., 2005, 2016). [START_REF] Bamber | A review of remote sensing methods for glacier mass balance determination[END_REF]. Il existe en revanche des méthodes permettant, à partir d'observations de longueurs ou de surfaces glaciaires, d'estimer un volume de glace. Plusieurs observations répétées peuvent ainsi permettre de faire une estimation du SMB sur une ou plusieurs années [START_REF] Bahr | The physical basis of glacier volume-area scaling[END_REF]. Cette méthode, nommée scaling theory, se base sur l'hypothèse que le volume du glacier (V), la surface (S) et la longueur (L) sont liés par une loi de type puissance, telle que :

 L S  (3)  S V  (4)
où  et  sont des coefficients établis empiriquement (Bahr, 1997). Cette méthode repose sur un modèle de continuité [START_REF] Nye | The response of glaciers and ice-sheets to seasonal and climatic changes[END_REF] se basant sur la réponse dynamique du glacier à une perturbation de son SMB pouvant engendrer une variation de longueur du glacier [START_REF] Hoelzle | Secular glacier mass balances derived from cumulative glacier length changes[END_REF]. Le temps de réponse associé du glacier est exprimé par l'équation 5 (Jóhannesson et al., 1989a(Jóhannesson et al., , 1989b)): [START_REF] Bahr | Width and length scaling of glaciers[END_REF]. Ces mesures étant complexes et coûteuses, elles sont rares et les coefficients ainsi dérivés sont parfois peu significatifs [START_REF] Bahr | Characteristic massbalance scaling with valley glacier size[END_REF]. De plus, les coefficients obtenus permettant d'estimer le volume du glacier à partir de sa longueur suggèrent une forme en «V» quasi systématique, ce que l'on sait aujourd'hui être peu probable [START_REF] Bahr | Width and length scaling of glaciers[END_REF]. Des méthodes similaires sont davantage utilisées afin d'estimer le volume de glace global que pour estimer une variation de volume (Grinsted, 2013;[START_REF] Huss | Distributed ice thickness and volume of all glaciers around the globe[END_REF]; Ice Thickness Models Intercomparison eXperiment (ITMIX) : Farinotti et al., 2017[START_REF] Farinotti | A consensus estimate for the ice thickness distribution of all glaciers on Earth[END_REF]. De plus, les variations de front peuvent être difficilement observables par télédétection chaque année et parfois complexes et erronées dans les zones couvertes de débris [START_REF] Paul | On the accuracy of glacier outlines derived from remote-sensing data[END_REF], indépendamment de la résolution du capteur [START_REF] Paul | Comparison of TM-derived glacier areas with higher resolution data sets[END_REF].

snoot a b h max   (5)
En conclusion, les mesures de variation de front peuvent représenter un indicateur climatique intéressant pour de longues périodes de temps (i.e.

plusieurs décennies, siècles), supérieures au temps de réponse des glaciers.

Mais malgré une étude ayant permis de quantifier le bilan d'environ 1000 glaciers sur le globe à partir d'observation de variations géométriques [START_REF] Hoelzle | Secular glacier mass balances derived from cumulative glacier length changes[END_REF] (faible AAR donc faible  ), alors qu'une année à bilan positif mènera à une [START_REF] Dumont | Monitoring spatial and temporal variations of surface albedo on Saint Sorlin Glacier (French Alps) using terrestrial photography[END_REF]Sirguey, 2009). Plusieurs études ont ensuite permis de valider le lien entre le SMB annuel des glaciers de montagne et  à l'échelle de glaciers individuels [START_REF] Colgan | Surface albedo as a proxy for the mass balance of Greenland's terrestrial ice[END_REF]Dumont et al., 2012;Sirguey et al., 2016;deux glaciers pour Brun et al., 2015;[START_REF] Wang | Spatial and temporal variations of albedo on nine glaciers in western China from 2000 to 2011[END_REF] trois glaciers pour [START_REF] Zhang | Annual Glacier-Wide Mass Balance (2000-2016) of the Interior Tibetan Plateau Reconstructed from MODIS Albedo Products[END_REF], ainsi qu'à l'échelle saisonnière (Sirguey et al., 2016).

Une limitation de cette méthode repose sur sa dépendance à des données de calibration obtenues par mesures in situ, véritable frein pour quantifier le SMB de glaciers ou de chaînes de montagne où pas ou peu de mesures existent. Une piste permettant de s'affranchir des mesures in situ a été présentée par Rabatel et al. (2017) [START_REF] Aa | Topographic control of equilibrium-line altitude depression on reconstructed "Little Ice Age" glaciers, Grovabreen, western Norway[END_REF][START_REF] Benn | Mass balance and equilibrium-line altitudes of glaciers in high-mountain environments[END_REF][START_REF] Gross | Methodische Untersuchungen über die Schneegrenze in alpinen Gletschergebieten[END_REF][START_REF] Torsnes | Modern and Little Ice Age equilibrium-line altitudes on Outlet Valley glaciers from Jostedalsbreen, western Norway: An evaluation of different approaches to their calculation[END_REF]. En utilisant la méthode de l'AAR, par exemple, il est possible, connaissant l'ancienne extension glaciaire de part une analyse géomorphologique (i.e. érosion, moraines, lichens), d'en déduire la ELA en se basant sur une hypothèse de l'AAR. La ELA étant un très bon indicateur climatique, son estimation permet, après application d'un modèle, de pouvoir quantifier des écarts de températures et de précipitations par rapport à aujourd'hui (e.g. [START_REF] Dahl | Paleoclimatic implications based on equilibrium-line altitude depressions of reconstructed Younger Dryas and Holocene cirque glaciers in inner Nordfjord, western Norway[END_REF][START_REF] Hughes | Little Ice Age glaciers in the Balkans: low altitude glaciation enabled by cooler temperatures and local topoclimatic controls[END_REF][START_REF] Klein | Modern and last local glacial maximum snowlines in the Central Andes of Peru, Bolivia, and Northern Chile[END_REF]Rabatel et al., 2008a).

Méthodes 31

La ELA est de plus fortement corrélée avec le SMB dont elle est un très bon indicateur (Braithwaite, 1984;[START_REF] Kuhn | The Response of the Equilibrium Line Altitude to Climate Fluctuations: Theory and Observations, Glacier Fluct[END_REF]. L'arrivée des images aéroportées puis des images satellites a permis d'observer la ligne de neige et de l'utiliser comme marqueur des variations de volume du glacier (LaChapelle, 1962;Meier and Post, 1962;[START_REF] Østrem | Erts Data in Glaciology-An Effort to Monitor Glacier Mass Balance from Satellite Imagery[END_REF]. n'étant pas sujets à la formation de glace surimposée [START_REF] Cogley | Glossary of glacier mass balance and related terms, IHP-VII technical documents in hydrology[END_REF], il est communément admis que la SLA est représentative de la ELA [START_REF] Lliboutry | Traité de glaciologie[END_REF]. (Braithwaite, 1984;Rabatel et al., 2005[START_REF] Rabatel | Can the snowline be used as an indicator of the equilibrium line and mass balance for glaciers in the outer tropics ?[END_REF] Benn and Lemkhul, (2000). [START_REF] Chandrasekharan | Quantification of annual glacier surface mass balance for the Chhota Shigri Glacier, Western Himalayas, India using an Equilibrium-Line Altitude (ELA) based approach[END_REF][START_REF] Kuhn | Mass Budget Imbalances as Criterion for a Climatic Classification of Glaciers[END_REF][START_REF] Mayo | Glacier Mass Balance and Runoff Research in the U. S. A[END_REF]Rabatel et al., 2005).

Les équations de cette méthode ainsi qu'une explication plus complète sont présentées dans l'article Rabatel et al. (2005) ainsi que dans la Sect. 4.

En utilisant la méthode décrite par Braithwaite, Dedieu [START_REF] Dedieu | Apport des données SPOT et Landsat TM pour le suivi de la fusion nivale et des bilans glaciaires dans les Alpes françaises[END_REF] puis Rabatel purent reconstruire le bilan de masse annuel de glaciers dans les Alpes françaises [START_REF] Rabatel | Reconstitution des fluctuations du bilan de masse du Glacier Blanc (Massif des Ecrins, France) entre 1985 et 2000, par télédétection optique (imagerie Spot et Landsat)[END_REF](Rabatel et al., , 2005(Rabatel et al., , 2008b[START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF], dans les Andes tropicales [START_REF] Rabatel | Can the snowline be used as an indicator of the equilibrium line and mass balance for glaciers in the outer tropics ?[END_REF] ainsi que [START_REF] Chandrasekharan | Quantification of annual glacier surface mass balance for the Chhota Shigri Glacier, Western Himalayas, India using an Equilibrium-Line Altitude (ELA) based approach[END_REF] sur un glacier himalayen, sans calibration nécessaire avec des mesures in situ excepté pour le gradient de bilan de masse. Cette méthode permet de calculer le bilan de masse de glaciers pour lesquels aucune donnée in situ n'a été recueillie.

D'autres approches se développèrent récemment en parallèle, ne se basant pas sur la méthodologie développée par Braithwaite, mais utilisant par exemple l'AAR comme indicateur du bilan de masse de surface (e.g. Dyurgerov et al., 2009;[START_REF] Hock | Deriving glacier mass balance from accumulation area ratio on Storglaciären, Sweden[END_REF]Mernild et al., 2013), où la TSL comme marqueur de fonte pendant la saison d'ablation (e.g. Barandun et al., 2018;[START_REF] Chinn | Glacier snow line variations in the Southern Alps, New Zealand[END_REF][START_REF] Kaser | Changes of the equilibrium-line altitude in the tropical Cordillera Blanca, Peru, 1930-50, and their spatial variations[END_REF]LaChapelle, 1962;Meier and Post, 1962;Mernild et al., 2013;[START_REF] Naegeli | Change detection of bare-ice albedo in the Swiss Alps[END_REF][START_REF] Østrem | Erts Data in Glaciology-An Effort to Monitor Glacier Mass Balance from Satellite Imagery[END_REF]Pelto, 2011;Shea et al., 2013;[START_REF] Spiess | MODISderived interannual variability of the equilibrium-line altitude across the Tibetan Plateau[END_REF][START_REF] Tawde | Estimation of Glacier Mass Balance on a Basin Scale:An Approach Based on Satellite-Derived Snowlines and a Temperature Index Model[END_REF] et al., 2002;[START_REF] Echelmeyer | Airborne surface profiling of Alaskan glaciers[END_REF].

La méthode géodésique permet, quant à elle, l'estimation du SMB à l'échelle régionale, tout en ayant une résolution adéquate pour quantifier le SMB de glaciers de petites tailles (> 0.5 km2, [START_REF] Menounos | Heterogeneous Changes in Western North American Glaciers Linked to Decadal Variability in Zonal Wind Strength[END_REF]. De nouvelles chaînes de traitement automatisées (e.g., Ames Stereo Pipeline ASP; [START_REF] Beyer | The Ames Stereo Pipeline: NASA's Open Source Software for Deriving and Processing Terrain Data[END_REF]Shean et al., 2016) Rabatel (e.g., Braithwaite, 1984 ;Rabatel et al. 2005Rabatel et al. , 2008Rabatel et al. , 2016)). Cette méthode, exceptée pour l'estimation du gradient de SMB, ne nécessite aucune calibration issue de données in situ. Le principale verrou d'une application de cette approche à une échelle régionale et multitemporelle réside dans l'automatisation de la détection de la SLA sur les images satellites même si une étude récente propose une méthode automatique pour déterminer la SLA, jamais encore appliquée à l'échelle régionale [START_REF] Rastner | On the Automated Mapping of Snow Cover on Glaciers and Calculation of Snow Line Altitudes from Multi-Temporal Landsat Data[END_REF]. Cette détection doit ensuite être combinée avec l'estimation du SMB pluri-décennal calculé grâce à la méthode géodésique. Ce jeu de données unique nous a permis de poursuivre des travaux menés précédemment sur un petit nombre de glaciers (e.g. [START_REF] Huss | Distributed ice thickness and volume of all glaciers around the globe[END_REF]Rabatel et al., 2013[START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF], dans le but d'étudier le lien entre les variables morpho-topographiques et le SMB. Une analyse des bilans moyens pluriannuels a permis de conforter (i.e. [START_REF] Brun | Heterogeneous Influence of Glacier Morphology on the Mass Balance Variability in High Mountain Asia[END_REF][START_REF] Fischer | Surface elevation and mass changes of all Swiss glaciers 1980-2010[END_REF] Les Andes tropicales sont le plus grand réservoir de glace des tropiques, lequel est essentiel pour la ressource en eau locale, aussi bien pour l'eau potable, l'agriculture, la production hydro-électrique, ou les aléas glaciaires (e.g. [START_REF] Bury | Glacier recession and human vulnerability in the Yanamarey watershed of the Cordillera Blanca, Peru[END_REF][START_REF] Mark | Climate Change and Tropical Andean Glacier Recession: Evaluating Hydrologic Changes and Livelihood Vulnerability in the Cordillera Blanca, Peru[END_REF]. Dans le contexte climatique actuel, une importante récession est observée sur les glaciers des Andes tropicales (e.g. Rabatel et al., 2013;[START_REF] Vuille | Climate change and tropical Andean glaciers: Past, present and future[END_REF], pouvant avoir des conséquences sur la saisonnalité et leurs contributions à l'hydrologie des bassins versants (Baraer et al., 2012;[START_REF] Huss | Global-scale hydrological response to future glacier mass loss[END_REF]. En revanche, peu de données existent pour quantifier les variations de masse des glaciers dans cette région, et les incertitudes restent importantes (e.g. [START_REF] Zemp | Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016[END_REF].

Bénéficiant de données de SMB in situ ainsi que de lignes de neiges délimitées manuellement dans la Cordillère Royale et Blanche [START_REF] Loarte | Determinación de la variación espacio-temporal de la altura de la línea de equilibrio de los glaciares de la Cordillera Blanca, Perú[END_REF][START_REF] Rabatel | Can the snowline be used as an indicator of the equilibrium line and mass balance for glaciers in the outer tropics ?[END_REF]Soruco et al., 2009), nous avons pu utiliser l'algorithme de détection semi-automatique de la ligne de neige et la méthode associée pour estimer le SMB annuel dans cette région, malgré une saisonnalité très différente des Alpes européennes [START_REF] Rabatel | Can the snowline be used as an indicator of the equilibrium line and mass balance for glaciers in the outer tropics ?[END_REF], Univ. Grenoble Alpes, CNRS, IRD, Institut des Géosciences de l'Environnement (IGE), F-38000 Grenoble, France; antoine.rabatel@univ-grenoble-alpes.fr; yves.arnaud@ird.fr; lucas.davaze@univ-grenoble-alpes.fr; jean-pierre.dedieu@univgrenoble-alpes.fr Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 61 Abstract: Glaciers are one of the terrestrial essential climate variables (ECVs) as they respond very sensitively to climate change. A key driver of their response is the glacier surface mass balance that is typically derived from field measurements. It deserves to be quantified over long time scales to better understand the accumulation and ablation processes at the glacier surface and their relationships with inter-annual changes in meteorological conditions and long term climate changes. Glaciers with in situ monitoring of surface mass balance are scarce at global scale, and satellite remote sensing provide a powerful tool to increase the number of monitored glaciers. In this study, we present a review of three optical remote sensing methods developed to quantify seasonal and annual glacier surface mass balances. These methodologies rely on the multitemporal monitoring of the end-of-summer snow line for the equilibrium-line altitude (ELA) method, the annual cycle of glacier surface albedo for the albedo method, and the mapping of the regional snow cover at seasonal scale for the snow-map method. Together with a presentation of each method, an application is illustrated. The ELA method shows promising results to quantify annual surface mass balance and to reconstruct multi-decadal time series. The other two methods currently need a calibration on the basis of existing in situ data; however a generalization of these methods (without calibration) could be achieved. The two latter methods show satisfying results at the annual and seasonal scales, particularly for the summer surface mass balance in the case of the albedo method and for the winter surface mass balance in the case of the snow-map method. The limits of each method (e.g., cloud coverage, debris-covered glaciers, monsoonregime and cold glaciers), their complementarities and the future challenges (e.g., automating of the satellite images processing, generalization of the methods needing calibration) are also discussed.

région

INTRODUCTION

Glaciers and ice caps are one of the 50 essential climate variables recognized by the Global Climate Observing System (GCOS), an internationally coordinated network of observing systems designed to meet evolving national and international requirements for climate observations [START_REF][END_REF]. Among the main glacier characteristics (e.g., length, surface area, mass balance), the surface mass balance (SMB) is the most relevant in terms of climate as it is governed by processes of accumulation and ablation at the glacier surface, primarily driven by the atmospheric conditions.

Logistical, human and financial constraints associated with field measurements result in a comparably small number of glaciers being permanently monitored in different climate regions worldwide. Regarding the SMB, the World Glacier Monitoring Service (WGMS) lists data available for 440 glaciers monitored worldwide through direct field measurements, including uninterrupted series longer than 40 years, for only 40 glaciers [2]. The 440 monitored glaciers represent only 0.25% of the nearly 200,000 glaciers inventoried on Earth and included in the Randolph Glacier Inventory (RGI) and GLIMS databases [3][4].

This small sample limits our understanding of the relationship between climate and glacier changes at regional to global scales, as well as our understanding of the contribution of glaciers to water resources and to the biodiversity of high-altitude watersheds.

Consequently, there is a need for monitoring approaches that are complementary to field measurements to quantify the SMB of individual glaciers at a regional scale. For this purpose, satellite remote sensing has proven to be a very useful tool. The increase in the number and availability of satellite sensors over nearly four decades, as well as their improved performance, i.e., their improved spatial, radiometric and temporal resolutions, has made remote sensing approaches vital for the monitoring of glaciers. Many remote sensing methods have been developed since the 1970s to determine Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 63 the physical characteristics of glaciers [e.g., 5]. While variables such as length and areal extent, as well as their variations, could be retrieved relatively easily from sparse or intermittent imagery, the increase of satellite stereo imagery and advances in space-borne photogrammetry have allowed glacier topography to be repeatedly mapped, and, in turn, geodetic mass balance to be derived. For example, digital elevation models (DEMs) can be derived from medium-to-high resolution optical satellite images (e.g., SPOT-5/HRS or ASTER) images or radar satellite data (e.g., SRTM) but the vertical precision (~5 m) prevents an accurate estimation of the mass balance of individual glaciers typically smaller than ~10 km². New generations of satellites submetric optical stereo-images (e.g. Worldview, Pléiades) show promises to provide accurate glacier topography with a vertical precision of ±1 m and even ±0.5 m on gently-sloping glacier tongues and retrieve annual and seasonal glacier-wide SMB using the DEM differencing method [e.g., [6][7]. However, the deployment of this method at a regional scale is still hampered by the limited availability, footprint (e.g., 20 km x 20 km for Pléiades) and the cost of these commercial stereo-images at the right time of the year. In addition, it is worth noting that the geodetic method needs a density value for the conversion from volume change to mass change. Because this value and its spatio-temporal variations at the glacier scale are hardly quantifiable, this method currently provides mass balance at a temporal resolution typically of five to ten years.

Alternatively, changes in the state of a glacier surface have been used to infer its surface mass balance. In particular, variations in the respective area of the accumulation zone and the ablation zone of glaciers via the mapping of the end-of-summer snow line support estimates of the equilibrium-line altitude (ELA), and/or the accumulation area ratio (AAR) that are used in turn to infer annual SMB. In this context, improvements and systematic application of the ELA method, as well as emerging approaches based on the regular monitoring of snow and ice surfaces have confirmed the potential to confidently estimate annual and seasonal SMB of glaciers [e.g., 8-9]. In this study, we present an overview of three remote sensing techniques developed recently to retrieve the annual and seasonal SMB of mountain glaciers using optical satellite images of high (SPOT, Landsat, ASTER), medium (MODIS) and low (SPOT/VGT, Proba-V) resolution. These methods share a common theme that is to characterize the state of the glacier surface (i.e., snow line altitude, glacier surface albedo, seasonal snow cover) and its spatio-temporal changes.

Each of the three methods is presented in section 4.2 with an application. Section 4.3 discusses the limits and complementarities of the methods as well as the future challenges together with possible further improvements.

METHODS TO RETRIEVE ANNUAL AND SEASONAL GLACIER-WIDE SURFACE MASS BALANCES USING OPTICAL SATELLITE IMAGERY

In order to simplify the text, the method using the snow line detection as a proxy of the ELA from high resolution optical satellite images is called the "ELA method"; the one that uses the quantification of the glacier-wide surface albedo from medium resolution optical satellite images is called the "albedo method"; and the one that uses the mapping of the average regional snow cover extent from low resolution optical satellite images is called the "snowmap method".

The ELA method

Principles and history of the method For glaciers where the formation of superimposed ice [10] is negligible, it is commonly accepted that the position of the glacier snow line at the end of the ablation season is representative of the equilibrium line [11]. Indeed, in the glacier reaches located above the altitude of the end-of-summer snow line, accumulation is greater than the ablation and this part therefore constitutes the accumulation zone. Conversely, in the glacier reaches located below the altitude of the end-of-summer snow line, ablation is greater than the Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 65 accumulation and this part therefore constitutes the ablation zone. Thus, the SMB at the altitude of the end-of-summer snow line is zero, and the snow line therefore matches the equilibrium line. Moreover, at annual scale, the altitude of the equilibrium line (ELA) exhibits high correlation with the SMB [e.g., 12-14].

The potential for mapping the end-of-summer snow line in glaciological studies was considered as early as the 1950-60s based on aerial photographs [15][16],

then since the early 1970s using the first civilian optical satellite imagery [17].

In the early 1980s, Braithwaite used in situ measurements of glacier SMB and ELA to propose a method to quantify the glacier SMB from its equilibrium line [12]. The method was based on the altitude difference between the ELA for a given year and the balance-budget ELA (called ELA0 and assumed to be constant by Braithwaite). To compute the annual SMB of the considered year, this difference is multiplied by the "effective balance gradient" which represents a time-and space-average of the SMB gradient with altitude; a value considered constant by Braithwaite for a given region. Braithwaite yet recognized that there are large regional variations of the "effective balance gradient" and pointed out the potential error of using a universal value.

Moreover, Braithwaite postulated that this value is representative of the mass balance gradient. However, when compared with measured SMB gradients available today, one can note that Braithwaite's "effective balance gradient" values were 2 to 7 times higher than the SMB gradients.

Building upon these pioneering works by Braithwaite, a method to quantify the glacier SMB from the equilibrium line has been developed with the idea to be applied with remote-sensing data only, thus without any reliance on in situ measurements [13,[18][19][20]. This method has been validated for mid-latitudes alpine glaciers [13,20], and glaciers in the outer tropics [21].

It is worth noting that other studies took advantage of the detection of the end-of-summer snow line from airborne or satellite data to document its yearto-year variability [22] or the changes in altitude of the snow line over the summer period as a proxy of ablation [23]. Other studies quantified the annual glacier wide SMB by combining the altitude of end-of-summer snow line with the glacier hypsography and the SMB gradients above and below the ELA [24],

or with a temperature-index model [25]. The method proposed by Shea et al.

[24] depends on the availability of the SMB gradients in the ablation and accumulation zones, and uncertainties associated with these gradients and their use to compute the glacier-wide SMB can lead to significant errors. On the other hand, the method proposed by Tawde et al. [25] uses meteorological data and therefore depends on their availability in the vicinity of the studied glaciers. In this method, the errors related to the extrapolation of the meteorological variables can significantly affect the computed SMB.

The strength of the method proposed by Rabatel et al. [13] is due to the combination of the geodetic method and the temporal changes in the ELA. The geodetic method using DEM differencing ensures an accurate quantification of the changes in glacier mass over the entire study period (i.e. the value of the glacier-wide average annual SMB), and the changes in the ELA allow computing the annual departure in the glacier-wide SMB from its average computed by the geodetic method. In addition, in Rabatel et al. [13] the value used for the SMB gradient is the one in the vicinity of the ELA, i.e. in an elevation range where assuming this gradient as a constant over time and space is reasonable. Note that this assumption is also supported by Mernild et al. [23].

Consequently, in the method described in Rabatel et al. [13], the glacier-wide annual SMB series is reconstructed in two steps. First, for each year, the difference between altitude of the snow line at the end of the ablation season, SLAi (in m, used as a surrogate of ELAi), and the theoretical altitude of the equilibrium line if the glacier had a balanced budget (glacier wide mass balance = 0) over the study period, ELAeq (in m), is calculated. ELAeq is computed as:

Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 67

       n i eq z b B SLA n ELA 1 / 1 , ( 1
)
where B is the average glacier-wide annual SMB over the study period (in m w.e. yr -1 ) quantified by the DEM differencing method, and ∂b/∂z is the SMB gradient in the vicinity of the ELA. This gradient, estimated from in situ data, is assumed constant across a glacierized region, for example 0.78 m w.e. 100 m -1 for the French Alps according to Rabatel et al. [13]. Rabatel et al. [13] demonstrated that even if ∂b/∂z for a specific glacier departs slightly from this value, using an average value at a regional scale does not compromise the results. The glacier-wide annual SMB Ba (in m w.e. yr -1 ), computed for each year over the study period can then be expressed as:

z b SLA ELA B i eq a      ) ( . (2)
Applying the ELA method requires the following:

 Two DEMs covering the entire region of interest (RoI); one at (or close to) the beginning of the study period, one at (or close to) the end.

 One cloud-free satellite image for each glacier under study and for each year of the study period. In temperate latitudes, the image must have been acquired at the end of the summer season, and in the outer tropics, during the dry season.

 An estimate of the mass-balance gradient in the vicinity of the ELA.

Application of the method on 30 glaciers in the French Alps

In a recent study, Rabatel et al. [9] Rabatel et al.

[9] also considered the impact of using an intermediate DEM when the study period is long (e.g. several decades). Indeed, to quantify the average annual SMB over the entire period by DEM differencing, the mean surface area of the glacier is generally taken as the arithmetic mean of the initial and final areas. This assumes that surface area and SMB trends are linear, which is not necessarily true. Accordingly, the authors concluded that for study period of several decades with potential periods of glacier mass gain and loss, using an intermediate DEM is recommended to reduce errors related to non-linear changes in surface area.

The albedo method

Principles and history of the method

The albedo, α, of a surface measures the ratio of solar radiation being reflected by this surface. In the case of snow and ice targets, albedo often controls the surface energy budget that in turn governs glacier ablation [e.g., 28-29]. On temperate glaciers, this control accentuates during the ablation season as relatively high reflective snow (α 0.8) progressively gives way to more absorptive ice surfaces (α 0.4 or less), concurrently with a significant increase in daily shortwave radiation [30][31]. The rationale behind the albedo method thus stemmed from the capacity of glacier surface albedo to integrate processes of accumulation (e.g., albedo rising with fresh snowfall) and ablation

Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 71 (decrease of albedo as snow transforms and/or glacier is increasingly exposed) combined with its significant physical control on net energy fluxes.

In attempting to estimate the mass balance of Vatnajökull using AVHRR data, pioneers of this approach De Ruyter De Wildt et al. [32] faced difficulties in resolving the equilibrium-line due to ambiguities in the transition between ice, firn, and snow, further complicated by cloudiness. Having converted spectral reflectance to broadband albedo, these authors averaged the albedo over the glacier surface, hereby named glacier-wide albedo, 𝛼 ̅. The sum of absorptance (i.e., 1-𝛼 ̅) weighted by the extraterrestrial solar radiation and integrated over a summer was postulated and found to correlate well with the surface mass balance. Greuell and Oerlemans [33] confirmed and refined this approach on a transect in Greenland by using the satellite-derived albedo in an approximation of the energy budget instead of an empirical regression. The method was then applied more widely by Greuell et al. [8] to resolve interannual anomalies in SMB of 18 Svalbard glaciers found to be consistent with observations. Nonetheless, Greuell et al. [8] stressed discrepancies between observations and satellite estimates that impair the usefulness of the method when not calibrated by absolute values of SMB. Dumont et al. [34] revisited this approach by building upon progress in the atmospheric and topographic corrections of optical imagery [35], combined with an improved retrieval method of snow and ice albedo [36]. This allowed time series of glacier surface albedo of mountain glaciers to be derived from MODIS data. Under the simplest assumption that, at the end of the ablation season, a glacier of surface 𝐴 , accounts for a relative share between an accumulation zone of area 𝐴 𝑎𝑐𝑐 , covered by snow with an average albedo of 𝛼 𝑠𝑛𝑜𝑤 , and an ablation zone of area 𝐴 𝑎𝑏𝑙 , exhibiting exposed glacier ice with an average albedo of 𝛼 𝑖𝑐𝑒 , then, assuming a linear mixing of the reflectance, the glacier-wide albedo can be obtained as

𝛼 ̅ = 𝛼 𝑠𝑛𝑜𝑤 𝐴 𝑎𝑐𝑐 +𝛼 𝑖𝑐𝑒 𝐴 𝑎𝑏𝑙 𝐴 . ( 3
) Using 𝐴𝐴𝑅 = 𝐴 𝑎𝑐𝑐 𝐴 = 1 - 𝐴 𝑎𝑏𝑙 𝐴
, it comes that:

𝛼 ̅ = 𝛼 𝑖𝑐𝑒 + 𝐴𝐴𝑅(𝛼 𝑠𝑛𝑜𝑤 -𝛼 𝑖𝑐𝑒 ). (4)
In other words, the glacier-wide albedo is a case of linear spectral mixing that approaches directly the AAR with the sensitivity of 𝛼 ̅ to AAR being controlled by the contrast in albedo between the exposed ice and the snow at the end of the ablation season. Note that this also requires that the albedo of snow and ice at the end of the ablation season remain constant from one year to another.

By leveraging frequent MODIS observations, Dumont et al. [34] used the minimum reached by 𝛼 ̅ at the end of the ablation season, namely 𝛼 ̅ min to capture the AAR. In turn, and consistently with studies showing AAR to be a good predictor of SMB [e.g., 14,37-38], Dumont et al. [34] showed that 𝛼 ̅ min could be used as a surrogate for the annual SMB on Saint-Sorlin glacier (French Alps) using a simple linear model, namely

𝐵 𝑎 = 𝑠 * 𝛼 ̅ min + 𝑜, (5)
where s and o are the slope and offset of the linear regression, respectively.

Using a similar approach, Wang et al. [39] also obtained significant correlations between 𝛼 ̅ at the end of summer and annual SMB on two glaciers in western China. Brun et al. [40] further tested the method on Chhota Shigri (India) and Mera (Nepal) glaciers in Himalaya allowing annual SMB series to be reconstructed over the period of MODIS observations. The correlation between 𝛼 ̅ min and 𝐵 𝑎 on Mera Glacier was found to be less robust because of persistent cloud cover and snow accumulation during the monsoon season that prevented the seasonal variations of 𝛼 ̅ to be fully resolved despite near-daily imaging with MODIS. As a summer accumulation type glacier, Mera Glacier also exhibited a relatively narrow range of 𝛼 ̅ min that contributed to a weaker linear model. A linear correlation between surface mass balance and albedo over Greenland's terrestrial ice was also shown by Colgan et al. [41] and used to estimate monthly surface mass balance. Finally Sirguey et al. [42] successfully applied the albedo method on Brewster Glacier in New Zealand, Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 73 allowing the relatively short record of in situ measurements to be reconstructed over 37 years.

Application of the albedo method on New Zealand Southern Alps glaciers Only a few direct SMB studies have been made in the Southern Alps of New Zealand, bringing relatively short term observations [START_REF] Chinn | ANNUAL ICE VOLUME CHANGES 1976-2008 FOR THE NEW ZEALAND SOUTHERN ALPS[END_REF]. In lieu of direct The study documented the relationship between AAR, ELA, and SMB, and supported the validity of linear approximations within range of variations of ELA or AAR. Sirguey et al. [42] confirmed the potential of the albedo method in New Zealand with 𝛼 ̅ min found to be a suitable proxy to estimate both the annual and summer SMB (Fig. 3). It enabled the SMB series of Brewster Glacier to be reconstructed over the period of MODIS acquisition and confirmed the potential of frequent monitoring of glacier surface properties as a method capable of sustained systematic observations of glacier SMB signal. Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 75 Therefore, the albedo method can capture a SMB signal comparatively relevant to that of the snow line aerial survey program (analog to the ELA method).

Moreover, the albedo method has the advantage of being available from satellite remote sensing and a potential of generalization to monitor an entire glacierized region with less logistical burden. However, the conversion of 𝛼 ̅ min to SMB values has, so far, been approached only via calibration with in situ data. It is possible to propose a development of the albedo method inspired from the ELA method to quantify the SMB readily from 𝛼 ̅ min , namely a gradientbased albedo method. First, we note that the SMB-elevation gradient 𝜕𝑏/𝜕𝑧, used by Rabatel et al. [13] can be converted into a SMB-area gradient 𝜕𝑏/𝜕𝐴𝐴𝑅, using the hypsometric curve of a glacier using:

𝜕𝑏 𝜕𝐴𝐴𝑅 = 𝜕𝑏 𝜕𝑧 × 𝜕𝑧 𝜕𝐴𝐴𝑅 . (6)
For example on Brewster Glacier, Cullen et al. [38] found that equation ( 4) yields estimates 𝐴𝐴 ̂𝑅𝑖 that match the values obtained from reanalysis of in situ measurements (Fig. 6b). With regards to the ELA method, the AAR of a glacier in equilibrium can be estimated by:

𝐴𝐴𝑅 𝑒𝑞 = 1 𝑛 ∑ 𝐴𝐴 ̂𝑅𝑖 + 𝐵 ̅ 𝜕𝑏 𝜕𝐴𝐴𝑅 . ( 7
)
Using the values above, we find 𝐴𝐴𝑅 𝑒𝑞 = 47% for Brewster Glacier which is close to the estimate based on the hypsometric curve (46%) reported Cullen et al. [38]. Finally, yearly values of surface mass balance 𝐵 𝑎 , can be estimated by:

𝐵 𝑎 = (𝐴𝐴𝑅 𝑒𝑞 -𝐴𝐴 ̂𝑅𝑖 ) 𝜕𝑏 𝜕𝐴𝐴𝑅 . (8)
Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 77 Figure 6c illustrates the good agreement between SMB values of Brewster Glacier derived by the gradient-based albedo method and those initially obtained via direct regression with in situ measurements. Although this example only serves as a proof of concept, it shows that the approach of the ELA method can be generalized to the albedo method successfully. Similar to Rabatel et al. [24] and Chinn et al. [START_REF] Chinn | ANNUAL ICE VOLUME CHANGES 1976-2008 FOR THE NEW ZEALAND SOUTHERN ALPS[END_REF] who postulated that 𝜕𝑏 𝜕𝑧 could be extrapolated to glacier (sub-)regions to apply the ELA method, a similar assumption could be made along with a specific glacier hypsometric curve to extrapolate the gradient-based albedo method to glaciers with no in situ measurements.

Figure 6: (a) Hypsometric curve of Brewster Glacier. (b) Comparison between AAR estimated from a linear spectral mixture model of glacier-wide surface albedo and reference values from reanalysis of in situ glaciological observations. (c) Comparison between mass balance obtained with the gradient-based albedo method and estimated from the re-analyzed glaciological method (2005-2013, [38]) and calibrated albedo method (2000-2004, [42]).

Finally, the albedo method proved to support more than the annual balance only. Sirguey et al. [42] hypothesizes that, in addition to 𝐵 𝑎 , the temporal variation of 𝛼 ̅ could also support a signal related to summer and winter mass balance 𝐵 𝑠 and 𝐵 𝑤 , respectively. 𝛼 ̅ min was found to correlate well with 𝐵 𝑠 which is consistent with glaciers whereby the variability of annual SMB is predominantly controlled by the variability of summer SMB [38,44].

Furthermore, Sirguey et al. [42] proposed the cumulative winter albedo 𝐴 𝑤 , namely the integration of 𝛼 ̅ over the accumulation season when 𝛼 ̅ exceeds a threshold 𝛼 𝑡ℎ large enough to capture conditions favorable to winter balance:

𝐴 𝑤 = ∫ { 𝛼 ̅(𝑡), 𝛼 ̅(𝑡) ≥ 𝛼 𝑡ℎ 0, o/w 𝑑𝑡 𝑤𝑖𝑛𝑡𝑒𝑟 . (9)
A significant correlation was found between winter SMB and 𝐴 𝑤 that allowed 𝐵 𝑤 to be estimated using a linear model. The concurrent evolution between 𝐴 𝑤 and 𝐵 𝑤 illustrated in Figure 7 supports that a signal of winter balance can be captured by monitoring variables of the glacier surface that integrate events that are favorable to winter accumulation. This result resonates with the following section of the paper whereby the inter-annual variability of the snow cover around glaciers during summer and winter is also shown to explain variability of the seasonal SMB.

Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 79

The snow-map method Principles and history of the method

This method is based on regional snow maps from daily optical satellite images (e.g., the daily 1 km resolution SPOT/VGT and PROBA-V images) to quantify glacier seasonal SMB. The method hypothesizes that the seasonal mean snow altitude around a glacier is a proxy for the seasonal glacier SMB.

The first step consists of calculating the normalized difference snow index (NDSI) from the satellite images to produce daily snow maps for the winter and summer seasons of the study area during the period covered by the satellite data (e.g., 1998-2014 for SPOT/VGT, 1 km pixel resolution). The NDSI was introduced by Crane and Anderson [45] and Dozier [46] for the Landsat

𝑁𝐷𝑆𝐼 = ( 𝐵0+𝐵2 2 -𝑆𝑊𝐼𝑅) ( 𝐵𝑂+𝐵2 2 + 𝑆𝑊𝐼𝑅 ⁄ ). (10)
A cloud mask is applied in order to flag cloudy pixels and to avoid overestimating the snow coverage. A temporal interpolation is then computed for the "cloudy" pixels. Maps of winter/summer NDSI are produced by averaging all NDSI syntheses (10-day syntheses for SPOT/VGT images)

included between 1 October and 30 April for each winter, and between 1 May and 30 September for each summer (Fig. 8).

Figure 8: Example of regional seasonal snow maps over the European Alps for two contrasted hydrological years produced by averaging all NDSI syntheses included between 1 October and 30 April for the winter season, and between 1 May and 30 September for the summer season.

Then, the seasonal snow maps are overlaid on a DEM (e.g. SRTM30) to derive the seasonal altitudinal distribution of NDSI in an optimized area surrounding each studied glacier called the windows of snow monitoring (WOSM). The Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 81 inter-annual variability of the altitudinal distribution of NDSI within the WOSM surrounding the glacier is thus obtained for the studied period (16 years for SPOT/VGT; Fig. 9). From the intersection between an optimized NDSI value (specific to each glacier) and each curve of the NDSI seasonal altitudinal distribution, a "mean regional" altitude of snow (Z) can be deduced for summers and winters of the study period.

Finally, a linear regression between the mean regional snow altitudes Z and the seasonal SMB is quantified (Fig. 10). This linear regression allows the estimation of the glacier seasonal SMB, as a function of Z inferred from satellite images (see Eq. ( 11), generalized for both seasons):

𝐵w/s_VGT =∝ w/s * 𝑍w/s + 𝛽w/s. ( 11
)
Bw/s_VGT is the winter/summer SMB estimated for the year y (in m. w.e. yr -1 ), Zw/s is the winter/summer "mean regional" altitude of snow (in m), the slope coefficient αw/s (in mm w.e. yr -1 m -1 ) represents the sensitivity of a glacier winter/summer SMB towards Zw/s, and βw/s is the winter/summer intercept term (in mm w.e. yr -1 ). To apply the snow-map method, the following are needed:

 Several optical satellite images for each season of the studied period, covering the region of interest (ROI).

 One DEM of the ROI acquired during the studied period.

 Observed seasonal SMB available over the studied period covered by the satellite (or for sub-period of the study period).

Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 83 Application of the method on 55 glaciers in the European Alps This method has been applied for 55 glaciers in the European Alps with SPOT/VGT images by Drolon et al. [51]. Promising linear relationships between the regional mean snow altitude and observed seasonal SMB have been found

over the calibration period 1998/1999-2008. The explained variance in winter is statistically significant for all glaciers (R²=0.84 on average) and the averaged RMSE (161 mm w.e. yr -1 ) is below the usual error Eobs associated with glaciological SMB measurements (typically included between ±200 and 400 mm w.e. yr -1 ).

For the summer period, the explained variance is lower although still statistically significant (R²=0.73 on average) with the averaged RMSE (318 mm w.e. yr -1 ) still in the range of Eobs. An individual cross-validation of the 55 linear regressions allowed validating the temporal stability and the robustness of 95% of the relationships in winter and of around 60% of the relationships in summer.

Beside the individual calculation of each glacier SMB presented in Drolon et al.

[51], we illustrate here a synthetic result through the mean reconstructed and observed MB time series averaged for the 55 studied glaciers (Figure 11). In winter (Fig. 11a), <Bw_VGT> is generally in close agreement with <Bw_obs>. The Regarding the snow-map method, the effect of pixels contaminated by undetected clouds (e.g. low-altitude clouds located in the valleys) can partly explain the reduced performance for summer. The error in a pixel NDSI value caused by a cloud is smaller for a snow-covered pixel than for a snow-free pixel. As there is less snow in summer than in winter, undetected clouds may impact NDSI and snow detection more in summer.

Glacier size and elevation

Regarding the ELA method, the main limiting factor for its application is related to the uppermost elevation reached by the glacier. Indeed, the uppermost elevation must be high enough to allow preserving an accumulation zone at the end of the ablation season. The albedo method leverages high frequency acquisitions to resolve the albedo cycle with sufficient details to retrieve minimums representative of the glacier surface condition at the end of the summer. In remote sensing, this typically involves imagery obtained at coarser spatial resolution. For example, the application presented here has relied on near-daily imaging from MODIS with albedo mapped at 250 m spatial resolution. This nominal footprint is compromised of the panoramic distortion of the across track scanner with about a third of observations acquired under view angle greater than 45˚. In effect, the relatively large pixels favor spectral mixtures with nearby terrain that hinder the quality of the albedo signal [42].

Furthermore, the area from which the glacier-wide albedo is calculated needs to account for pixels in both the ablation and accumulation areas to capture contrast and statistically reduce the noise inherent to the albedo retrieval method via averaging. It has so far been applied only to glaciers of area greater than 2 km 2 captured typically by more than 20 MODIS pixels chosen inside the glacier outlines to minimize spectral mixtures. Glaciers on and/or Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 87 surrounded by steep topography can also be exposed to severe shadowing that can compromise the albedo retrieval. This happens typically during winter months when relying on morning MODIS/Terra acquisition. Early afternoon imagery from MODIS/Aqua could be used to mitigate this limitation.

For the snow-map method, the uppermost elevation reached by a glacier or the minimum size of a glacier do not constitute an a priori limit of applicability.

However, in some rare cases, if the regional snow line rises above mountains the snow-map method suffers from same issue as the ELA method and does not allow the estimation of the seasonal SMB (after the calibration period). In those particular cases, the NDSI value optimizing the linear relationships between Z and SMB over the calibration period does not intersect the NDSI altitudinal distribution curves (the curves are below the NDSI value). Thus, no SMB can be retrieved, except if the WOSM size is increased (but at the cost of degrading the regression quality). Moreover, the snow-map method does not seem to perform well for large glaciers (e.g., Aletsch Glacier in the case of the application in the European Alps, see [51]) whose ablation area reaches relatively low elevation. Glacier snout can thus remain snow-free for a considerable part of the period used for winter SMB estimate (October to April).

Winter SMB data for those glaciers thus include both accumulation and melting, whereas the snow-map method for winter is optimized to represent snow accumulation.

Debris covered glaciers

None of the methods presented in this review has been applied on debriscovered glaciers. It is noteworthy that among the glaciers with long-term in situ monitoring of their surface mass balance that are used to calibrate/validate the remote-sensing methods, the large majority are free of debris, or the debris-covered parts are negligible in comparison with the total surface area of the glaciers. Regarding the ELA method, it is assumed that the glacier-wide mass balance is similar to the surface mass balance obtained in the vicinity of the ELAeq. Rabatel et al. [9] underlined that this assumption is reasonable because the mass balances obtain with the ELA method are normalized by the glacier-wide SMB of the entire study period, using satellite DEM differencing. In the case of heavily debris-covered glaciers this assumption deserves to be verified. Indeed, ablation processes on the debriscovered areas differ from the debris-free areas and can strongly influence the glacier-wide SMB [e.g., 54]. Studying a small debris-covered glacier in the Everest region, Vincent et al. [55] indicated that, for this particular glacier, the insulating effect of the debris cover has a larger effect on the total mass loss than the enhanced ice ablation due to supraglacial ponds and exposed ice cliffs.

The calibration of the albedo method has not yet been attempted on glaciers with extensive debris cover, but does prove successful when the area covered by debris is small [40]. As the relative share of the debris-covered area increases, a disconnection between SMB of the debris-free area and the glacier wide SMB may compromise the relationship between AAR and surface mass balance, and in turn the albedo method. Nonetheless, the suggested development of the gradient-based albedo method will benefit from the same assumption and limitation of the ELA method as it will take similar advantage of the average glacier-wide SMB.

Regarding the snow-map method, it is based on the variability of the seasonal regional mean snow altitude to retrieve the SMB. For glaciers with extensive debris cover, the SMB is controlled by other processes than the snow accumulation/ablation. Hence, the snow-map method may be less relevant to retrieve seasonal SMB of debris-covered glaciers.

Polar and monsoon-regime glaciers

The presented methods have been applied on temperate glacier of the midlatitudes, and also on glaciers located within the outer-tropics of the Andes in the case of the ELA method [21].

Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 89

The ELA method requires that the elevation of the snow line at the end of the hydrological year must be representative of the ELA. In the case of polar glaciers where the formation of superimposed-ice can represent a significant part of the accumulation, such a requirement is not verified and the ELA method cannot be applied. It is also the case for equatorial glaciers where accumulation at the glacier surface is almost constant all year long [e.g., 56]

or glaciers located in mountain ranges where a monsoon regime co-exists with a winter precipitation regime.

The albedo method initially relies on a linear relationship between AAR and SMB, which has been shown to be a common case [37]. It leverages a well defined seasonal albedo-cycle whereby the minimum albedo is reached at the end of the ablation season. Glaciers with less seasonal contrast in surface albedo, due for example to a mixed ablation/accumulation regime such as during monsoon, may complicate this relationship. Nevertheless, Brun et al.

[40] found that minimum glacier-wide albedo could still capture valuable variance of the SMB on a monsoon influenced glacier (Mera) in the Himalayas.

The snow-map method hypothesizes that the seasonal mean snow altitude around a glacier is a proxy of the glacier seasonal SMB. It works well for temperate glaciers whose winter/summer SMB is mainly controlled by snow accumulation/melting. The method relies on the variability of the seasonal snow cover, easily detected when there is a seasonal partitioning of the snow accumulation and ablation. The snow-map method will therefore be less efficient for equatorial or monsoon-regime glaciers.

Spatial distribution of surface mass balance

The presented methods increase the number of monitored glaciers for their glacier-wide annual to seasonal surface mass balances, but none of the methods allows retrieving the spatial distribution of surface mass balance.

In the case of the albedo method, further work is needed in this area and to link the albedo value (annual minimum or integrated albedo value over the winter or summer season) of each pixel to its annual or seasonal surface mass balance. However, the 250-m spatial resolution of the currently used MODIS images might be a limitation because with such a large pixel size, the albedo value of the pixel could hardly be related to the SMB at one point (i.e., a stake) within this pixel. In the case of relatively small mountain glaciers, individual pixels remain potentially affected by varying pixel footprints and spectral mixtures that are compounded by the uncertainties inherent to the topographic and atmospheric corrections and albedo retrieval algorithm.

However, the Sentinel-2 images, available at a 10-m spatial resolution every

five days, open interesting perspectives to resolve albedo cycles from smaller distinct areas of a glacier from which local SMB information may be obtained.

Future challenges

Automating the data processing

Regarding the ELA method, the next step will be automating the method in order to apply it globally or at least for glaciers where the end-of summer snow line can be considered as a proxy of the ELA. Such a global application still needs to tackle some problems such as automating the detection of the endof-summer snow line and the computation of its elevation. Different spectral band combinations or band ratios can be tested [21], associated with cloud and shadows masks. Another important issue for the long-term application of the ELA method is to find DEMs of sufficient accuracy to compute the geodetic balance from the 1970s or 1980s, which for some mountain ranges might not be easy.

Regarding the albedo method, by relying on the systematic processing of frequent image acquisitions, the method requires automating. This has been achieved via the MODImLab MATLAB toolbox whereby the downloading and processing of thousands of MODIS images can be efficiently batched.

MODImLab also provides algorithms to compute isotropic and anisotropic

Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 91 albedos [34] for every MODIS pixels, allowing end-users to easily obtain albedo maps.

Regarding the snow-map method, the data processing is partly automated.

Indeed, a chain of MATLAB routines allow processing the downloaded SPOT/VGT images and the DEM in order to estimate the seasonal SMB.

Use of other images

Regarding the ELA method, the application presented here used highresolution satellite images (i.e. SPOT, Landsat, ASTER). However, as mentioned earlier the method of mapping the end-of-summer snow line can be based on different sources, from very high resolution aerial photographs PROBA-V also provides images at 300 m and 100 m spatial resolution. A study is underway to assess the impact of the PROBA-V image resolution on the accuracy of the seasonal SMB estimated with the snow-map method. In fact, the 100 m resolution PROBA-V images could improve the seasonal SMB estimates, specially the summer ones. The summer mean regional snow altitude detected with the PROBA-V 100 m images could be closer to the snow line notion. The snow-map method applied for the summer period could certainly also take advantage of higher resolution images (like Landsat-8 and/or Sentinel-2) that would enable the inter-annual monitoring of residual small snow-covered surfaces (snow patches). More research is needed in this area.

Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 93

Generalization of the calibrated methods

The methods based on the albedo and regional seasonal snow maps are currently used on glaciers where in situ SMB data are available for calibration.

However, in both cases, a generalization can be made to apply these methods on glaciers where no SMB data are available, as in the case of the ELA method.

Section 2.2.2 proposed a way forward to generalize the albedo method to avoid reliance on current calibration with in situ data. It is suggested that such an approach could mimic the ELA method by using long term geodetic estimates of average glacier wide SMB and the use of the mass-balance gradient. A SMB-area gradient can be made analog to the SMB-elevation gradient via the use of glacier morphometry. The regionalization of such a gradient could then allow estimates of SMB for unmonitored glaciers.

Regarding the snow-map method, a generic relationship showed good performance for winter SMB, even on non-monitored glaciers (without calibrated relationships established between Z and SMB). Indeed, for a dataset of 15 glaciers of the Alps, the winter generic relationship showed a satisfying RMSE (490 mm w.e. yr -1 ). The winter generic model can thus be applied at the scale of the European Alps to sufficiently estimate winter SMB of glaciers for which no in situ data are available. The application of this generic model to other mountain ranges requires further developments.

CONCLUSIONS

We presented a review on three existing methodologies based on optical satellite images to quantify the glacier-wide seasonal to annual surface mass balances. These methods show promising results to considerably increase the number of time series of surface mass balance; a glacier variable necessary to better understand the relationship between climate and glacier changes as well as the contribution of glaciers to the hydrological regime of glacierized catchments. It is worth noting that the best method will vary by region, depending on, for example, cloudiness, glacier size or maximum elevation.

Even if the albedo method and the snow-map method have currently been applied on glaciers where in situ SMB data were available for calibration purposes, both methods can be generalized to be applied on unmonitored glaciers as in the case of the ELA method.

Regarding the seasonal SMB, it is noteworthy that the albedo method performs better for the summer SMB which can be efficiently quantified from the minimum average albedo value of the glacier surface area, or the integration of this value over the entire summer period. On the other hand, the snow-map method shows better performance for the quantification of the winter SMB because the method is based on the variations of the seasonal altitude of the snow cover distribution which are significantly representative of the winter SMB; whereas in summer the variations of the seasonal altitude of the snow cover distribution are lower and limited to representing the high year-to-year variations of summer SMB.

As a consequence, a further exercise could consist of applying the three methods simultaneously to the same glaciers to take advantages of the better performance of each method. In that sense, the winter SMB could be quantified from the winter snow maps, the summer SMB could be quantified from the albedo method and then, a cross validation between the annual SMB quantified from the ELA method and from the sum of the two seasonal components could be done.

Finally, it is worth noting that the three methods were developed, calibrated and validated on temperate glaciers of the mid-latitudes (and outer tropics)

with no or limited debris-coverage. We suggest that for polar glaciers or glaciers in the equatorial and monsoon-regime regions these methodologies might not be efficient and the methods need to be adapted or alternative methods need to be found. This study is divided in three sections; [START_REF][END_REF] the description of the required data;

(2) the principle of snowline detection and the algorithm step-by-step description;

(3) the validation of the semi-automatically derived SLA on 44 glaciers located in the western European Alps.

DATA

The algorithm presented in the current study requires four types of input data: satellite images, a digital elevation model, the glaciers' outline and central flow line which are described in the following sub-sections.

A semi-automated method to estimate the glacier end-of-summer snowline altitude from optical satellite images 101

Satellite images

As the snowline is the limit between two areas of the glacier with contrasted spectral characteristics (snow-covered and bare ice) that can be visually interpreted, optical satellite images have been widely used to map it and track its elevation (Mernild et al., 2013;Pelto, 2011;Rabatel et al., 2005Rabatel et al., , 2008aRabatel et al., , 2013a;;[START_REF] Winsvold | Regional Glacier Mapping Using Optical Satellite Data Time Series[END_REF]. Because of the typical size of mountain glaciers A semi-automated method to estimate the glacier end-of-summer snowline altitude from optical satellite images 103

METHOD

The flowchart presented in Fig. S2 shows the processing and post-processing steps to retrieve the SLA. Each section of the flowchart is developed in this section. To facilitate understanding, some of the steps are illustrated Fig. S3.

The algorithm is applied and validated on the western European Alps where SLAs were manually delineated on 44 glaciers by Rabatel andothers (2013, 2016) but the algorithm is designed to be applied to other glacierized regions.

Principles of the detection

The detection method is based on the bulk processing of all the available images from the previously mentioned sensors, and can be extended to other optical sensors. For each image, the algorithm detects the snowline on the glacier as the steepest gradient in surface spectral signature corresponding to the local transition between ice and snow. This transition, visible on glaciers during the ablation season, is named transient snowline (TSL and TSLA for its altitude). At the end of the ablation season, as the snowline migrates to its highest elevation, it corresponds to the end-of-summer snowline. In its current version, the algorithm is not able to map any snowline contours confidently but rather to capture a range of elevation in which the snowline is. Rectangles represent the required datasets except for the last one that correspond to the product of the algorithm ("One SLA per year and per glacier"). Diamonds represent the processing and post-processing steps computed automatically and the A semi-automated method to estimate the glacier end-of-summer snowline altitude from optical satellite images 105 circle, the manual step of the algorithm. Steps coloured with light blue and numbered with letters (e.g. A,B…) are illustrated Fig. 3.

Processing of satellite images

This section presents the eight processing steps of the algorithm, required to estimate the glacier TSLA and SLA.

Fig. S3 Illustration of some of processing steps presented in the flowchart, Fig. 2. Letters correspond to particular steps from the flowchart (Fig. 2). This example is focused on Glacier Blanc located in the Ecrins mountain range, France. The background image is from Landsat-8 and was taken the 6 th September 2015.

Removal of clouds

In order to avoid erroneous detections induced by cloud coverage on the glacier, a cloud mask is produced at the image scale and is clipped and resampled at the glacier scale (illustrated Fig. S3A). Cloud masking takes advantage of the contrasting spectral signatures of snow and ice surfaces and optically thick clouds (Fig. S1). A Normalized Difference Snow Index (NDSI, Dozier, 1989) is computed to detect the presence of clouds. After manual visualization on our set of studied glaciers, pixels showing a NDSI smaller than 0.25 are considered to be cloud-covered. The chosen threshold is more conservative than usual cloud cover NDSI threshold (i.e. 0.57 in [START_REF] Choi | Cloud detection in Landsat imagery of ice sheets using shadow matching technique and automatic normalized difference snow index threshold value decision[END_REF] to not exclude any pixels of relatively debris-rich ice.

Completely cloud covered images at glacier scale are removed at this stage.

Cosine topographic correction

Optical satellite radiance in the reflective domain is modulated by the geometry of the observed surface [START_REF] Richter | Correction of satellite imagery over mountainous terrain[END_REF], relative to the sun direction (e.g. slope, aspect, …). Over mountainous regions, steep slopes can introduce local variability in observed radiances. This variability can be reduced by correcting the observed radiance using the cosine topographic correction model.

Equation 1 provides the local illumination angle IL of an observed pixel and is obtained using the Euler's rotation theorem. It is computed with the surface slope and aspect angle θsurf, φsurf, the solar zenith and azimuth, θsun and φsun respectively. θsurf and φsurf are derived from the DEM at the same resolution as the processed images, while θsun and φsun are obtained from the image metadata.

) cos( sin sin cos cos cos

surf sun surf sun surf sun IL          (S1)
The illumination angle (cos IL) is used to correct the surface observed radiance to represent a horizontal pixel under a Lambertian assumption (also-called cosine correction) as shown in Equation S2.

IL R R sun hor cos cos  (S2)
R stands for radiance of the input satellite pixel while Rhor stands for radiance of the pixel if the surface was horizontal. To prevent aberrant corrections for A semi-automated method to estimate the glacier end-of-summer snowline altitude from optical satellite images 107 low incident angles (cos IL ~ 0), known as a weakness of this correction model (e.g. Meyer and others, 1993), all pixels with cos IL < 0.2 have been discarded

from the images to avoid artefacts in corrected surface reflectance [START_REF] Davaze | Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data[END_REF]Sirguey, 2009).

Band ratios

Already in the 1970s, the idea of combining spectral bands (e.g. NIR/SWIR) to infer surface properties from optical remote sensing images was introduced (Valovcin, 1976). Several studies have therefore used band ratios to derive surface properties at regional scales, such as global snow cover from NDSI [START_REF] Hall | Development of methods for mapping global snow cover using moderate resolution imaging spectroradiometer data[END_REF], or use of vegetation indices [START_REF] Rouse | Monitoring vegetation systems in the Great Plains with ERTS[END_REF] to assess global change in land use [START_REF] Defries | NDVI-Derived Land Cover Classification at a Global Scale[END_REF]. To identify ice/snow transition, Winsvold and others (2016) recommended the use of the Red/SWIR ratio as the most discriminating. By looking at the different spectral reflectance of snow, ice and clouds (displayed Fig. S1), testing manually several band ratios and using results from Rabatel and others, (2012), we concluded that an enhanced ratio of the one proposed by Valovcin (NIR²/SWIR) was the most suitable to discriminate ice, snow, clouds and extra-glacial features (illustrated Fig. S3B and Fig. S4). Furthermore, it has the advantage to allow topographic corrections on input bands, while using a simple ratio would simplify this correction (refers to Equation S2).

Clipping and resampling

Topographically corrected ratios are resampled to match the DEM resolution (i.e. 30 m) and clipped at the glacier contour, manually modified if necessary, from the RGI 6.0 (as displayed Fig. S3C).

Removal of shadows

Shaded areas can lead to errors in SLA estimation. To mitigate this problem, self and cast shadows corresponding to the time of image acquisition are computed using the DEM. Cast shadows are detected by estimating the sky Performing this classification improves the snowline detection by mitigating the risk of detecting small features as SLA (e.g. snow patch, image artefacts…) due to a dilation-erosion processing. An example is provided Fig. S3F and Fig. S3G).

Estimating TSLAs and deriving the SLA

The glacier is divided in elevation bands of 100 m using the DEM (as displayed Fig. S3F). For each valid image (i.e. manually checked), the classified mode of each band is computed. Subsequently, the algorithm targets the band where the mode changes from the classified value of ice (i.e. 1) to snow (i.e. 2). This step is achieved by computing the first derivative of mode values per altitude band. The first derivative is maximized at the altitude where the mode changes from snow to ice cover. This altitude is subsequently considered as the altitude of the TSL.

Once the elevation band of the TSL is identified, the process is repeated by segmenting three 100 m bands around the first guessed band, into 20 m bands to refine the result (illustration provided Fig. S3G). This yields a single TSLA being identified per glacier and per image, and its altitude is defined as the mid-elevation of the 20 m elevation band where the first derivative is maximal.

To retrieve the annual SLA, the highest altitude reached by the TSL at the end of the ablation season is selected (between 15 July and 20 October in the European Alps).

Identification and filtering of outliers

To mitigate the impact of erroneous detection, we performed a filtering to remove SLA outliers. As the algorithm has been developed to estimate the SLA of unmonitored glaciers (where no data are available for validation), the statistical approach used to detect and remove outliers has to be independent from validation data. The chosen approach is based on a comparison of the retrieved SLA of each individual glacier with SLAs from surrounding glaciers.

The underlying hypothesis is that a common regional response of glaciers to To detect and remove outliers, the SLA of each individual glacier is compared to the regional median SLA of all monitored glaciers from the surroundings.

According to Vincent and others (2017) or [START_REF] Pelto | How Unusual Was 2015 in the 1984-2015 Period of the North Cascade Glacier Annual Mass Balance?[END_REF], the common variance of mountain glaciers SMB is a function of their relative distances. For example, using results from Vincent and others (2017), for two glaciers 100 km apart, 97.5% of the SMB variance is common and of climatic origin. As the SLA is a good proxy of the SMB, we hypothesize that the same relationship applies. To mitigate the impact of morpho-topographic features which could impact the SLA (slope, aspect…), SLA time series are standardized as:

SLAstandardized = SLA 𝑡,𝑔 -Median(SLA 𝑔 ) NMAD 𝑔 (S4)
Where SLA 𝑡,𝑔 stands for SLA at the year t and for the glacier g and Median(SLA 𝑔 ) for the median of all the SLAs of glacier g. NMAD represents the normally distributed median absolute deviation.

We also define a regional SLA signal, corresponding to the median of all the glacier standardized SLAs at year t except for the glacier g in order to be independent of the considered glacier (leave-one-out cross-validation). SLAs at year t and for glacier g are discarded if they fulfill at least one of the two following conditions:

 The standardized SLA is considered as an extreme value regarding the distribution of all the computed standardized SLA (probability density function equal to zero in between the value and the rest of the distribution, SLA at more than 3  from the median).

 For each glacier, a comparison between its standardized SLA time series and the regional SLA signal is performed, along with the 95% confidence interval that takes into account both the SLA uncertainty and the dispersion of the regional SLA signal. Then, a regression is computed following a robust linear model with support for the M-estimators [START_REF] Susanti | [END_REF]. Individual SLAs are removed if lying out of the 95% confidence interval (i.e. too different from the regional SLA signal).

Gap-filling of the SLA timeseries

Removed SLAs are reconstructed using the Lliboutry linear approach [START_REF] Lliboutry | Multivariate Statistical Analysis of Glacier Annual Balances[END_REF]) also applied to SLA by Rabatel and others (2013). This approach (described in Equation S5) takes into account inter-annual (regional effect, according to glaciers from the same region) and spatial (site effect, according to the average of each individual glacier) variabilities to reconstruct missing SLAs:

gt t g gt SLA       (S5)
Where gt SLA stands for the snowline altitude for year t and glacier g, g  is a variable corresponding to the average SLA for glacier g (function of its climatic and topographic context). t

 is a variable specific to the year t, related to meteorological conditions which are represented by the average of all the glaciers for this year, and gt  is the linear model residual, i.e. a random variable with a zero-mean Gaussian distribution [START_REF] Lliboutry | Multivariate Statistical Analysis of Glacier Annual Balances[END_REF]. In other words, the method relies on numerous observations of SLAs providing information on the spatial and temporal variabilities of the SLA. For a long enough period of monitoring (e.g. >10 yr), the behavior of a glacier compared to the regional signal is expected to be conserved every year. Reconstruction of SLA has been only performed when a sufficient number of non-filtered SLAs were available per glacier. A threshold of 31.8% of available SLAs on the whole time series per glacier has been chosen. We consider that if more than 68.2% (±  1 ) of the SLAs were missing, reconstruction of missing SLAs was not robust enough to be performed.

A semi-automated method to estimate the glacier end-of-summer snowline altitude from optical satellite images 113 5.3 RESULTS AND VALIDATION

Errors and uncertainty of the retrieved SLAs

The overall error in the retrieved SLAs is defined as the discrepancy between manually and semi-automatically derived SLAs. For seven glaciers (Blanc, Casset, Saint-Sorlin, Arcelin, Sassière, Mer de Glace, Pré de Bar), SLAs have been manually derived using the same detection footprint and the same images as the algorithm. These glaciers have been chosen because of their various morpho-topographic features (latitude, orientation, slope, altituderange, area), representing the variability of our glaciers dataset and assessing the different sources of errors related to the processing of the images, detection method and the DEM used for the retrieval.

Errors in retrieved SLAs compared to manually derived SLAs can arise from two sources:

 The different detection footprint between manually derived SLAs (detected on the whole glacier width) and semi-automatically derived SLAs (detected on a buffer on the central part of the glacier).  Detection of a steeper gradient than the transition snow/ice due to local topographical context (e.g., snow patches, changing slope, crevasses)

In order to quantify the impact of these sources, we first assessed the impact of detecting the SLA on a buffer along the central part of the glacier. As illustrated Fig. S5A, SLAs have been manually derived within a buffer and on the total glacier width to study the impact of different detection footprint.

Widths of delineated buffers have been calculated using Eq. S3. On the 118 compared SLAs, no significant mean bias can be observed ( 1 

Bias

). This observation implies that the detection of the snowline on a central buffer does not tend to overestimate or underestimate its altitude, at least on the seven studied glaciers. There is a statistically significant relationship at a confidence

Validation of the retrieved SLA

After the validation on seven glaciers, the general performance of the algorithm has been assessed using the set of annual SLA from 44 glaciers of the western European Alps derived manually by Rabatel and others (2013).

Figure S6A compares manual and semi-automatic SLA estimations on 44 glaciers over the 18 study years, without filtering and gap-filling the derived SLAs. For this comparison, images used for the manual delineation and by the algorithm can be different. One can note the presence of outliers, with discrepancies of the same order than the comparison previously made on seven glaciers. The observed mean bias of +36 m, close to the +44 m computed previously, is small and not statistically significant compared to the estimated errors in retrieved SLA (𝛥SLA=±123 m), with causes of this bias being the same than enounced previously. There is a significant correlation at the 1% risk level between manually and semi-automatically retrieved SLA, with 35% of common variance. The RMSE of 158 m is of the same order of magnitude than the uncertainty derived previously ( 𝛥SLA=123m), even if slightly higher. Comparisons have been made at glacier scale and r², RMSE and mean biases have been computed in Table 2 for each of the 44 glaciers.

For the raw values of retrieved SLAs (i.e. before the filtering step), we observed an average r² of 0.3 (ranging from 0 to 0.75), an average RMSE of 140 m (ranging from 48 m to 296 m) and a mean bias of about 40 m. These statistics reveals that the algorithm is able to derive SLAs and to capture the inter-annual variability for a large share of the 44 glaciers (refers to Table S2 for list of involved glaciers). However, there are for some glaciers (e.g., Rouies, Mahure, Des Glaciers…) and for some years, outliers leading to poor correlation coefficients at glacier scale.

A semi-automated method to estimate the glacier end-of-summer snowline altitude from optical satellite images 117 As these outliers do not reflect real fluctuations of SLAs but are artefacts of the algorithm, the developed filtering and gap-filling method has been applied. Filtering and gap-filled SLA are then compared to the same manually derived SLAs. The comparison has been made on 34 glaciers, as the filtering step led to the removal of 10 glaciers on which the method is not applicable because of peculiar morpho-topographic conditions making the detection very challenging (e.g. presence of snow-patches, narrow, steep or highly crevassed glaciers with important seracs falls). As illustrated Fig. S6B, filtering and gapfilling led to the removal of a large share of observed outliers in Fig. S6A and improved the correlation (r 2 =0.41 and RMSE = 144 m). However, it did not solve the mean observed bias that remains almost equal (from +36 to +35 m), as both overestimated and underestimated SLAs have been filtered. At glacier scale (Table S2), the filtering and gap-filling steps improved the correlation with the reference dataset, rising the average r² from 0.30 to 0.34.

Lowering of the average RMSE and mean bias is also observed, from 140 to 97 m and from +40 to +33 m respectively, which indicates that a large share of outliers due to detection errors have been removed at this stage. As displayed Fig. S6B, the filtering step is not able to remove some of the outliers, which led to the observed positive bias.

Finally, on the initial set of semi-automatically retrieved SLAs corresponding to 725 values, considering 44 glaciers over 18 years, and 67 missing SLAs due to cloud coverage, 284 have been filtered and 104 reconstructed using the Lliboutry approach. 180 SLAs were not reconstructed because of the 10 glaciers discarded by the filtering process. In addition to allow gap-free SLA time series to be reconstructed on mountain glaciers, filtering and gap-filling retrieved SLA significantly improved the agreement with manually delineated SLAs.

Table S2. r², RMSE and mean bias of the comparison between manually and automatically derived SLAs. Glaciers are listed according to the mountain range they belong. Coefficients have been computed for retrieved SLAs and gap-filled SLAs. For glaciers removed in the filtering steps, "-" indicates that no coefficients have been computed as the glacier has been discarded. On an indicative basis, for an 18-yr period (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016), the correlation is significant at the 1% risk level for r²>0.34, and for r²>0.21 at the 5% risk level. For initial and filtered-gap-filled time series, possible missing or filtered SLAs changes the number of samples of the comparison. This also changes the threshold of statistical significance, increasing when the number of samples decreases. 120 5.3.3 Temporal variability of the retrieved regional SLA

To assess the ability of the algorithm to capture the SLA regional variability, we computed for both manual and semi-automatic retrieval, the annual median and distribution of SLAs for the 34 glaciers remaining from the previous steps. The temporal comparison between both datasets is illustrated These conditions make the automatic detection more challenging and can explain the general overestimation of the SLA during this year. Same conclusions can be made for 2003, which was an exceptional warm year with very high SLAs recorded in the European Alps [START_REF] Paul | On the impact of glacier albedo under conditions of extreme glacier melt: the summer of 2003 in the Alps[END_REF]Rabatel et al., 2013a). In 2003, the SLA was close to the glacier top for glaciers with a rather low maximum altitude (~3,300 m a.s.l.), for which the algorithm fails to derive the SLA out of the detection footprint (100 m below glacier top). This inability led the semi-automatic retrieved SLAs to be underestimated for such exceptional warm years, caused by estimation of late TSLAs rather than SLAs.

Adapting the choice of glaciers considering their maximum elevation would solve this problem but could also reduce the number of studied glaciers. On the other hand, the variability is more pronounced for manually delineated SLAs than for semi-automatically derived SLAs as illustrated by the boxplots in Fig. S8A and S8B (representing manual and semi-automatic 122 derived SLAs respectively). It is noteworthy that the manual method delineates the SLA on the whole glacier width, while the algorithm estimates a range of altitude where the SLA lies. Manual delineation is therefore able to monitor a greater variability in monitored SLA, especially for exceptionally high and low SLA (i.e. 2001SLA (i.e. , 2003SLA (i.e. and 2009)). Likewise, years with exceptionally high or low SLA are not well captured by the algorithm because of the previously enounced reasons.

Small discrepancies between manually and semi-automatically derived

Fig. S8 Annual variability of manually and semi-automatically retrieved SLAs (Fig. 8A and 8B respectively). For each time series, boxplots represent, from top to bottom; maximum, quartile 75%, median, quartile 25% and minimum SLA for a given year.

A semi-automated method to estimate the glacier end-of-summer snowline altitude from optical satellite images 123

Lines display the temporal variability of the annual median SLA, also illustrated Fig. 7.

Aware of these limitations, the proposed approach allowed semiautomatically deriving continuous SLA time series for 34 mountain glaciers over an 18-yr period.

LIMITATIONS OF THE METHODS

Despite the overall performances of the proposed algorithm, some limitations may affect its success. First, if no snowline can be identified for a given image (undetected clouds or complete snow cover), the algorithm will still derive the steepest gradient, leading to an erroneous estimate of a TSLA. This is currently mitigated by the visual check prior to the computation of the SLA. Such source availability [START_REF] Lamare | Optical satellite retrievals of snow physical properties in complex terrain[END_REF]. We hence encourage future endusers of our algorithm to consider this limitation on glaciers prone to inaccurate topographic correction (e.g., very steep and crevassed glaciers, glacier subject to large DEM errors).

Finally, our method differs from the recently published method from [START_REF] Rastner | On the Automated Mapping of Snow Cover on Glaciers and Calculation of Snow Line Altitudes from Multi-Temporal Landsat Data[END_REF], which is based on the differentiation of ice and snow by applying a dynamic threshold on the reflectance of Landsat NIR band on each considered image. The method from [START_REF] Rastner | On the Automated Mapping of Snow Cover on Glaciers and Calculation of Snow Line Altitudes from Multi-Temporal Landsat Data[END_REF], even if able to derive a snow cover area on the glacier and to spatialize the derived SLA on 26 glaciers from glaciers in the Ötztaler Alps, has never been applied at regional scale.

INTRODUCTION

Beyond their iconic role of climate change [START_REF] Mackintosh | Reconstructing Climate from Glaciers[END_REF], mountain glaciers are considered as natural "climate-meter" [START_REF] Vaughan | Observations: Cryosphere[END_REF] because of their high sensitivity to climate change. Present in ~26% of the main hydrological basins, about one-third of the Earth inhabitants directly or indirectly rely on water from glacier melt [START_REF] Beniston | Climatic Change in Mountain Regions: A Review of Possible Impacts[END_REF]. Because of the rapid warming, amplified in mountainous regions [START_REF] Research | Elevationdependent warming in mountain regions of the world[END_REF], mountain glaciers are losing mass worldwide [START_REF] Zemp | Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016[END_REF], contributing to a third of the observed sea-level rise (Church et al., 2013), and are expected to reach their peak of ice mass loss from 2040 to 2100, function of the considered climate scenario (Marzeion et al., 2012). The glacier-wide surface mass balance (SMB) is considered as the most comprehensive metric to assess both the glacier mass and dynamics changes, and interactions with local climate. Quantifying the SMB of individual glaciers at regional scale remains today challenging, especially using the "logistically expensive" glaciological method from in situ data. Recent methods have taken advantage of the multiplicity of remote-sensing sensors, leading to regional estimates of glacier multi-decadal average mass balance, using gravimetry (e.g. Jacob et al., 2012), geodesy (e.g. [START_REF] Brun | A spatially resolved estimate of High Mountain Asia glacier mass balances, 2000-2016[END_REF][START_REF] Menounos | Heterogeneous Changes in Western North American Glaciers Linked to Decadal Variability in Zonal Wind Strength[END_REF], altimetry (e.g. [START_REF] Bolch | Mass loss of Greenland's glaciers and ice caps 2003-2008 revealed from ICESat laser altimetry data[END_REF]Gardner et al., 2013) or glacier length change [START_REF] Hoelzle | Secular glacier mass balances derived from cumulative glacier length changes[END_REF]. However, these methods fail in estimating the inter-annual variability of the SMB, providing critical insights of the impact of climate change at a regional scale, and valuable dataset for hydrological studies in glacierized catchments, or for model parameterization.

In this study, we propose an automation and a regional application of the SLA-method (Braithwaite, 1984;Rabatel et al., 2005[START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF]Rabatel et al., , 2017)), based on the estimation of the end-of-summer snowline altitude (SLA), known as a proxy of the equilibrium-line altitude (ELA, [START_REF] Lliboutry | Traité de glaciologie[END_REF], and used to compute the annual SMB time series. Obstacles of a regional application of this method were relying on the automatic delineation of the SLA and the availability of geodetic mass balance estimates used to quantify the multidecadal average mass balance needed to retrieve the annual SMB in the SLAmethod (Rabatel et al., 2017). The recent development of algorithms to automate the processing of satellite images (e.g. [START_REF] Beyer | The Ames Stereo Pipeline: NASA's Open Source Software for Deriving and Processing Terrain Data[END_REF][START_REF] Noh | Automated stereophotogrammetric DEM generation at high latitudes: Surface Extraction with TIN-based Search-space Minimization (SETSM) validation and demonstration over glaciated regions[END_REF]Shean et al., 2016;Sirguey, 2009), the increasing storage and computing capacities and the multiplicity of optical remote sensing sensors, suitable to monitor glacial features (e.g. Landsat, Sentinel-2, ASTER) have

Region-wide annual glacier surface mass balance for the European Alps from 2000 to 2016 129 encouraged the development of automatic methods (Berthier et al., 2016;[START_REF] Brun | A spatially resolved estimate of High Mountain Asia glacier mass balances, 2000-2016[END_REF][START_REF] Rastner | On the Automated Mapping of Snow Cover on Glaciers and Calculation of Snow Line Altitudes from Multi-Temporal Landsat Data[END_REF], enabling to overcome these obstacles.

Despite their minor potential impact on sea-level rise [START_REF] Huss | Distributed ice thickness and volume of all glaciers around the globe[END_REF][START_REF] Zemp | Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016[END_REF] and the relatively low ice extent (2,092 km2, RGI 6.0, RGI [START_REF] Consortium | [END_REF], we applied the SLA-method on the European Alps for the period 2000-2016 because of: i) the availability of a high number of in situ monitored glaciers with SMB data (23 in this study), among which nine are considered as reference glaciers by the World Glacier Monitoring Service (WGMS, see Table S3 for detailed list); and ii) manually delineated snowline from 44 glaciers (Rabatel et al., 2013a[START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF] used to validate our estimates.

Besides providing a new estimate of the mean European Alps mass balance for the period 2000-2016, we used these data to: i) propose and assess a methodological framework to derive the annual SMB times series for individual glaciers at regional scale from optical remote sensing data; ii) investigate the spatial and temporal patterns of SMB fluctuations across the European Alps;

iii) assess the possible climatic and morpho-topographic drivers affecting the SMB variability at sub-regional scale; and iv) determine the representativeness of the available annual SMB time series quantified from in situ data and gathered by the WGMS.

STUDY AREA

We define the European Alps as the region located between the coordinates 44° to 48°N and 5° to 14°E. The European Alps are the second most glacierized region of Europe after Scandinavia (2,092 vs. 2,949 km2 according to the Randolph Glacier Inventory 6.0, RGI Consortium, 2017). Two sets of glaciers within the European Alps have been chosen in the present study. The first set corresponds to glaciers selected to compute the geodetic mass balance using the ASTER images archive (see section 1.3.1.); larger than 0.5 km2 as recommended by [START_REF] Menounos | Heterogeneous Changes in Western North American Glaciers Linked to Decadal Variability in Zonal Wind Strength[END_REF] who applied the same method to quantify the geodetic mass balance. The second set corresponds to glaciers selected to compute the annual SMB time series using the SLA-method.

In order to be able to detect the SLA for each year of the studied period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016), glaciers of the second set were selected using the following criteria: [START_REF][END_REF] glaciers had to have a high enough maximum elevation to allow observation of the snowline every year; (2) glaciers had to be larger than 0.5 km2 because of the dependency to geodetic mass balance to compute the annual SMB. The total considered glacier number and glacierized area by dataset are presented Table 1. Glaciers for which the annual SMB time series have been computed

are distributed over 10 sub-regions, on a 1x1 degree grid. In order to assess the glacier mass change over the 17 years of the study period, we computed both the geodetic mass balance corresponding to a multidecadal average and the annual SMB time series from the SLA-method at the glacier scale. Berthier et al., 2016;Gardelle et al., 2013) and allowed to refine their contribution to sea-level rise (Gardner et al., 2013;[START_REF] Zemp | Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016[END_REF].

Following the recent studies using the large ASTER archive (e.g. Berthier [START_REF] Consortium | [END_REF].

In order to estimate the decadally-averaged glacier-wide mass balance, surface elevation changes are converted into mass changes applying the consensual volume-to-mass conversion of 850±60 kg.m-3 (Huss, 2013).

Reconstructing annual surface mass balance time series from optical satellite images

The geodetic approach provides spatially resolved mass balance estimation over decadal time spans, but fails in estimating the annual fluctuations for mountain glaciers. In this study, we used the SLA-method, based on the identification of end-of-summer snowline altitude (SLA), representative of the equilibrium-line altitude [START_REF] Lliboutry | Traité de glaciologie[END_REF] which can in turn be used as an indicator of the annual SMB (Braithwaite, 1984;[START_REF] Kuhn | The Response of the Equilibrium Line Altitude to Climate Fluctuations: Theory and Observations, Glacier Fluct[END_REF]. The SLA-method has been proposed by Braithwaite in 1984, and later applied by, among others, [START_REF] Dedieu | Apport des données SPOT et Landsat TM pour le suivi de la fusion nivale et des bilans glaciaires dans les Alpes françaises[END_REF], Rabatel et al. (2005Rabatel et al. ( , 2008[START_REF] Rabatel | Can the snowline be used as an indicator of the equilibrium line and mass balance for glaciers in the outer tropics ?[END_REF][START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF] and [START_REF] Chandrasekharan | Quantification of annual glacier surface mass balance for the Chhota Shigri Glacier, Western Himalayas, India using an Equilibrium-Line Altitude (ELA) based approach[END_REF]. The first step consists in estimating the SLA.

Because of the large number of glacier considered in this study (239), we semi-automatically identified the transient snowline altitude (TSLA) using the workflow detailed in the supplementary material (Sect. 5 of the manuscript).

We used multi-sensor optical satellite images of decametric resolution (i.e. LANDSAT 5,7,8,Sentinel2,ASTER) and identified the TSLA using a band ratio (NIR 2 /SWIR), enhancing the snow/ice contrast. The ASTER GDEM v2 has been used in the present study to compute the snowline altitude. For each image, the workflow allows detecting the TSLA on the glacier as the steepest gradient in the surface band ratio signature corresponding to the local transition between ice and snow. For each glacier we then considered the highest annual derived TSLA to be representative of the glacier annual SLA. Missing SLAs, due to cloud cover or image unavailability are gap-filled using a statistical linear approach [START_REF] Lliboutry | Multivariate Statistical Analysis of Glacier Annual Balances[END_REF] taking into account inter-annual (regional effect, according to glaciers from the surrounding region) and spatial (site effect, according to the average of each individual glacier) variability to reconstruct missing SLAs. Following the conclusion by [START_REF] Vincent | Common climatic signal from glaciers in the European Alps over the last 50 years[END_REF] and Pelto (2018), we considered that glaciers 100 km apart have a significant common variance (97.5% in the European Alps for SMB data, [START_REF] Vincent | Common climatic signal from glaciers in the European Alps over the last 50 years[END_REF]. As the SLA is a good proxy of the SMB, we assumed that the same relationship applies.

Following the methodology detailed in the supplementary material, 239 SLA time series have been derived. Using the method presented in Rabatel et al. (2016), the glacier-wide SMB has been reconstructed for the period 2000-2016. In summary, this approach computes for each year the departure between the annual SLA and ELAeq, the theoretical altitude of the equilibrium line if the glacier had a balanced budget for the study period. ELAeq is computed using the mean glacier-wide mass balance over the period, obtained with the geodetic method, and the mass-balance gradient (∂b/∂z) in the vicinity of the ELA. Although peculiar to each glacier, several studies [START_REF] Benn | Mass balance and equilibrium-line altitudes of glaciers in high-mountain environments[END_REF][START_REF] Chandrasekharan | Quantification of annual glacier surface mass balance for the Chhota Shigri Glacier, Western Himalayas, India using an Equilibrium-Line Altitude (ELA) based approach[END_REF][START_REF] Kuhn | Mass Budget Imbalances as Criterion for a Climatic Classification of Glaciers[END_REF][START_REF] Mayo | Glacier Mass Balance and Runoff Research in the U. S. A[END_REF]Rabatel et al., 2005[START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF] demonstrated that choosing a regional ∂b/∂z does not depreciated the results and lead to similar results that individual glacier ∂b/∂z estimates. We therefore selected a mass-balance gradient of 0.0078 m w.e. m-1, according to Rabatel et al. (2005).  the uncertainty related to the mass balance gradient in the vicinity of the ELA, considered as 10% of db/dz as suggested by Rabatel et al. (2005).

 the uncertainty in the computed geodetic mass balance presented above.

 the time difference between the image acquisition used to estimate SLA, and the last day of the hydrological year. [START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF] estimated this uncertainty at ±0.09 m.we. a-1 in average on 30 glaciers in the French Alps for a 30-yr study period. This estimates made for glaciers in the French Alps is considered here at the scale of the European Alps.

Finally, the uncertainty in the computed annual SMB with the SLA-method equals ±0.26 m w.e. a -1 on average, ranging from 0.22 to 0.4 mw.e. a -1 depending on the year and the glacier.

Remotely-sensed SMB validation

Using the SLA-method, the glacier-wide annual SMB time series of 239 glaciers were reconstructed for the period 2000-2016. The European Alps being one of the most monitored glacierized regions worldwide, we could validate our estimates with 23 glacier-wide annual SMB time series estimated from in situ data and computed with the glaciological method (full list presented Table S3 in the supplementary material).

As displayed in Figure 1A, annual SMB from remote sensing and from in situ data are in general agreement, with 33% of common variance, no observed mean bias and RMSE (0.6 m w.e.) of the same order of magnitude than the computed remotely-sensed SMB uncertainty (0.26 m w.e. in average).

There is however a small underestimation for positive and highly negative SMB, resulting in the slope of the linear regression to be smaller than 1 (a = 0.4).

When comparing the standardized glacier-wide annual SMB of the 23 glaciers with in situ estimates (Figure 1B), the interannual variability agrees quite well except in between 2004 and 2008. There are three possible explanations for this discrepancy: i) errors in SLA estimation from the SLA-method; ii) representativeness of the computed glacier-wide SMB from sparse in situ data; and iii) difference of glacier extent used to compute the glacier-wide SMB in the two methods. Indeed, when comparing manually delineated and semiautomatically derived SLA (see supplementary material, Sect. 5 of this manuscript), the algorithm fails in reproducing extremely high or low SLA due to the position of the SLA itself, either too close to the glacier head or snoot.

Therefore, retrieved SMB from the SLA-method would not well capture positive or extremely negative years leading to the observed discrepancy. Furthermore, Region-wide annual glacier surface mass balance for the European Alps from 2000 to 2016 135 most of the SLA estimates between 2003 and 2008 were based on Landsat 7 images, for which the Scan Line Corrector (SLC) failed in May 2003, resulting in stripped images reducing drastically the number of usable images to detect the SLA. On the other hand, glacier-wide estimates from in situ data using the glaciological method is dependent from the point measurements distribution on the glacier (often biased toward the ablation area due to access facilities), the interpolation method, and glacier area change during the monitored period (e.g., [START_REF] Basantes-Serrano | Slight mass loss revealed by reanalyzing glacier mass-balance observations on Glaciar Antisana 15<span class="italic">α</span> (inner tropics) during the 1995-2012 period[END_REF][START_REF] Rabatel | Estimation of Glacier Thickness From Surface Mass Balance and Ice Flow Velocities: A Case Study on Argentière Glacier, France[END_REF][START_REF] Thibert | Glaciological and volumetric mass-balance measurements: error analysis over 51 years for Glacier de Sarennes, French Alps[END_REF].

Geodetic calibrations can reduce the related errors when photogrammetry acquisitions are performed [START_REF] Zemp | Reanalysing glacier mass balance measurement series[END_REF]. In addition, in situ SMB can show surprising values, as for the north-facing Corbassière Glacier in 2006, with a measured annual SMB of -3.56 m w.e. a -1 and a corresponding ELA at 3,875 m a.s.l., which is more than 450 m above the uppermost manually delineated SLA from neighbouring French glaciers for the period 1984-2010 (Rabatel et al., 2013). Such negative SMB values from in situ measurements can be questionable and a source of discrepancy with remote-sensing estimates. Finally, the different glacier surface areas considered in the two methods can explain part of the differences. Indeed, glacier surface areas used for the SLA and geodetic methods correspond to the glacier outlines mostly from 2003 regarding the European Alps (RGI Consortium, 2017) while the year of glacier extent delineation used for the glaciological method is variable from glacier to glacier, and dependent of the local inventories. Different glacier surface areas can therefore have an impact on the retrieved annual SMB values but is challenging to estimate.

Despite the observed discrepancies, computed annual SMBs from the SLA-method are able to reproduce the monitored SMB from in situ measurements together with the interannual variability (Fig. 1B) and allow applying the SLA-method at regional scale.

Figure 1 : Annual SMBs computed from remote-sensing data using the SLA-method versus annual SMBs from in situ data using the glaciological method. In Fig. 1A, the red line illustrates a robust linear regression with support for the M-estimators (Susanti and others, 2014), together with the 95% confidence interval, illustrated with grey lines. Statistics are displayed to assess the correlation between both datasets; a represents the slope of the linear regression and b the intercept. Heatmap style colours are given as an indicator of the point density with one point representing one compared annual SMB. In Fig. 1B, the median standardized interannual variability is represented for in situ and remote-sensing estimates for the 23 glaciers monitored in situ and detailed Table S3. 6.4.3 Glacier morpho-topographic features, one of the multi-decadal average SMB drivers?

Previous studies showed that multi-decadal average mass balances of individual glaciers can be related to glacier morpho-topography (e.g. [START_REF] Brun | Heterogeneous Influence of Glacier Morphology on the Mass Balance Variability in High Mountain Asia[END_REF][START_REF] Fischer | Surface elevation and mass changes of all Swiss glaciers 1980-2010[END_REF][START_REF] Paul | Spatial variability of glacier elevation changes in the Swiss Alps obtained from two digital elevation models[END_REF]Rabatel et al., 2013[START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF]. In addition, according to Huss (2012) and Abermann et al.

(2011) differences in specific mass balance between neighbouring glaciers can often be more controlled by glacier morpho-topography than climate variability, and observed mass balance response to similar climate forcing can differ by up to a factor of four from one glacier to another. However, assessing the impact of glacier morpho-topography on the average glacier mass balance was challenging due to the limited number of in situ measurements, not necessarily representative of a regional glacier diversity A B

Region-wide annual glacier surface mass balance for the European Alps from 2000 to 2016 137

We willingly computed the mean glacier slope rather than the largely used glacier tongue slope to also take into account accumulation basin geometries.

For the six studied features, three can be significantly correlated to computed multi-decadal average mass balance (with a statistical risk of error <5% according to a Student's test); (1) the median altitude (r 2 =0.24), ( 2) the mean glacier slope (r 2 =0.19) and ( 3) the maximum altitude (r 2 =0.12). Using a least absolute shrinkage and selection operator (LASSO) regression analysis, the three variables explain 36% of the observed geodetic multi-decadal average mass balance variance, with the slope and the median altitude being the (Gardner et al., 2013). Therefore, geodetic estimates of multi-decadal mass changes for individual glaciers at regional scale provide a lever to investigate the impact of glacier geometry on the decadal mass balance variability from one glacier to another. We tested the relationship between the geodetic multi-decadal average mass balance for 361 glaciers for the period 2000-2016, and several glaciers morpho-topographic features (i.e. mean glacier slope, minimum, median and maximum altitude, aspect and area, displayed dominant terms. Adding other variables (area, minimum altitude and mean aspect) only rises the explained variance to 39%, which underlines their relative low impact on the multi-decadal average mass balance at regional scale. This result strengthens previous observations made within the European Alps (e.g. [START_REF] Fischer | Surface elevation and mass changes of all Swiss glaciers 1980-2010[END_REF][START_REF] Huss | Distributed ice thickness and volume of all glaciers around the globe[END_REF][START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF] and in other mountain regions (e.g., [START_REF] Brun | Heterogeneous Influence of Glacier Morphology on the Mass Balance Variability in High Mountain Asia[END_REF], even if we did not find a significant statistical correlation between the glacier area and the multi-decadal average mass balance as suggested by previous studies [START_REF] Huss | Distributed ice thickness and volume of all glaciers around the globe[END_REF][START_REF] Paul | Spatial variability of glacier elevation changes in the Swiss Alps obtained from two digital elevation models[END_REF]. This difference could also rest on the absence of small glaciers (<0.5km 2 ) in our dataset.

It is finally worth noting that correlations are positive for the three significantly correlated morpho-topographic variables, suggesting that steep and high glaciers with a high accumulation basin experience the less negative multi-decadal average mass balances.

6.5 DISCUSSION (2000-2009 and 2009-2016) as done in the study by [START_REF] Menounos | Heterogeneous Changes in Western North American Glaciers Linked to Decadal Variability in Zonal Wind Strength[END_REF] and [START_REF] Gurriaran | Annual glacier-wide mass balances of 82 glaciers in the tropical Andes over a 18-yr period from optical remote sensing data[END_REF]. Results for the two sub-periods are not displayed because no contrasted pattern is observed between the two sub-periods, as also illustrated by the cumulative annual SMB (Fig. 6). Regarding the inter-annual SMB variability, similar patterns can be observed for some sub-regions (Fig. 4). To investigate these patterns, we computed the Mahalanobis distance [START_REF] Mahalanobis | On the generalized distance in statistics[END_REF] between each of the sub-regions and represented their respective distances in Figure 5. This distance allows seeing whether spatial patterns emerge and if sub-regions show similar temporal variability in their SMB time series.

mean mass balance over 2000-2016 computed with the geodetic method, with white slices corresponding to the percentage of glacierized pixels not monitored because of large uncertainties. Numbers (1 to 10) denominate regions where annual SMB time series have been computed using the SLA-method. Dashed green and continuous orange polygons illustrate the two main groups of sub-regions identified from the analysis of inter-annual SMB time series.

Figure 4 : Median SMB time series (in m.w.e. a-1) for each of the ten studied subregions (displayed Figure 3). The confidence interval in light grey illustrates the annual dispersion of the SMB estimates for individual glaciers in the sub-region, computed using the normalized median absolute deviation (NMAD). The number of glaciers is displayed on the up-left corner of each subplot (n=). The red histograms sort each annual SMB as a function of its value, in m w.e. a -1 .

Analysis of possible climatic drivers

In order to understand the observed heterogeneity of the computed annual SMB time series, we investigated different climatic variables during summer (May to September included) and winter (October to April included) from the ERA5 reanalysis. For each 1x1 degree grid cell, we computed the summer average of the 2-m air temperature (T2ms), sea-level pressure (SLPs), zonal and meridian temperature advection (Tus and Tvs), and the winter average of precipitations (Pw), sea-level pressure (SLPw), precipitable water over the entire atmospheric column (Pwatw) and the zonal and meridian humidity advection (Quw and Qvw). We then performed, per grid cell, a LASSO regression analysis between these climatic variables and the median annual SMB time series over the 2000-2016 period. For each comparison, a maximum of three variables explaining most of the observed variance are identified and reported Table 2, with red variables being the first, salmon the second and blue the third variable in terms of explained variance. It has to be noted that this statistical analysis is conducted over a 17-yr time period, which from a climatic point of view is rather short. Therefore, the interpretation of the series, while a high distance (i.e. blue) indicates heterogeneous pattern between the compared sub-regions.

Region-wide annual glacier surface mass balance for the European Alps from 2000 to 2016 143 results has to be considered with caution and as a preliminary study of the relationship between the SMB and climate inter-annual variability at the European Alps scale.

Interestingly, from the LASSO analysis, the two groups of sub-regions identified from the SMB time series also emerge. For the sub-regions #1, #5, #6 and #7, the only climatic variable statistically significant and explaining more than 30% of the observed inter-annual SMB variance is the precipitable water during winter, which is largely impacted by the winter air temperature (maximum water content in humid air depending on temperature). This is in agreement with observations made for glaciers in these regions (i.e. Austrian Alps, [START_REF] Abermann | Climatic controls of glacier distribution and glacier changes in Austria[END_REF] suggesting that glaciers in this area could be maintained because of "above-average accumulation".

Sub-regions #3, #4, #9 and #10 have the same first explicative variable (Quw, except for #3 where it is second), corresponding to eastward advection of humidity and seems also highly impacted by summer air temperature and summer meteorological conditions (Tvs, SLPs and T2ms), mainly dominated by southward fluxes. These sub-regions are all in the western European Alps and could therefore experience similar meteorological forcing during zonal or meridian meteorological events.

Sub-region #2 seems to be highly impacted by summer meteorological conditions and could not be assimilated to any other sub-region. In contrary to sub-region #3, #2 seems not impacted by temperature or humidity advection and is too distant from #3 to be grouped with, despite their common high dependence to summer air temperature. Finally, the sub-region #8 is dominated by summer meteorological conditions, winter eastward humidity advection and precipitable water. The predominating drivers and observed inter-annual SMB are not identical from the surrounding regions despite their close geographical distances, and our analysis does not allow to conclude on the reasons for this peculiar behaviour.

Surprisingly, winter precipitation is not identified as an explicative variable for any of the monitored regions. One hypothesis could rest on the poor capacity of ERA5 to spatialize precipitation regarding the relatively low resolution of the used DEM (i.e. 30 km), particularly important in mountainous terrain. Precipitable water is often considered as more accurate as it does not perform a spatial estimation of precipitation but rather a budget of precipitable water.

In situ monitored SMB representativity

Historically, glacier-wide SMB is monitored in situ and glaciers have often been chosen because of their access, morpho-topographic features enabling their instrumentation or political and historical reasons [START_REF] Kaser | A manual for monitoring the mass balance of mountain glaciers[END_REF]. As the "ideal" glacier for mass balance investigation does not represent the global glacier diversity, biases can occur when considering such glacier as representative of the regional scale glacier behaviour. Indeed, negative biases have already been observed when comparing spatially resolved SMB from geodetic estimations to interpolations of sparse glaciological measurements, not representative for large glacierized regions (Gardner et al., 2013).

Recently developed geodetic methods allow overcoming this issue of representativeness due to temporal and spatial resolutions of the computed SMB (e.g. [START_REF] Zemp | Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016[END_REF][START_REF] Menounos | Heterogeneous Changes in Western North American Glaciers Linked to Decadal Variability in Zonal Wind Strength[END_REF]. 6). As shown in Figure 2, the gentler the glacier slope, the higher the overall mass loss. This relationship has already been observed, either sporadically on a small set of glaciers (e.g. Rabatel et al., 2013[START_REF] Huss | Distributed ice thickness and volume of all glaciers around the globe[END_REF], or at a larger scale (e.g. [START_REF] Brun | Heterogeneous Influence of Glacier Morphology on the Mass Balance Variability in High Mountain Asia[END_REF], Fischer et al., 2015). A possible interpretation of this observation is that the glacier mean slope and the glacier response time are related with steeper glaciers having a faster dynamics and adapting quickly to changing climate conditions compared to glaciers with low slopes. This hypothesis would imply that steeper glaciers faster reach a balanced-state with less negative SMB while glaciers with low slope remain in imbalance for a longer time, experiencing more negative SMB to adapt to changing climatic conditions (corroborated by [START_REF] Brun | Heterogeneous Influence of Glacier Morphology on the Mass Balance Variability in High Mountain Asia[END_REF].

Another hypothesis would rest on an increase in precipitation over the period 2000-2016 that would benefit steeper glacier because of their faster mass turnover. The studying period is however too short to conclude on a significant increase in precipitation (illustrated Fig. S9, and this hypothesis is rather speculative and would require a longer observation period and assessment of Region-wide annual glacier surface mass balance for the European Alps from 2000 to 2016 147 the glacier dynamics response to increasing precipitation with numerical models.

On the other hand, the in situ SMB time series fails in representing glaciers from class 3 and 4 (slopes greater than 20.5°), while these classes account for 45% of the total European Alps glacierized area. We therefore suggest that in situ SMB time series may illustrate a negatively biased picture of the regional behaviour because of their poor representativeness for the steeper glaciers representing almost half of the glacier coverage.

A similar analysis has been conducted for the glacier median altitude but the only class standing out from the crowd (with median altitude greater than 3,207 m a.s.l. versus 3,103 m a.s.l. in average for in situ monitored glaciers) concerns about 20% of the total European Alps glacierized area and glacier from the glaciological sample are spanning over the four defined classes, not permitting to conclude on significant relationship or mis-representativeness of the glaciological sample. For the maximum altitude, when divided our glacier subset into four classes as a function of their maximum altitude, glaciers with higher maximum elevations seem to experience less negative SMB. It is important to highlight that SLAs are rarely observable for the entire studied period on glaciers with low maximum elevations (<3,400 m a.s.l.). Low elevation glaciers have not been selected because SLAs were not identifiable every year on these glaciers. Therefore, our subset is not representative of the total European Alps population (12% in our subset vs. 46% for the entire European Alps, in terms of corresponding area for glaciers with maximum altitude <3,471 m a.s.l.).

We therefore suggest that because of a low representativeness for steep glaciers, in particular with slopes greater than 20°, studies estimating regionscale SMB from in situ extrapolation could be substantially biased, if morphotopographic features are not considered. This could be a source of differences of our regional mass balance estimates (-0.74±0.20 m w.e. a -1 ), ~25% lower than the estimate from [START_REF] Huss | Distributed ice thickness and volume of all glaciers around the globe[END_REF] and ~13% lower than the one from [START_REF] Zemp | Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016[END_REF]. We hypothesize that a part of the observed difference could be attributed to our inclusion of steep glaciers, and our spatially resolved estimated.

152

INTRODUCTION

Mountain glaciers represent only 3% of the ice volume on the Earth but contribute significantly to sea level rise (e.g. Church et al., 2013;Gardner et al., 2013;Jacob et al., 2012). In addition, millions of people partly rely on glaciers, either for drinking water or agriculture or due to related glacier hazards (Baraer et al., 2012;[START_REF] Chen | ESTIMATION OF ALPINE GLACIER WATER RESOURCES AND THEIR CHANGE SINCE THE 1870S[END_REF][START_REF] Immerzeel | CLIMATE CHANGE WILL AFFECT THE ASIAN WATER TOWERS[END_REF]Kaser et al., 2010;[START_REF] Sorg | CLIMATE CHANGE IMPACTS ON GLACIERS AND RUNOFF IN TIEN SHAN (CENTRAL ASIA)[END_REF][START_REF] Soruco | CONTRIBUTION OF GLACIER RUNOFF TO WATER RESOURCES OF LA PAZ CITY, BOLIVIA (16• S)[END_REF]. The surface mass balance (SMB) of glaciers is directly driven by the climate conditions; consequently, glaciers are among the most visible proxies of climate change [START_REF] Dyurgerov | TWENTIETH CENTURY CLIMATE CHANGE: EVIDENCE FROM SMALL GLACIERS[END_REF][START_REF] Haeberli | CLIMATE CHANGE AND ITS IMPACTS ON GLACIERS AND PERMAFROST IN THE ALPS[END_REF][START_REF] Oerlemans | GLACIERS AND CLIMATE CHANGE[END_REF]Stocker et al., 2013). Measuring and reconstructing glacier SMB therefore provides critical insights into climate change both at global and regional scales [START_REF] Oerlemans | QUANTIFYING GLOBAL WARMING FROM THE RETREAT OF GLACIERS[END_REF].

Systematic SMB monitoring programmes began in the late 1940s and early 1950s in most of the European countries (e.g. France, Norway, Sweden,

Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data 157 Switzerland). Gradually, more glaciers have become monitored, reaching the present worldwide figure of 440. However, this represents only a small sample of the nearly 250000 inventoried glaciers worldwide (Pfeffer et al., 2014).

Among the existing methods to quantify changes in glacier SMB, the wellestablished glaciological method has become a standard, widely used worldwide, yielding most of the reference datasets (World Glacier Monitoring Service, WGMS; Zemp et al., 2015). Based on repeated in situ measurements, this method requires intensive fieldwork. However, this method is unable to reconstruct the SMB of unmonitored glaciers. The Global Terrestrial Network for Glaciers (GTN-G) aims at increasing substantially the number of monitored glaciers to study regional climate signal through changes in SMB. To reach this objective, the development of methods complementary to the ground-based glaciological method is therefore required. Since the 1970s, several methods have taken advantage of satellite imaging to compute changes in glacier volume [START_REF] Kääb | REMOTE SENSING OF GLACIER-AND PERMAFROSTRELATED HAZARDS IN HIGH MOUNTAINS: AN OVERVIEW[END_REF]Rabatel et al., 2017;[START_REF] Racoviteanu | OPTICAL REMOTE SENSING OF GLACIER CHARACTERISTICS: A REVIEW WITH FOCUS ON THE HIMALAYA[END_REF].

Several glacier surface properties have thus been used as proxies for volume fluctuations -changes in surface elevation from differencing digital elevation models (DEMs) (e.g. Belart et al., 2017;Berthier et al., 2016;Gardelle et al., 2013;[START_REF] Ragettli | Heterogeneous glacier thinning patterns over the last 40 years in Langtang Himal, Nepal[END_REF]Shean et al., 2016); end-of-summer snow line elevation from high spatial resolution optical images (Braithwaite, 1984;[START_REF] Chinn | USE OF THE ELA AS A PRACTICAL METHOD OF MONITORING GLACIER RESPONSE TO CLIMATE IN NEW ZEALAND'S SOUTHERN ALPS[END_REF]Meier and Post, 1962;Mernild et al., 2013;Rabatel et al., 2005Rabatel et al., , 2008[START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF]Shea et al., 2013); mean regional altitude of snow from low spatial resolution optical images [START_REF] Chaponniere | A COMBINED HIGH AND LOW SPATIAL RESOLUTION APPROACH FOR MAPPING SNOW COVERED AREAS IN THE ATLAS MOUNTAINS[END_REF]Drolon et al., 2016); or changes in the glacier surface albedo from high temporal resolution images (Brun et al., 2015;Dumont et al., 2012;Greuell et al., 2007;Greuell and Knap, 2000;Shea et al., 2013;Sirguey et al., 2016). Widely used over icecaps or large ice masses, satellite derived DEMs cannot yet be confidently used to compute annual or seasonal SMB of mountain glaciers, although recent studies have revealed promising results for determining SMB changes of large mountainous glacierized areas (Belart et al., 2017;[START_REF] Ragettli | Heterogeneous glacier thinning patterns over the last 40 years in Langtang Himal, Nepal[END_REF]. The method based on the correlation between the regional snow cover and glacier SMB has shown satisfying results to retrieve seasonal SMB, especially for the winter period. This method was used for the quantification of 55 glaciers SMB in the European Alps over the period 1998-2014 (Drolon et al., 2016). The method based on the identification on high spatial resolution optical images of the end-of-summer snow line altitude has shown encouraging results in the French Alps, multiplying by 6 the available long-term annual SMB time series [START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF], but needs to be automated to compute glacier SMB at regional scales. In addition, monitoring glacier surface properties on the daily or weekly basis and over large glacierized regions is still challenging with high spatial resolution images. The current study is based on the albedo method used in Dumont et al. (2012), Brun et al. (2015) and Sirguey et al. (2016). et al., 2014).

Studied glaciers have been selected following four criteria related to the availability of field data and remote sensing constraints, namely (i) the annual glacier-wide SMB for the study period had to be available, (ii) the glacier surface area had to be wide enough to allow robust multi-pixel analysis, (iii) the glacier had to be predominantly free of debris to allow remotely sensed observations of the albedo of snow and ice surfaces, and (iv) summer SMB records had to be available to consider summer variability. Finally, 11 glaciers have been selected in the Ecrins Range, 14 in Vanoise and 5 in Mont Blanc (Fig. 1, and listed Table 1).

MODIS satellite images

The MODIS sensor, on board the TERRA-EOS/AM-1 satellite has been acquiring near-daily images of the Earth since 25 February 2000. With 36 spectral bands ranging from 0.459 to 14.385µm, and spatial resolution ranging from 0.25 to 1km depending on the spectral band, MODIS is nowadays one of the most used optical sensors for land surface observations. Because of its short temporal revisit time, its long acquisition period and its moderate resolution, images from MODIS are the most suitable for the present work. We therefore rely on about 15000 MODIS calibrated Level 1B (L1B) swath images.

Surface mass balance data

In the French Alps, six glaciers allow both the summer and annual analyses to be conducted, due to the availability of summer SMB data (bs) obtained from in situ measurements with the glaciological method (unpublished data, LGGE internal report, listed in Table 1). In addition, glacier-wide annual SMB values of the 30 studied glaciers were computed by [START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF] using the end-of-summer snow line measured on optical remote-sensing images and the glacierwide mass change quantified from DEM differencing.

Figure 1. Map of the region of interest with the studied glaciers shown in red (numbers refer to Table 1). The four AWS used in the present study were set up on Saint-Sorlin Glacier (no. 20). Adapted from [START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF].

For the six glaciers where glacier-wide annual SMB is available from the two methods, i.e. in situ and satellite measurements, the average of the two estimates was used to calibrate and evaluate the albedo method, in order to

Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data 161 derive for each glacier a single relationship SMB vs. computed albedo. We do not discuss here the differences between the considered datasets because these differences have been investigated by [START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF]. sensor accuracy (±3%, Six et al., 2009), the calculated albedo at the AWS shows a ±10% accuracy (Kipp and Zonen, 2009;Dumont et al., 2012).

METHODS

MODImLab products

MODIS L1B images were processed using the MODImLab toolbox (Sirguey, 2009). Image fusion between MOD02QKM bands 1 and 2 at 250m resolution and MOD02HKM bands 3-7 at 500m resolution allows seven spectral bands at 250m resolution to be produced [START_REF] Sirguey | IMPROVING MODIS SPATIAL RESOLUTION FOR SNOW MAPPING USING WAVELET FUSION AND ARSIS CONCEPT[END_REF]. Then, atmospheric and topographic corrections are applied that include multiple reflections due to steep surrounding topography (Sirguey, 2009). Various products are derived from the corrected ground reflectance including snow and ice surface albedo (Dumont et al., 2012). As recommended by Dumont et al. (2012), the whitesky albedo (estimated value of the surface albedo under only diffuse illumination) is considered. The use of an anisotropic reflection model for snow and ice has been preferred to the isotropic case, due to its closer agreement with in situ measurements (Dumont et al., 2012). The MODImLab toolbox also produces sensor geometrical characteristics at the acquisition time such as the solar zenith angle (SZA) and the observation zenith angle (OZA) used for postprocessing the images (Sect. 7.3.4). The MODImLab cloud detection algorithm is more conservative than the original MODIS product (MOD35), and has been preferred as recommended in Brun et al. (2015).

Table 1. List of studied glaciers, characteristics and albedo/mass balance correlations over 2000-2015, except for summer coefficients (over 2000-2010). For localization, refer to Fig. 1. Bolded rows exhibit glaciers where annual and summer in situ glacierwide SMB data are available. The mask size is expressed in number of pixels. To obtain the glacier mask area in km 2 , one should multiply the mask size by 0.0625km 2 . Determination coefficients are expressed for each glacier. Note the units of r 2 (%), RMSE, P1 and P2 (m w.e.).

Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data 163

According to Dumont et al. (2012) and further assessed by Sirguey et al. (2016) the overall accuracy of MODImLab albedo product under clear-sky conditions is estimated at ±10%.

To mitigate the impact of shadows over the glaciers, MODImLab uses a DEM from the Shuttle Radar Topography Mission (SRTM -90m resolutionacquired in 2000) to estimate the sky obstruction by the surrounding topography and to correct the impact of shadows (see [START_REF] Sirguey | SUBPIXEL MONITORING OF THE SEASONAL SNOW COVER WITH MODIS AT 250M SPATIAL RESOLUTION IN THE SOUTHERN ALPS OF NEW ZEALAND[END_REF].

The algorithm implemented in MODImLab is fully described by [START_REF] Sirguey | SUBPIXEL MONITORING OF THE SEASONAL SNOW COVER WITH MODIS AT 250M SPATIAL RESOLUTION IN THE SOUTHERN ALPS OF NEW ZEALAND[END_REF] and was inspired by Dozier et al. (1981) and Dozier and Frew (1990) for the sky obstruction factor processing (Horizon and Vsky in Sirguey et al., 2016), and [START_REF] Richter | Correction of satellite imagery over mountainous terrain[END_REF] For one glacier in the Alps (Dumont et al., 2012), two in the Himalayas (Brun et al., 2015) and one in the Southern Alps of New Zealand (Sirguey et al., 2016), the summer minimum glacier-wide averaged albedo (  and glacier-wide SMB results from the fact that solar radiation is the main source of energy for melting snow and ice, both at the surface and within the first centimetres below the surface [START_REF] Van As | Warming, glacier melt and surface energy budget from weather station observations in the Melville Bay region of northwest Greenland[END_REF]. But this is not sufficient to explain why averaged surface albedo is suitable for monitoring glacier SMB.

If we consider a temperate glacier in the mid-latitudes, its surface is fully covered by snow in winter, leading to high and uniform surface albedo ( min  ≈0.8 in Cuffey and Paterson, 2010). During the ablation season, the accumulation area is still covered with snow conversely to the ablation area where the ice is exposed and sometimes covered by debris. The overall albedo of the glacier surface is therefore decreasing over the course of the ablation season, providing information on the ratio of these two areas. The ratio between the size of the accumulation zone and the entire glacier, called the accumulation-area ratio (AAR) has often been used as a predictor of SMB both qualitatively (LaChapelle, 1962;Meier and Post, 1962;[START_REF] Mercer | THE RESPONSE OF FJORD GLACIERS TO CHANGES IN THE FIRN LIMIT[END_REF] or quantitatively (Dyurgerov et al., 2009). Therefore, assessing  has been directly correlated with available annual SMB data (listed in Table 1).

Following the work by Sirguey et al. (2016) on Brewster Glacier, a similar approach has been used in order to validate the method at a summer scale but only on six glaciers (within our sample of 30) for which the summer SMB are available. Conversely to Sirguey et al. (2016), the summer SMB bs has been compared to the integrated albedo signal Integrated summer albedos enable us to account for snowfall events that can occur during the ablation period (punctual high albedos). As an example, a strong summer snowfall event leading to a rather persistent snow coverage of the glacier will "feed" the integrated albedo, and physically reduces the glacier melting, which has an impact on the SMB [START_REF] Oerlemans | EFFECT OF SUMMER SNOWFALL ON GLACIER MASS BALANCE[END_REF]. The method therefore accounts for snowfall events to retrieve the glacier summer SMB. To compare each year together and remove the impact of the variable integration time period for each glacier, int s  has been divided by the number of integrated days.

Data filtering

MODIS offers the opportunity to get daily images, but retrieving daily maps of Earth surface albedo remains challenging. Indeed, various sources of error require filtering of the available images in order to only capture physical changes of the observed surface and not artifacts. Clouds are known to be a major problem in optical remote sensing of the Earth surface especially in the case of ice-and snow-covered surfaces. Even though some algorithms exist to differentiate clouds and snow-covered areas (e.g. [START_REF] Ackerman | DISCRIMINATING CLEAR SKY FROM CLOUDS WITH MODIS[END_REF][START_REF] Sirguey | SUBPIXEL MONITORING OF THE SEASONAL SNOW COVER WITH MODIS AT 250M SPATIAL RESOLUTION IN THE SOUTHERN ALPS OF NEW ZEALAND[END_REF], omission errors are difficult to avoid, leading to erroneous albedo of the surface.

In this study, all images with a presence of cloud greater than 30% of the total glacier surface area have been discarded. This threshold is higher than that chosen in Brun et al. (2015) on the Chhota Shigri Glacier (20%), and we thus discuss in Sect. 7.5.1 the impact of the computed cloud threshold on the derived albedo results. When determining min a  , 0% of cloud cover has been imposed as a condition, and a visual check for each year and each glacier has been performed. Snapshots from the fusion of MODIS bands 1-3 and from bands 4-6 [START_REF] Sirguey | SUBPIXEL MONITORING OF THE SEASONAL SNOW COVER WITH MODIS AT 250M SPATIAL RESOLUTION IN THE SOUTHERN ALPS OF NEW ZEALAND[END_REF] have been used to visually check the images, Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data 167 together with images from other satellites (mostly from the Landsat archive) and pictures and comments from mountaineering forums. This last step, although laborious when studying 30 glaciers, allowed the identification of the summer minimum to be improved. A visual check of these images also confirms that projected shadows of clouds are not affecting the albedo maps of summer minimum. Another source of error is the impact of the OZA. As mentioned in Sirguey et al. (2016), accuracy of the MODIS retrieved albedo strongly decreases for viewing angles above 45 • as pixel size increases from 2to 5-fold from OZA=45° to 66° (Wolfe et al., 1998). This phenomenon is accentuated when observing steep-sided snow/ice surfaces, surrounded by contrasted surfaces (rocks, forests, lakes). This distortion could lead to capturing the mean albedo of a glacier plus its surroundings. As a result of this, we decided to filter the images according to their OZA angle, as further described Sect. 7.4.1.

RESULTS

Retrieved albedo assessment

A quantitative evaluation of the retrieved albedo has been performed with AWS deployed on Saint-Sorlin Glacier. Measurements have been synchronized between punctual albedo for MODIS and a 2h averaged albedo around MODIS acquisition time for the AWS. It is worth recalling some differences between the in situ measured albedo data and the one retrieved using MODIS. The downward-facing pyranometer stands at around 1m above the surface, corresponding to a monitored footprint of ca. 300m 2 (theoretical value for a flat terrain) while the pixel area of MODImLab products matches 62500m 2 .

Quantified albedos from each method are therefore not representative of the same area. On the other hand, incoming radiation data are extremely sensitive to a tilt of the sensor located on the AWS, and maintaining a constant angle throughout the monitoring period remains challenging, especially during the ablation season. For instance, a tilt of 5 • of the pyranometer at the summer solstice can increase by 5% the error on the irradiance measurement [START_REF] Bogren | TILT ERROR IN CRYOSPHERIC SURFACE RADIATION MEASUREMENTS AT HIGH LATITUDES: A MODEL STUDY[END_REF]. No sensor tilt was deployed on the AWS, thus preventing the application of tilt correction methods (e.g. [START_REF] Wang | A RETROSPECTIVE, ITERATIVE, GEOMETRY-BASED (RIGB)[END_REF]. Nonetheless, regular visits allowed us to maintain the sensor horizontal and to limit errors in the irradiance measurements. In Fig. 3, the spread between MODIS and AWS albedos is higher for low albedos (i.e. ablation area). This is related to the footprint difference as described earlier, accentuating the albedo differences when monitoring heterogeneous surfaces (snow patches, melt pounds...), which are even more pronounced in summer. One can also note that MODIS albedo often overestimates the AWS albedo value. This overestimation could be explained by ( 1 (residual sum of square). Correlation coefficient per OZA classes are also provided, with r 2 , RMSE together with the number of compared measurements (Nb), and coefficients of the equation: MODISalbedo =P1 MODIS AWSalbedo +P2 MODIS . The continuous grey line illustrates the 1 : 1 relationship between AWS and MODIS retrieved albedo. Thin and dotted lines represent the combined uncertainties on both AWS-and MODISretrieved albedo (absolute value of 10% for each), only accounting for intrinsic sensor accuracy and not for errors related to the acquisition context, e.g. size of the footprint.

Finally, Fig. 3 shows substantial differences between OZA<10 • and other OZA classes. For OZA<10°, MODIS albedos better agree with AWS albedos than for the three other classes. Integrating MODIS images with OZA>10° substantially deteriorates the agreement with AWS albedos (in term of r 2 , RMSE and the slope P1 MODIS ), especially on "narrow" targets as alpine mountain glaciers. We therefore chose to prioritize images acquired with low OZA to avoid detection of non-glacierized surfaces. Therefore, four classes of images have been selected following the criteria presented in Table 2.

For the rest of the computation, the absolute ±10% accuracy per pixel estimated in Dumont et al. (2012) has been considered. We determined the uncertainty on α by accounting for the spatial variability of the albedo signal within the glacier and considering that our sets of pixels are independent from each other (Eq. 2): , (2) where σ stands for the standard deviation (SD) of the pixels albedo with N the number of pixels.

Temporal variability of the albedo signal

Figure 4. The ∼16-year albedo course for Saint-Sorlin Glacier. Glacier-wide averaged albedo is represented with the continuous black line. The green dots spot for each summer the minimum average albedo, and have been manually checked for all years and glaciers. Dashed red and blue lines stand for the beginning of the defined ablation and accumulation seasons (May and 1 October respectively).

Using the "step-by-step" filtering procedure explained in Sect. 7.3.4, the ∼16year albedo cycle of each of the 30 glaciers was obtained. Figure 4 illustrates the entire albedo time series for Saint-Sorlin Glacier over the period 2000-2015. We observed that the albedo decreases from the beginning of summer (dashed red line), reaching min a  in August/September and rising again at the end of September. This cyclicality is a proxy of surface processes. The snow cover decreases at the beginning of summer until reaching its lowest extent, and finally increases again with the first snowfall in late summer to reach its maximum extent in winter/spring. The periodicity of the albedo signal is however not so well defined for some of the studied glaciers. For instance, Argentière Glacier exhibits a severe drop of α in winter, reaching values as low as summer minima (  ∼0.4). The observed drop of albedo in winter occurs during more than 1 month centred on the winter solstice (21 December) and is observed for nine glaciers (Argentière, Baounet, Casset, Blanc, Girose, Pilatte, Vallon Pilatte, Tour and Sélé glaciers). These glaciers are located within the three studied mountain ranges and have the common characteristic of being very incised, with steep and high surrounding faces. We studied the albedo series as a function of the SZA to reveal possible shadowing on the observed surfaces. Figure 5 displays the same cycle as Fig. 4 for Argentière Glacier but providing information about SZA. As a reminder, the MODImLab white-sky albedo is independent of the illumination geometry but the computed albedo for each pixel can be subject to shadowing from the surrounding topography.

Two main observations stand out from the winter part of the cycle in Fig. 5: (i) most of MODIS  severely decrease under  = 0.6 for SZA greater than 60° corresponding to November to January images, and (ii) these drops are not systematic and we rather observe a dispersion cone than a welldefined bias. As there are no physical meanings to systematic change of the surface albedo during a part of the winter period and owing to the fact that this dispersion is only observed for topographically incised glaciers, these decreases in albedo have been considered as artifacts. These observations led us to perform a sensitivity study on the validity of the shadow mask produced by MODImLab, and to study the impact of these shadows on the retrieved glacier-wide albedo (see Sect. 7.5.2).  and ba. However, the linear correlation has no statistical significance for three glaciers with r 2 <0.25. A possible explanation is the high number of removed images in summer due to manually checked thin overlying clouds not detected by the MODImLab cloud algorithm.

Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data 173

Looking at the 27 glaciers for which significant relationships have been found, 2001 is regularly identified as an outlier. According to existing SMB datasets, 2001 is the only year of the period 2000-2015 for which the annual SMB has been positive for all the studied glaciers (+0.80mw.e.yr -1 on average).

To predict correctly the SMB values for the year 2001 using the albedo method, monitored minimum glacier-wide average albedo would need to be extremely high (often greater than 0.7, i.e. 0.83 and 0.95 for Rochemelon and Vallonnet glaciers, respectively), to match the regression line derived from other years of the time series (Table 1). Taking into consideration snow metamorphism during the summer period, melting at the surface and possible deposition of for Saint-Sorlin Glacier. Saint-Sorlin Glacier, together with the five other seasonally surveyed glaciers, showed a significant correlation between the two observed variables (from r2 =0.46 to r2 =0.94 with an error risk<5%, all statistics detailed in Table 1  has been significantly correlated with bs and is therefore a reliable proxy to record the ablation season.

DISCUSSION

In this section, we first discuss the impact of the threshold applied to the cloud cover fraction on the obtained results. Then, a sensitivity study focused on the algorithm correcting the shadows is presented. We finally express the main limitations and assessments of the albedo method.

Cloud coverage threshold

As stated in Sect. 7.3.4, a value of 30% of cloud coverage over the glacier mask has been defined as the acceptable maximum value for considering the albedo map of the day. We computed a sensitivity study on the impact of this threshold on the value of the obtained correlations between the integrated summer albedo and the in situ summer SMB. The summer period has been chosen as it represents the period when the albedo of the glacier is the most contrasted, between bare ice and snow/firn. The glacier-wide average albedo in this period is therefore more sensitive to possible shading of a part of the glacier. Figure 8 illustrates the results for the six seasonally surveyed glaciers.

The used value of the allowed cloud coverage appears not to have a substantial impact on the correlation. This observation implies that the MODImLab cloud product is reliable enough to only compute surface albedo and to avoid too frequent misclassification between the clouds and the surface. It also suggests that removing too many images because of partial cloud cover removes information about the glacier-wide average albedo variability. However, allowing all images, even when the glacier-wide average albedo is computed on only 10% of the glacier (90% of detected cloud coverage), does not reduce significantly the correlation for most of the six glaciers.

Nevertheless, hypothesizing that the glacier-wide average albedo of a small fraction of the glacier (e.g. 10%) is suitable to represent the entire glacierized surface is questionable. It therefore depends on the size of the observed glacier, where 10% of a glacier of 3 and 30km 2 do not have the same meaning, but also on the delineated mask (ablation area not entirely considered because of debris coverage). The summer-integrated albedo is also highly dependent on the time gap between valid useful images. In other words, if an image has an "anomalous" glacier-wide average albedo because of high cloud coverage, the impact on the integrated value will be smaller if "normal condition" albedos are monitored at nearby dates.

Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data 177 value is rather a compromise between the number of used images and the resulting correlation with glacier-wide SMB. We finally concluded that selecting a cloud coverage threshold of 30% presents the best determination coefficients between the integrated summer albedo and the summer balance for most of the six glaciers without losing too much temporal resolution.

Assessment of the impact of shadows on retrieved albedos

In light of the documented dispersion on  during some of the winter months on several studied glaciers (Sect. 7.4.2), sensitivity of the MODIS retrieved albedo against correction of shadows had been assessed. This work has only been conducted on the 250m resolution raster products and specifically on the cast shadow product, because self-shadow corrections can be considered as reliable enough as they are only related to the DEM accuracy. We thus defined a pixel as "corrected" when at least one of its sub-pixels was classified as shadowed. From then on, two glacier-wide albedos  have been defined: (i)

 non-cor computed on non-corrected pixels only, classified as non-shadowed, and (ii)  of both corrected and non-corrected pixels, equal to the glacier- wide average albedo. Figure 9 illustrates the difference between  non-cor and  as a function of the percentage of corrected pixels over the entire glacier.

The study was performed on Argentière Glacier (111 pixels) that exhibited large  artifacts in winter (Fig. 5). The inner diagram allows us to emphasize the annual "cycle" of modelled shadows, contrasted between nearly no cast shadows in summer and an almost fully shadowed surface in winter. We represent the 1SD of  , averaged by classes of 5% corrected pixels. In other words, it illustrates the mean variability of the glacier-wide surface albedo.

Therefore, for images with  non-cor - within the interval defined by 1SD of  , errors resulting from the correction algorithm are smaller than the spatial variability of the glacierwide albedo glacier. We also selected only significant values, following a normal distribution of the averaged  . Consequently, only values at ±1σ (68.2%) in terms of percentage of corrected pixels have been retained (i.e. when the relative share of corrected pixels ranged from 15.9 to 84.1%). Between 0 and 15.9%,  non-cor and α are not sufficiently independent

Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data 179 because of the low number of corrected pixels, and beyond 84.1%,  non-cor is computed over a too small number of pixels. As a consequence, even if the albedo correction in the shadowed parts of the glacier could be improved, most of the errors related to this correction do not depreciate the results. Above 80% of corrected pixels (December to early February), differences between  non-cor and  exceed the monitored spatial variability of  . These anomalies are at the root of the observed artifacts in Fig. 5 caused by the severe drops of albedos and described Sect. 7.4.2.

In addition, a seasonality in the albedo signal can be observed with  non-cor - > 0 in early spring (February to April) and  non-cor - <0 in summer and autumn (June to November). This could be explained by different localizations of shadowed area for a given ratio of corrected pixel.

As an example, a glacier could have in October a snow and shadow-free snout and a fresh snow-covered and shadowed upper section. This configuration would induce a negative difference, as we observe from June to November.

Conversely, this glacier could present in March (same ratio of corrected pixels as October) a complete snow coverage, leading to a smaller difference between  non-cor and  (<0.1), that could even result in a positive difference, as we observe from February to April.

Finally, observed albedo artifacts in winter are most likely due to the correction of shadows. On the other hand, correcting shadows accurately and consistently is extremely challenging. As illustrated by Fig. 9, a way to confidently consider the albedo signal is to exclude values with too large a share of corrected pixels. However, because of the interannual approach carried out in this study, such a systematic artifact is not depreciating the results but would be a major issue in studies focused on albedo values themselves (e.g. maps of snow extent).

Limits of the albedo method

In agreement with Dumont et al. (2012) and Brun et al. (2015), retrieving the glacier annual SMB from albedo summer minima proves to be an efficient method. Low correlations often result from high and persistent cloud coverage during summer, reducing the chance of spotting the albedo summer minimum.

For SMB reconstruction purpose, a future line of research could rest upon linking morpho-topographic features of the glacier, such as glacier surface area, mean altitude or slope, to the regression coefficients of both annual and seasonal SMB vs. albedo relationships, giving the opportunity to establish analogy between monitored and unmonitored glaciers. Tests have been carried

out, but no significant and satisfying results have been obtained, due to a presumably too heterogeneous dataset, where large glaciers (>10km 2 ) and/or

Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data 181 south-facing glaciers are largely under-represented. Larger-scale studies and multi-variable correlations in between morpho-topographic features could be for instance envisaged. Rabatel et al. (2017) recently proposed an alternative approach to reconstruct the annual mass balance of unmonitored glacier on the basis of the albedo method. This approach relies on the ELA method (Rabatel et al., 2005), but using the remotely sensed monitored min a  together with the AAR, the glacier hypsometry, and the regional SMB elevation gradient (which is the annual SMB gradient in the vicinity of the glacier ELA). For an exhaustive description of this approach, see Rabatel et al. (2017).

Using the albedo method for the summer period has shown promising results, with significant correlations found for the six seasonally monitored glaciers. There is still in this approach a step to retrieve the summer SMB of an unmonitored glacier with high confidence.

Table 3. Coefficients of determination for the relationship between the winter SMB bw and the integrated winter albedo, computed with and without the albedo threshold T  .

The winter period has also been considered in the framework of this study, but has not been presented in the main body of this publication because of underwhelming results. The albedo signal between 1 October and 30 April  vs. bw, with and without the use of an albedo threshold T  .

For Argentière and Mer de Glace glaciers, a significant correlation is found whatever the value of the albedo threshold T  . For the four other glaciers, using T  largely improves the correlation. However, T  is far from being uniform on the six glaciers (0.53≥ T  ≥0.76). In addition, for most of the considered glaciers, correlation coefficients abruptly deteriorate when changing this threshold, which does not allow us to use a "regional" threshold for all considered glaciers. On the other hand, Argentière and Mer de Glace without the use of T  provide the best correlation coefficients compared to the other four glaciers; it is noteworthy that they are by far the largest glaciers of our monitoring set (14.59 and 23.45km 2 for Argentière and Mer de Glace glaciers, respectively). With a glacier snout reaching 1600ma.s.l., the tongue Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data 183 of these glaciers can experience melting events (resulting in contrasted pixels in terms of albedo value), even during the winter season. Another difference between our study and Sirguey et al. (2016) is that their work focused only on Brewster Glacier, defined as a maritime glacier. These types of glaciers, even during the accumulation period, can experience strongly varying albedos in their lower reaches, which leads to similar behaviours in winter as for Argentière and Mer de Glace glaciers. We therefore reconsider the idea of Sirguey et al. (2016) to use a threshold as a representative value of fresh snowfall, as there is no physical reason that this threshold varies, at least within the same region. However, an interesting perspective would be to apply the method without threshold, on a set of other maritime or large glaciers (>10km 2 ).

An additional approach has been carried out, aiming at retrieving bw by deduction from the reconstructed ba and bs from the albedo signal. This approach, not using the winter albedo signal, is poorly correlated (r 2 <0.16) with in situ bw for the six seasonally monitored glaciers. Indeed, the result extremely depends on the quality of the correlations between ba, bs and the albedo signals. Saint-Sorlin Glacier is a good example, being one of the glaciers with the highest correlations for the annual (r 2 =0.86) and summer (r 2 =0.94) SMB. Subtracting bs from ba to computed bw leads to an average difference between computed and measured bw of ±0.41mw.e. for the 10 simulated years.

As a consequence, in the case of low correlations between SMB and albedo, errors in the computed winter SMB become exacerbated.

CONCLUSION

In this study, we used the so-called albedo method to correlate annual and summer SMB to glacier-wide average albedos obtained from MODIS images.

This method has been applied to 30 glaciers located in the French Alps, over the period 2000-2015. Image processing has been performed using the MODImLab algorithm, and filters on the images have been applied, removing images with more than 30% cloud coverage, and excluding images with satellite observation angles greater than 30 • . Quality assessment has been performed and close agreement has been found between albedos from AWS installed on Saint-Sorlin Glacier and MODIS retrieved albedo values. Annual SMB has been significantly correlated with the summer minimum albedo for 27 of the 30 selected glaciers, confirming this variable as a good proxy of the glacier-wide annual SMB. For the six seasonally monitored glaciers, summer SMBs obtained from the glaciological method have been significantly linked to the integral of the summer albedo. However, calculating the integral of the winter albedo to quantify the winter SMB as done by Sirguey et al. (2016) has shown underwhelming results. Monitoring winter glacier surface albedo may provide good insight into the frequency of snow accumulation at the surface of the glacier but is poor at quantifying the amount of accumulation. Glaciers that experience complete snow coverage during most of the winter season showed the lowest correlation (r 2 ≤0.33) while the two glaciers showing the best correlations are subject to some events of surface melting in their lower reaches. This approach should not be definitively forsaken, but it requires improvements in order to confidently retrieve winter SMB.

Sensitivity study on the impact of the considered cloud coverage has revealed a high confidence in the MODImLab cloud algorithm, limiting pixel misclassifications, and a rather high tolerance of the integrated signal to the number of partly cloud-covered images. This confidence on cloud filters is very promising to document unmonitored glaciers. Correction of shadows by the MODImLab algorithm has however revealed some limitations when a large share of the glacier is shadowed by the surrounding topography (around winter solstice). Despite this, severe and artificial drops of albedo in winter have not been identified as an obstacle for monitoring both summer and winter SMB.

Such systematic errors are not an issue for inter-annual studies, but would be a serious issue on studies focused on albedo values themselves. For future works, the MODIS archive together with albedo maps, cloud and shadow

Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data 185 masks processed with MODImLab, together with validation data from AWS, offer a unique dataset to monitor the temporal and spatial evolution of the surface albedo of glaciers at a regional scale. For instance, computing the absorbed solar radiation [START_REF] Bair | VALIDATING RECONSTRUCTION OF SNOW WATER EQUIVALENT IN CALIFORNIA'S SIERRA NEVADA USING MEASUREMENTS FROM THE NASA AIRBORNE SNOW OBSERVATORY[END_REF] by date and for each glacier would be an appropriate protocol to estimate the impact of a changing glacier surface albedo in terms of snow or ice melt. Quantifying albedo changes and resulting mass losses with such an approach would be of major interest to better understand the potential effects of possibly increasing dust content, glacier orientation or snow grain growth on glacier surface melt processes.

To conclude, the use of optical satellite images to estimate glacier surface processes and quantify annual and summer SMB from the albedo cycle is very promising and should be expanded to further regions. Using images from different satellites, combining high spatial and temporal resolution instruments, could substantially reduce uncertainties, especially for spotting the albedo summer minimum with more confidence, but also to improve the temporal resolution. This method could then in the short term become reliable for retrieving SMB values of monitored and unmonitored glaciers. Grâce à des lignes de neige délimitées manuellement sur 44 glaciers dans les Alpes françaises et italiennes (Rabatel et al., 2013a[START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF], nous avons pu valider l'aptitude de cet algorithme à identifier l'altitude des TSLA au cours de l'été, et d'en déduire la SLA, représentative de la ELA [START_REF] Lliboutry | Traité de glaciologie[END_REF]. Une limitation de notre méthode réside dans la vérification manuelle Question 2 : l'analyse de cartes d'albédo de surface permet-elle de quantifier sur les glaciers de montagne les processus d'accumulation et d'ablation à l'échelle régionale ?

Grâce au développement d'algorithmes permettant l'estimation de l'albédo de surface de la glace et de la neige [START_REF] Dumont | Monitoring spatial and temporal variations of surface albedo on Saint Sorlin Glacier (French Alps) using terrestrial photography[END_REF](Dumont et al., , 2012;;Sirguey, 2009) et Le second facteur impactant les bilans de masses calculés réside dans les morpho-topographiques des glaciers observés. En effet nous avons pu confirmer une dépendance de la perte de masse cumulée sur l'ensemble de la période d'étude aux variables morpho-topographiques de chaque glacier, confirmant des études précédentes [START_REF] Abermann | Climatic controls of glacier distribution and glacier changes in Austria[END_REF][START_REF] Brun | Heterogeneous Influence of Glacier Morphology on the Mass Balance Variability in High Mountain Asia[END_REF][START_REF] Fischer | Surface elevation and mass changes of all Swiss glaciers 1980-2010[END_REF][START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF] [START_REF] Fitzpatrick | Surface Energy Balance Closure and Turbulent Flux Parameterization on a Mid-Latitude Mountain Glacier, Purcell Mountains, Canada[END_REF]Six et al., 2009;[START_REF] Van As | Warming, glacier melt and surface energy budget from weather station observations in the Melville Bay region of northwest Greenland[END_REF], ou à l'accumulation et à la caractérisation de la variabilité spatiale du manteau neigeux lorsque couplées avec des mesures de type LiDAR ou Radar (e.g. Réveillet et al., 2018;[START_REF] Sold | Methodological approaches to infer end-of-winter snow distribution on alpine glaciers[END_REF]. Les estimations glaciologiques du C'est quelque chose que je me force à faire même si parfois j'utilise à tort ce genre de tournure, pour me faire comprendre par les élèves. Je travaille làdessus pour avoir un discours le plus rigoureux possible avec un minimum de jargon, ou en l'ayant défini auparavant.

 Il manque un bac dans l'expérience. Tu compares neige et neige sale en faisant varier la quantité de "saleté". Mais quand tu compares ces deux avec la glace bleue, tu fais varier deux paramètres, la glace et la couleur.

Peut-être rajouter un bac avec de la glace non colorée ?

C'est une très bonne remarque. Je vais essayer par la suite d'inclure un nouveau bac et de tester cela. Si cela fonctionne, la démonstration sera d'autant plus parlante/rigoureuse.

 Il serait intéressant d'enregistrer la température de surface de la neige pour souligner cet impact de l'albédo.

Par manque de moyen financier je n'ai pas eu accès à des thermomètres infrarouges mais cela serait pertinent.

Tu devrais laisser une trace aux élèves qu'ils puissent revenir plus tard sur ton intervention et ton débat. Par exemple, un texte à trou, ou un petit jeu qui leur permet, s'ils le souhaitent de remobiliser leurs connaissances par la suite.

 C'est vrai que je ne leur laisse qu'une feuille avec des photos de différentes parties du glacier mais rien ne reprenant mon intervention. C'est une excellente remarque complétement adaptée à un public scolaire. Cela a aussi participé à ma réflexion sur le fait qu'une intervention attrayante n'est pas forcément une intervention où l'on aura appris quelque chose, en tout cas de manière durable. Je vais travailler là-dessus pour essayer de préparer une feuille à laisser aux élèves.

Conclusion et perspectives de perfectionnement

Cette

i«

  Avoir la foi, c'est monter la première marche même quand on ne voit pas tout l'escalier » Martin Luther King ii iii

  étudier l'impact du climat local sur les glaciers de montagne ; (2) d'S'INTERESSE-T-ON AUX GLACIERS DE MONTAGNE ?

  7 présentent les articles évoqués précédemment, dont le stade de publication est précisé. Du fait des différents styles imposés par les journaux de publication, certaines redondances pourront être observées entre les articles présentés, je m'en excuse par avance. Finalement, la Sect. 8 clôt ce manuscrit par une synthèse du travail effectué et permet d'aborder les perspectives autour des méthodes et jeux de données développés énigme de l'écoulement des glaciers "[...] quand on agit sur un morceau de glace, qu'on le frappe, on lui trouve une rigidité qui est en opposition directe avec les apparences [...]. Peut-être que les expériences faites sur de plus grandes masses donneraient d'autres résultats ? " Le Chanoine Français Louis Rendu, Théorie des glaciers de la Savoie, 1840 L'intérêt des scientifiques pour les glaciers de montagne est relativement récent. On attribue souvent à Horace Bénédicte de Saussure (1740-1799) et Louis Agassiz (1807-1873) la paternité de la glaciologie même si des observations et des questionnements émergeaient déjà au XVIème siècle.

Figure 1 :

 1 Figure 1 : Figure et légende traduites et extraites de Walker and Waddington,(1988). "Abri construit par Hugi en 1842 sur la moraine médiane de l'Unteraar, sous la confluence du Lauteraar (à gauche) et du Finsteraar (à droite). Cet abri, construit à la base de l'éperon rocheux au centre de l'image a atteint l'endroit représenté en 13 années. Cette moraine a été l'endroit où Agassiz installa l'Hôtel des Neuchâtelois pendant ses recherches en 1840."

2. 1 . 2

 12 Le bilan de masse glaciaire, une variable ad hoc pour observer les glaciers Comme indiqué précédemment, on peut attribuer les premières mesures glaciologiques à De Saussure, Agassiz, Forbes et leurs collaborateurs dans la première moitié du XIXème siècle, même si celles-ci n'étaient pas orientées vers la quantification des variations de volume. Dès 1839, Agassiz met en lien le SMB et la dynamique glaciaire : "[A] glacier is always descending; if the

Figure 2 :

 2 Figure 2 : Reproduction d'un panorama des Glaciers des Bossons et de Taconnaz peint par Samuel Birmann depuis le sommet du Brévent, en 1823. Extrait de Nussbaumer and Zumbuhl, 2012 ("N:454. Samuel Birmann. Ad. nat. f. au sommet du Brévent 1823."; pencil and watercolour; 46.7 × 225.5 cm; Kunstmuseum Basel, Kupferstichkabinett, Inv. Bi.417; photograph by H. J. Zumbühl, with kind permisson of Kunstmuseum Basel)

Figure 3 :

 3 Figure 3 : Figure illustrant quatre glaciers instrumentés par l'IGE dans les Alpes françaises dans le cadre du programme d'observation les GLAciers, un Observatoire du CLIMat (GLACIOCLIM). Sur chacun des glaciers (b-e) le réseau de balises est représenté, avec les triangles rouges et bleus désignant les balises localisées en zone d'ablation et d'accumulation respectivement. Le trait noir désigne le contour de chacun de ces glaciers, et les lignes grises indiquent les courbes d'altitude au pas de 50 m. (a) Carte illustrant la répartition géographique des quatre glaciers représentés (en rouge). (b) Glacier de Saint-Sorlin. (c) Glacier de Gébroulaz. (d) Mer de Glace ainsi que deux glaciers tributaires, Leschaux et Talèfre (désormais déconnecté). (e) Glacier d'Argentière ainsi que deux glaciers tributaires, le Glacier des Améthystes et le Glacier du Tour Noir (TN).Figure extraite de Réveillet et al. (2017).

  SMB durant la période t (e.g. année hydrologique, hiver ou été) sur un glacier de surface tot S , comme la somme des bilans it b mesurés aux balises, et représentative d'une zone de surface i S  (surface relative à la surface totale).

Figure 4

 4 Figure 4 : Schéma illustrant le principe de mesure des variations de surface aux balises. a) représente le protocole utilisé en zone d'ablation, et b) en zone d'accumulation. Schéma dessiné par Horst Machguth.

Figure 5 :

 5 Figure 5 : Schéma d'un glacier hypothétique où plusieurs balises sont représentées. Les balises bleues et rouges sont implantées en zone d'accumulation et d'ablation respectivement. Ces balises sont représentatives de zones délimitées par un trait noir et ont une surface relative à la surface totale du glacier i S  .Figure issue de Letreguilly (2018).

Figure 6 :

 6 Figure 6 : Estimation du bilan de masse de surface à partir de données gravimétriques GRACE sur les hautes montagnes d'Asie, pour la période 2003-2010, issue de Jacob et al. (2012). La figure a correspond à la topographie de la région d'étude. Les points et croix colorés correspondent aux cellules dans lesquelles le changement de masse est estimé (appelé mascon). La figure b exprime les taux de changement de masse de surface pour la période 2003-2010 corrigés des fluctuations hydrologiques et de l'ajustement isostatique de la lithosphère.

  cible et à la pénétration souhaitée dans le milieu) et à mesurer le temps de trajet de cette onde. En connaissant la position exacte

Figure 7 :

 7 Figure 7 : Carte du Sud-est du Groenland montrant les variations de l'altitude de surface (dz/dt) sur la période 2003-2007, exprimées en mètre par an, issue d'une étude de Howat et al. (2008). Les traces rectilignes illustrent la couverture spatiale des données altimétriques acquises par le satellite ICESat. Dans cette étude, ces données ont été complétées par des données issues du capteur ASTER, et sont représentées dans les encadrés illustrant les glaciers exutoires de la calotte.

Figure 8 :

 8 Figure 8 : Taux de variation de l'altitude de surface exprimé en m an-1 sur la période 2000-2009 (A), 2009-2018 (B) et 2000-2018 autour du Mont Reinier (South Cascades, Etats-Unis). Les lignes noires délimitent les contours glaciaires et les zones grisées correspondent à des zones où les données sont manquantes. Figure adaptée de Menounos et al. (2019).

Figure 9 :

 9 Figure 9 : Illustration du système de stéréoscopie ASTER grâce aux deux capteurs ayant des angles de vue différents. Cette combinaison permet une acquisition au nadir (à droite du schéma), puis une deuxième acquisition visant la même fauchée que l'image précédente mais B=361.8 km et 64 seconde plus tard (à gauche du schéma). Cette acquisition stéréoscopique permet d'avoir une représentation en trois dimensions de la surface observée. Schéma issu de Toutin (2002).

Figure 10 :

 10 Figure 10 : Carte d'albédo obtenue grâce au traitement d'une image MODIS du 25 Février 2012 avec l'algorithme MODImLab, permettant de calculer un albédo de surface du Glacier Brewster. Les pixels d'une couleur orange à jaune clair peuvent être attribués à une surface englacée alors que les pixels allant du vert clair au bleu sont représentatifs d'une surface couverte de neige. Carte extraite de Sirguey et al. (2016).

  ELA, la ligne de neige, désigne une ligne visible sur la majorité des glaciers de montagne, et représentant la limite aval du couvert nival sur le glacier. A la fin de la saison d'ablation (exemple Fig. 11), on appellera l'altitude de cette ligne, SLA (pour end-of-summer snowline altitude). Au cours de la saison d'ablation, la ligne de neige migre du front du glacier à sa position maximum en fin d'été, la SLA. Lors de cette remonté graduelle, elle est appelée ligne de neige transitoire (TSL pour transient snowline et TSLA pour l'altitude de la TSL). Cette migration peut être observée grâce à l'imagerie satellitaire optique, comme illustré Fig. 12 sur le Glacier Blanc. Pour les glaciers

Figure 11 :

 11 Figure 11 : Photo du Glaciar Adela Superior et du Glaciar Torre, Argentine. La ligne rouge délimite l'altitude de la ligne de neige à la fin de la saison d'ablation (SLA), prise fin février 2014 (été austral). Elle peut ainsi être considérée comme représentative de la ligne d'équilibre glaciaire.

Figure 12 :

 12 Figure 12 : illustration de la migration graduelle de la TSL du 6 Juillet au 12 Octobre 2016, sur le Glacier Blanc, situé dans le massif des Ecrins, France. Les images ont été acquises par les satellites LANDSAT 8 (*) et Sentinel 2A (+). La couleur provient d'une composition de trois bandes spectrales (bleue, rouge et proche infrarouge) faisant apparaître la glace en bleu foncé, la neige en bleu très clair, et le substratum rocheux ainsi que la végétation dans des couleurs plus naturelles (vert et marron respectivement). L'altitude maximum atteinte par la TSL a été atteinte entre le 2 et le 8 septembre 2016, et son altitude entre ces dates peut être considérée comme la SLA.

3. 2

 2 DEVELOPPEMENT D'UN ALGORITHME PERMETTANT L'IDENTIFICATION AUTOMATIQUE DE LA SLA SUR DES GLACIERS DE MONTAGNE, A PARTIR DE DIFFERENTS CAPTEURS SATELLITAIRES OPTIQUES Article 2 : An automated method to estimate the glacier end-of-summer snowline altitude from optical satellite images Lucas Davaze, Antoine Rabatel, Pascal Sirguey, Yves Arnaud. Soumis et refusé à Journal of Glaciology, annexe de l'article 3, en préparation 1 .Suites aux conclusions et perspectives de l'étude présentée précédemment(Rabatel et al., 2017; article 1), un algorithme capable d'estimer la TSLA et la SLA à l'échelle de plusieurs centaines de glaciers de montagne a été développé et est ici présenté. Les SLAs reconstruites peuvent ensuite être utilisées comme indicateur de la ELA et permettre la quantification du SMB annuel grâce à la méthode ELA. Elles peuvent aussi permettre des interprétations climatiques précieuses, par exemple pour des études hydrologiques à l'échelle de bassins versants englacés ou pour la paramétrisation de modèles glaciologiques.L'algorithme développé repose sur le traitement d'images satellites optiques de résolutions décamétriques (10-30m), adaptées à l'observation de glaciers de montagne (e.g., LANDSAT, ASTER, Sentinel-2). Afin de pallier la topographie de surface parfois complexe et la présence de nuages, différentes corrections (e.g., topographiques, radiométriques…) sont réalisées. Une classification non-supervisée permet ensuite de classer l'état de surface de chacun des pixels (i.e. glace, neige, autre). Après avoir divisé le glacier en bandes d'égale altitude, la TSLA est identifiée, pour chaque image traitée, comme étant l'altitude où la surface du glacier change de majoritairement englacée à majoritairement couverte de neige. La répétitivité des capteurs utilisés permet l'estimation d'une à plusieurs TSLAs par an, la plus haute étant considérée comme la SLA. Finalement, les séries temporelles de SLA déduites sont filtrés des potentielles fausses détections en se basant sur la valeur de chaque SLA par rapport au signal moyen régional des SLAs. Finalement, les SLAs manquantes et filtrées sont comblées en utilisant l'approche bi-linéaire de[START_REF] Lliboutry | Multivariate Statistical Analysis of Glacier Annual Balances[END_REF] permettant d'exploiter à la fois les variabilités spatiales et temporelles de notre jeux de données. Une application de l'algorithme est proposée sur les Alpes françaises et une validation des SLAs reconstruites est effectuée grâce à l'existence de SLAs manuellement tracées. L'algorithme développé permet de reproduire le signal interannuel moyen observé manuellement à l'échelle des Alpes françaises (r 2 =0.83 ; RMSE=82m), pour 34 des 44 glaciers choisis. L'incertitude des SLAs reconstruites automatiquement est de ±123 m. Ces résultats lèvent ainsi le principal verrou identifié (l'automatisation) dans l'article 1 vis à vis de cette méthode et ouvre la porte à l'estimation du SMB annuel à l'échelle régionale. Ce travail a tout d'abord été soumis à Journal of Glaciology le 29 janvier 2019 dans l'idée de présenter une méthode « clefs en main » pour estimer l'altitude de la SLA, développée sur les Alpes françaises. Le 29 mars 2019, notre étude a malheureusement été rejetée à la suite de commentaires proposés par deux des trois rapporteurs. La principale critique reposait sur le manque de disponibilité de l'algorithme (e.g. mise en ligne sur une plateforme partagée), rendant la méthode peu convaincante pour une diffusion et une utilisation par tous. Nous avons ainsi choisi, en prenant les commentaires et corrections en compte, d'inclure ce travail dans une autre étude, visant à quantifier le bilan de masse pluriannuel et annuel de 239 glaciers dans les Alpes européennes à partir d'images satellites optiques. Dans cette étude, nous avons utilisé des estimations automatiques de SLA afin d'estimer le bilan de masse annuel, en utilisant la méthode de la ligne de neige. Ce travail sera donc inclus, en tant que Supplementary Material dans l'article 3, présenté Sect. 6, tout en suivant des recommandations faites par les rapporteurs, de transformer notre premier projet en une étude appliquée sur une région englacée. Par soucis de simplicité, ce travail sera présenté en amont de l'application à l'échelle régionale. des bilans aux balises, et pouvant mener en fonction de la distribution spatiale des balises sur le glacier et des glaciers, à des estimations biaisées.

  dans laquelle la méthode a été conçue et validée. Ce travail a été effectué par Léna Gurriaran, que nous avons accueilli pour un stage de 5 mois pendant son année de césure à l'ENS. Une application sur 82 glaciers répartis sur la Cordillère Blanche et Royale, sur la période 2000-2017 a permis de valider la méthode et l'algorithme sur des glaciers tropicaux, à partir d'images satellites (i.e. ASTER, LANDSAT, Sentinel-2) où 59% de variance commune est en moyenne observée, entre les SLAs dérivées manuellement et semi automatiquement.Une comparaison des SLAs à l'échelle régionale révèle en revanche un biais, menant en moyenne, à une surestimation des SLAs d'environ 100 m. Ce biais n'a malheureusement pu être étudié plus en détails du fait de la durée du stage de Léna. La variabilité interannuelle étant représentée de manière satisfaisante, ce biais n'impacte pas la quantification du bilan de masse dans la mesure où on travaille en écart à l'altitude de la ELAeq, laquelle est calculée à partir de la moyenne des SLA observée.Les séries temporelles de SLA obtenues ont permis de déduire les variations annuelles de SMB grâce à l'estimation du SMB pluriannuel effectuée par Inès Dussaillant à l'aide de la méthode géodésique, basée sur l'imagerie ASTER[START_REF] Gurriaran | Annual glacier-wide mass balances of 82 glaciers in the tropical Andes over a 18-yr period from optical remote sensing data[END_REF]. Le SMB régional moyen pour la période 2000-2017 est -0.40 m.w.e. a -1 , avec des bilans plus déficitaires pour la Cordillère Royale (-0.65 m.w.e. a -1 ) que pour la Cordillère Blanche (-0.33 m.w.e. a -1 ). Les SMBs individuels moyens affichent des bilans majoritairement déficitaires, allant de +0.10 m.w.e. a -1 (Glacier de Safuna) à -1.29 m.w.e. a -1 (GLIMSId : G291837E16266S). Une analyse de ces bilans confirme le lien entre l'index ENSO 1+2 (0-10°S ; 80-90°W) durant les mois de décembre, janvier et février et le SMB annuel dans cette région, déjà établi lors d'étude précédentes sur un nombre plus faible de glaciers (e.g. Vuille et al., 2008). Le pourcentage de variance commune est de 64% et 46% pour la Cordillère Blanche et Royale respectivement. Une étude de l'impact des variables morpho-topographiques sur la période (2000-2017) révèle une forte dépendance du SMB moyen à l'altitude médiane du glacier dans la Cordillère Royale (r 2 =0.54) même si cela n'a pu être observée dans la Cordillère Blanche, notamment à cause d'une forte dispersion des valeurs de SMB pour les glaciers de basse altitude (altitude médiane < 5150 m a.s.l.). Une analyse complémentaire sur les deux cordillères (illustré Fig. 15, révèle une corrélation significative (r 2 =0.23 ; r 2 =0.18, n=82) entre le SMB moyen et respectivement l'altitude maximale et minimale du glacier. Il n'y a pas de corrélation significative entre le SMB et la pente moyenne du glacier (r 2 =0.05) même si une tendance peut être observée, où les glaciers plus raides semblent avoir un SMB moins négatif. Cette étude permet de souligner la capacité de l'algorithme et de la méthode de la ligne de neige à quantifier la variabilité spatiale et temporelle du SMB, à l'échelle régionale, dans une région où la saisonnalité est très différente des Alpes européennes. Le travail de Léna a aussi permis de considérablement améliorer les performances de l'algorithme notamment en 4 ANNUAL AND SEASONAL GLACIER WIDE SURFACE MASS BALANCE QUANTIFIED FROM CHANGES IN GLACIER SURFACE STATE: A REVIEW ON EXISTING METHODS USING OPTICAL SATELLITE IMAGERY Antoine Rabatel 1, *, Pascal Sirguey 2 , Vanessa Drolon 3 , Philippe Maisongrande 3 , Yves Arnaud 1 , Etienne Berthier 3 , Lucas Davaze 1 , Jean-Pierre Dedieu 1 , Marie Dumont 4 1

  applied the ELA method to quantify the annual SMB of 30 glaciers in the French Alps over the period 1983-2014. For this period, only five glaciers in the French Alps had continuous annual SMB time series quantified from in situ measurements. As a consequence, this study multiplied by six the number of long term (> 30 years) SMB time series in the French Alps. Figure 1 illustrates an example of the snow line identification on the Glacier d'Argentière in the Mont-Blanc Massif from a Landsat image.

Figure 1 :

 1 Figure 1: Examples of the snow line identification (in yellow) on the Glacier d'Argentière. The images have been acquired on 18 August 2002 by Landsat-7 (left) and 9 September 2013 by Landsat-8 (right). These illustration uses a spectral band combination involving the bands #2 (0.52-0.60 μm), #4 (0.77-0.90 μm) and #5 (1.55-1.75 μm) for Landsat-7, and the bands #3 (0.53-0.59 μm), #5 (0.85-0.88 μm) and #6 (1.57-1.65 μm) for Landsat-8. 50-m elevation contour lines are shown in gray.

Figure 2 illustrates

 2 the cumulative SMB time series. One can note the very good consistency in the inter-annual mass changes between the glaciers related to homogeneous meteorological conditions and inter-annual variability at the regional scale. This result is in good agreement with the homogeneity in interannual glacier mass changes quantified from in situ point surface mass balance data at the scale of the European Alps[26] and in other regions of the world Annual and seasonal glacier wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery 69[27]. On the other hand, Figure2shows that the cumulative mass loss at decadal time scales varied considerably from one glacier to another, and the authors showed that the changes in glacier surface area and length are not representative of changes in mass balance at this decadal time scale. On the other hand,Rabatel et al. [9] emphasized the shift in the temporal trend of glacier mass loss in the early 2000s, with a 3.25-fold increase in the average mass losses by the 30 glaciers between the periods 1983-2002 and 2003-2014 (-0.47 and -1.52 m w.e. yr -1 , respectively).

Figure 2 :

 2 Figure 2: Cumulative glacier-wide mass balance time series for 30 glaciers in the French Alps (1983-2014), quantified from the end-of-summer snow line. The black curve on the graph is the average of the 30 glaciers and the gray area shows the ± 1 standard deviation interval. Numbers in brackets after the glacier name refer to the map of the right. Adapted from Rabatel et al. [9].

  glaciological monitoring, annual aerial surveys have been conducted since 1977. Oblique photographs captured at the end of the austral summer are used to document the position of the snow line (end-of-summer snow line, corresponding to the SLAi in the previously described ELA method) on 50 "index" glaciers which sample several transects along and across the alpine range. The photographic record supports estimates of SLAi as the annual ELA of glaciers, whose inter-annual variability is used as a surrogate to glacier mass-balance changes [22]. Cullen et al. [38] recently consolidated the longest in situ record of a glacier SMB in New Zealand, namely Brewster Glacier.

Figure 3 :

 3 Figure 3: Relationships between (a) annual and (b) summer SMB and minimum glacier-wide albedo calculated across 26 MODIS pixels on Brewster Glacier (from [42]).

Figure 4 :

 4 Figure 4: (a) Park Pass Glacier with footprint of MODIS 250m resolution pixels used to retrieve glacier surface albedo. (b) Average seasonal albedo cycle of Park Pass Glacier for the 2000-2015 period. Retrieval of albedo is compromised during winter months as the glacier is almost fully in the shade at the time of MODIS/Terra acquisition.

Figure 5 :

 5 Figure 5: (a) Evolution of 𝛼 ̅ 𝑚𝑖𝑛 and 𝑆𝐿𝐴 𝑖 from the snow line aerial survey program for Park Pass glacier over the 2000 -2015 period. (b) Comparative relationship between 𝛼 ̅ 𝑚𝑖𝑛 and 𝑆𝐿𝐴 𝑖 for Park Pass and Brewster glaciers.

  75 m w.e. 100 m -1 across 𝐸𝐿𝐴 𝑒𝑞 . From the hypsometric curve shown in Figure 6a, 𝜕𝑧 𝜕𝐴𝐴𝑅 is found to be ca -2.1 m % -1 . The value of 𝜕𝑏 𝜕𝐴𝐴𝑅 can be estimated as -36.7 mm w.e. % -1 . From Cullen et al. [38], the average annual SMB over the period 2005-2013 is estimated to 𝐵 ̅ = -0.201 m w.e. yr -1 (note that the latter could be obtained from the geodetic method that uses DEM differencing). Using 𝛼 ̅ min values derived by Sirguey et al. [42] and assuming 𝛼 𝑠𝑛𝑜𝑤 = 0.62 (old snow) and 𝛼 𝑖𝑐𝑒 = 0.48 similar to observations made by Cullen and Conway [31],

Figure 7 :

 7 Figure 7: Temporal variations of winter mass balance and cumulative winter albedo for Brewster Glacier.

  sensor. The index has since been widely used with different sensors [e.g.,47- 49]. The NDSI value is proportional to the snow cover fraction inside the pixel and allows an efficient monitoring of the spatial and temporal snow cover variations[50]. InDrolon et al. [51], an NDSI adapted to the SPOT/VGT sensor was proposed, inspired byChaponnière et al. [50], computed from the mean of the red (B0) and blue (B2) channels and from the SWIR band:

Figure 9 :

 9 Figure 9: Example of altitudinal distribution of NDSI for each year since 1998. The NDSI was averaged in a square window centered on Griesgletscher, central Swiss Alps. The red horizontal line represents the NDSI value from which the mean regional snow altitude Z (represented by the red vertical line) is inferred for each year. (a)

Figure 10 :

 10 Figure 10: Observed (a) winter and (b) summer SMB of Griesgletscher, central Swiss Alps, as a function of the mean regional snow altitude Z for each year of the calibration period represented by colored dots. Dashed thin lines represent the 95% confidence intervals for linear regression (solid line). From Drolon et al. [51].

85 Figure 11 :

 8511 Figure 11: Time series of observed SMB (red) and VGT SMB (blue) over the period 1998-2008, averaged for the 55 glaciers. The pink curves represent the ± 1 standard deviation for all glaciers. (a) Winter SMB, (b) Summer SMB. From Drolon et al. [51].

  [e.g., 15-16 and more recently 22] to medium resolution satellite images like MODIS data[e.g., 24,57-58] when the studied ice bodies are wide (i.e. hundreds to thousands of km²).Adaptation of the albedo method to other sensors (currently only MODIS) is also needed. There is good prospect for sustainable observations of glacier surface by upcoming satellite programs such as Sentinel-3 and JPSS constellations, both embarking instruments (e.g, OLCI and VIIRS) capable of providing high temporal revisit frequency at the regional scale similar to MODIS. Image processing currently done for MODIS to retrieve horizontal equivalent reflectances and retrieve surface albedo of snow and ice needs to be adapted to the new spectral resolution. Further adaptation of the albedo retrieval to Landsat-8 and Sentinel-2 data are also needed to research the potential of resolving the spatial distribution of SMB as well as mitigating spectral mixture affecting albedo retrieval from currently coarse pixels. The combination of observations from Landsat-8 and Sentinel-2 may achieve sufficient redundancy to resolve albedo cycle albeit at a higher cost in computing resources.Applying the snow-map method with images from other optical sensors is possible as far as the temporal resolution of image acquisition makes possible the calculation of seasonal snow maps. Such an exercise has already been initiated with observations from the PROBA-V satellite, launched in May 2013 in order to continue the SPOT/VGT mission (which ended in May 2014). Note that a similar exercise is still needed with MODIS and Sentinel-3 data.Regarding PROBA-V, an inter-calibration of the NDSI computed from the PROBA-V 1 km resolution images and from the SPOT/VGT images has been achieved byDrolon [59] over the acquisition overlap period of the two sensors(Oct. 2013 -May 2014). A consistency has been found between winter SMB estimations from the two sensors, in 2013-14. A "near real-time" seasonal SMB monitoring is thus conceivable with PROBA-V or any other daily resolution optical imagery for a large sample of Alpine glaciers, if an inter-calibration of the satellite data can be achieved. Low resolution optical satellite images are available two to three days after the acquisition dates, implying that seasonal snow cover and seasonal SMB could be estimated a few weeks after the end of each season.

( 1 -

 1 100 km 2 ), decametric spatial resolution has been preferred to hectometric to kilometric imagery (e.g. MODIS, SPOT/VEGETATION). Although, our algorithm is suitable for coarser resolution images (e.g. MODIS) when tracking snowline on large glaciers, we developed and applied it to the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), Sentinel-2 A&B and Landsat-5 to -8 images over the 1999-2016 time period. These sensors have been chosen because of their open data policy, revisit times, spatial and radiometric resolutions suitable to monitor mountain glaciers at regional scale. Sensor main features, band denominations and radiometric extents are shown TableS1and Fig.S1. Only orthorectified and radiometrically calibrated images have been used in this study.5.1.2 Digital elevation modelsA DEM is used at multiple stages in our algorithm. Because of its worldwide availability and no penetration into the snowpack (on the contrary of SRTM), ASTER GDEM v2 has been chosen in this study, but one can use any available DEM of sufficient accuracy (  30 m). The ASTER GDEM v2 results in an aggregate of stereo-pair images collected by the ASTER instrument between 2000 and 2009. This last feature could be a source of uncertainty especially on "fast melting" mountain glaciers, as it corresponds to a mean altitude over 2000-2009 and has been therefore taken into account in this study.5.1.3 Glacier outlines and central flow linesGlacial outlines are used in the algorithm to clip the images at the glacier extent. The Central Europe outlines from the Randolph Glacier Inventory 6.0(Consortium RGI, 2017) have been used in the present study.

Fig. S1 .

 S1 Fig. S1. Glacier ice, coarse and fine-grained snow and optical thick clouds spectral reflectance. The bands used in the algorithm are shown for each sensor. TM designates the Thematic Mapper sensor, on-board Landsat 5, ETM+ the Enhanced Thematic Mapper Plus on-board Landsat 7, OLI the Operational Land Imager on-board Landsat 8.

Fig. S2 .

 S2 Fig. S2. Flowchart of the algorithm to estimate SLAs from optical satellite images.Rectangles represent the required datasets except for the last one that correspond to the product of the algorithm ("One SLA per year and per glacier"). Diamonds represent the processing and post-processing steps computed automatically and the

A

  semi-automated method to estimate the glacier end-of-summer snowline altitude from optical satellite images 111 climatic signal exists as demonstrated byVincent and others (2017) and Pelto (2018) by comparing in situ SMB dataset across the European Alps and the North Cascades glaciers, respectively. However, the response to climate forcing can differ for some glaciers in relation with specific morpho-topographic features, dynamics and/or local meteorological conditions, which justify the development of automatic method to monitor large dataset of glacier fluctuations.

  level of 99% (r 2 =0.88) between the detected SLAs on the two different footprints. This comparison indicates that detecting snowlines on a central buffer should not depreciated the results, but introduces an uncertainty in retrieved SLAs of 55 m, corresponding to the root mean squared error (RMSE) of the comparison illustrated in Fig. S5A. The same exercise has been carried on but comparing manually delineated SLAs on whole glacier width and semi-automatically derived SLAs on a central buffer. This comparison is illustrated in Fig. S5B. Deriving the SLA with the algorithm introduces outliers arising from detection errors. Looking carefully at these outliers with altitude differences up to 700 m, it indicates that these discrepancies are caused by the presence of crevasses or change of slopes that have not been sufficiently mitigated by the algorithm. These outliers are mostly detected far above the real SLA, in steep part of the glacier uppermost reaches. Outliers detected close to the glacier snoot are often not kept as SLA if others TSLAs are detected higher during the same summer season. This leads to an observed mean bias of +44 m, largely impacted by the presence of extreme outliers. Despite these outliers, there is a statistically significant relationship at a confidence level of 99% (r 2 =0.56) between manually and semi-automatically derived SLAs. With an RMSE of 123 m for the comparison between manually delineated SLAs on the total glacier width and semi-automatically derived SLAs on a central buffer, we can estimate the uncertainty in retrieved SLAs by the algorithm of 𝛥SLA=± 123 m with 55 m related to difference in detection footprints. Semi-automatically derived 𝛥SLA is in the same order of magnitude than the uncertainty estimated by Rabatel and others, (2013) for manually delineated SLAs (ranging from ±15 to ±170 m depending of the year and the glacier) As manually and semi-automatically derived SLAs are based in the same DEM (ASTER GDEM v2) the encompassed comparison did not account for the DEM accuracy in the estimation of 𝛥SLA. However, regarding the computed uncertainty by the ASTER GDEM validation team (17 m at 95% confidence altitude from optical satellite images 115 level -Tachikawa and others, 2011), it could be neglected in comparison with 𝛥SLA.

Fig. S5 .

 S5 Fig. S5. Error assessment for seven glaciers in the French Alps. Fig. S5A illustrates the comparison between manually delineated SLAs on a central buffer and on the total glacier width. Fig S5B illustrates the comparison between semi-automatically derived SLAs on the same central buffer and manually derived SLAs delineated on the total glacier width. Nb stands for the number of SLAs compared. Heat-map style colours are given as an indicator of the point density, with one point representing one compared SLA.

Fig. S6 .

 S6 Fig. S6. Semi-automatically versus manually derived SLA. Fig S6A illustrates the comparison of SLAs without the filtering and gap-filling steps on 44 glaciers located in the western Alps. Fig. S6B illustrates the comparison for filtered and gap-filled SLA on 34 of the 44 studied glaciers. Statistics are displayed to assess the correlation between semi-automatic and manually derived SLAs. Heat-map style colours are given as an indicator of the point density, with one point representing one compared SLA.

Fig. S7 ,

 S7 Fig. S7, and the annual variability of the retrieved SLAs is presented Fig. S8.As displayed in Fig.S7, the algorithm is able to capture, in a certain extent, the temporal variation of the manually delineated SLAs, with 83% of common variance. The RMSE of this comparison (82 m) is in the range of the previously estimated 𝛥SLA =123 m and in the range of the uncertainty of manually derived SLA (ranging from ±15 to ±170m). The mean observed bias of +48 m indicates that the algorithm tends to overestimate the regional SLA. This was particularly true for 2001 where the filtering and gap-filling steps did not succeed to reproduce the exceptionally low median SLA. 2001 was the only year with positive SMB over the 1999-2016 period with some glaciers experiencing SLA close to the glacier front. In such conditions, the glacier is almost entirely snow-covered, with the SLA lying in some cases beyond the lower boundary of the detection footprint (i.e. 50 m above the terminus).

  SLAs can also be observed for 2007 and 2012 and in a lesser extent for 2006 A semi-automated method to estimate the glacier end-of-summer snowline altitude from optical satellite images 121 and 2013. These discrepancies can arise from different missing data in the two compared dataset, caused by sensor availability and cloud coverage during the acquisitions. Missing SLAs from different mountain ranges can introduce significant bias depending on which mountain range SLAs have not been derived (Rabatel et al., 2013a). Furthermore, most of the discrepancies are observed between 2003 and 2008, period during which SLA estmiates were mostly based on Landsat 7 images, suffering of a failure of the Scan Line Corrector (SLC) resulting in stripped images, reducing drastically the number of available images.

Fig. S7 .

 S7 Fig.S7. Comparison between manually and semi-automatically retrieved regional SLAs (corresponding to the annual median of all SLAs for the two methods). Solid grey line represents the annual median of semi-automatically retrieved SLAs on the 34 glaciers kept after filtering. The solid blue line illustrates the annual median SLA on the same set of glaciers but for manually delineated SLAs.

  detection errors could also occur for exceptional positive SMB years where the SLA is close to the glacier tongue (as observed in 2001) and close to the glacier top (i.e. 2003), beyond the detection footprint. For the same reason, this algorithm must not be applied to derive the SLA on glaciers where the snowline may reach and virtually exceed their uppermost elevation during severe ablation years, but can be used to partially track the TSLA during the ablation season. Another limitation comes from the topographic correction model. The algorithms developed for such a correction can mitigate illumination variations on relatively smooth surfaces but become sensitive to artefacts over steep and rough topographies (steep hanging glaciers, crevasses …) not well captured by the DEM. Improvement of this correction could reduce artefacts due to changing illumination that can compromise the detection. No consensus could however be found in the current literature about the most capable algorithm that could reinforce the robustness of this correction (Minnaert model for Richter and others (2009); slope matching method for Mishra and others (2010); processing scheme for standardized surface reflectance for Ediriweera and others (2013); Meyer's model for Ke and others (2014); sun canopy sensor + cosine correction for Wu and others (2016)). Promising algorithms have been recently developed but are not yet operational and lack in open

6. 3 . 1

 31 Reconstructing multi-decadal mass balance from digital elevation models Digital elevation models (DEMs) have been extensively used to compute multidecadal mass changes of mountain glaciers for entire glacierized regions (e.g.

  et al., 2016;[START_REF] Brun | A spatially resolved estimate of High Mountain Asia glacier mass balances, 2000-2016[END_REF] Dussaillant et al., Under review;[START_REF] Menounos | Heterogeneous Changes in Western North American Glaciers Linked to Decadal Variability in Zonal Wind Strength[END_REF] and thanks to the recent development of two open-source packages (i.e. the Ames Stereo Pipeline (ASP),[START_REF] Beyer | The Ames Stereo Pipeline: NASA's Open Source Software for Deriving and Processing Terrain Data[END_REF] Shean et al., 2016) and the Ohio State University's Surface Extraction from triangulated irregular network-based Search-space Minimization (SETSM;Noh and Howat, 2015, 2017), ~73060 DEMs have been produced from the 30-m ASTER stereoscopic imagery. Produced DEMs are coregistered over stable terrain (excluding ice and lakes) following theNuth and Kaab (2011) approach and using the Global Region-wide annual glacier surface mass balance for the European Alps from 2000 to 2016 131 DEM (GDEMv2,[START_REF] Tachikawa | ASTER Global Digital Elevation Model Version 2summary of validation results[END_REF] as a reference. Because of the ASTER revisit time and global coverage, several DEMs can be produced each year, adding the opportunity to go beyond conventional geodetic approaches that differentiate only two DEMs from two different epochs (e.g.[START_REF] Bamber | A review of remote sensing methods for glacier mass balance determination[END_REF] Gardelle et al., 2013). For each 30-m pixel, a linear regression is fitted to the ASTER-derived DEMs to estimate the temporal trend in elevation change (dh/dt), following the methodology extensively described inBerthier et al. (2016) and[START_REF] Brun | A spatially resolved estimate of High Mountain Asia glacier mass balances, 2000-2016[END_REF]. The overall surface elevation change for each considered individual glacier is then derived using the RGI 6.0 outlines (RGI

6. 3 . 3 1

 331 Climatological data In order to evaluate the role of climate forcings on the observed glacier mass changes, we analyzed the climatic conditions during the period 2000-2016 using data from the ERA5 reanalysis products (Copernicus Climate Change Service (C3S), 2017). ERA5 is a global atmospheric reanalysis produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) that consists in 137 levels from the surface up to a height of 80 km, on a 30 km grid and available since 1979. ERA5 combines a large amount of surface-and satellite-based observations into an estimate at global scale using data assimilation and advanced modelling. In this study we analyzed glacierrelevant variables from ERA5 (i.e. air temperature 2 meter above the grid Region-wide annual glacier surface mass balance for the European Alps from 2000 to 2016 133 mean elevation, total precipitation of the atmospheric column, sea-level pressure, precipitable water, zonal and meridian component of the water vapor and temperature advection into the grid cell). Uncertainty in the annual and decadal retrieved glacier-wide mass balance The uncertainty on individual glacier multi-decadal average mass balance is related to three independent sources: the uncertainty on the glacier surface area, the uncertainty on the rate of elevation change and the uncertainty on volume-to-mass conversion factor. The uncertainty on the rate of elevation change integrated over the area of each glacier is estimated following Rolstad et al. (2009) with spherical variogram of 500 meter range. The standard deviation of stable terrain is drawn in 1x1 degree tiles and applied to glaciers in the corresponding tile, with an average value of 0.5 m a-1. The uncertainty on the glacier surface area is conservatively defined as 5% for clean-ice glacier[START_REF] Paul | On the accuracy of glacier outlines derived from remote-sensing data[END_REF], and the uncertainty on the volume-to-mass conversion factor is estimated as ±60 kg m -3(Huss, 2013). The overall uncertainty is estimated for each glacier, with a mean value of 0.20 m w.e. a -1 , ranging from 0.06 to 0.71 m w.e. a -1 depending of the sub-region.The overall uncertainty in the retrieved annual SMB can be computed as the cumulative of the uncertainty of each variable used in the calculation (equations 1 and 2 from[START_REF] Rabatel | Spatiotemporal changes in glacier-wide mass balance quantified by optical remote sensing on 30 glaciers in the French Alps for the period 1983-2014[END_REF]. It therefore corresponds to the quadratic sum of the uncertainty from:  the uncertainty in the semi-automatically retrieved SLA, equal to ±123 m. It takes into account: i) the DEM resolution and uncertainty; ii) the pixel size and resampling of the satellite images used for the detection and; iii) uncertainties related to the detection method and potential false detections caused by the presence of crevasses, undetected clouds or change of slope not mitigated by the algorithm. A full description is presented in the supplementary materials (Sect. 5 of this manuscript).

Figure 2

 2 Figure 2).

Figure 2 :

 2 Figure 2 : Geodetic multi-decadal average mass balance for the period 2000-2016 as a function of six morpho-topographic variables. Correlation coefficient is provided for each of the variables. A regression has been computed following a robust linear model with support for the M-estimators[START_REF] Susanti | [END_REF], and is represented by a red line, excepted for the glacier aspect, for which a linear regression was not judicious.

6. 5 . 1 Figure 3 :

 513 Figure 3: Map of European Alps illustrating the glacierized and the ice-free (framed cells) 1x1 degree grid cells. For each glacierized cell, the pie charts illustrate the

  The multi-decadal average mass balances for the European Alps over the period 2000-2016 is -0.74±0.20 m w.e. a -1 , corresponding to an average annual loss of -1.54±0.38 Gt, which is less negative than estimates from previous studies: -0.99±0.04 m w.e. a -1 in Huss (2012) computed from in situ data interpolation for the period 2000-2010, or -0.87 ± 0.07 m w.e. a -1 in Zemp et al. (2019) computed from both in situ interpolation and local geodetic estimates for the period 2006-2016.

Figure 5 :

 5 Figure 5 : Mahalanobis distance of the median SMB between each sub-region. A low Mahalanobis distance (i.e. red) corresponds to sub-regions having similar SMB time

Figure 6 :

 6 Figure 6 : Cumulative SMB of glaciers from the European Alps as a function of their mean slopes. The dashed line represents the median cumulative SMB of 23 glaciers obtained from in situ data. Classes 1, 2, 3 and 4 illustrate the cumulative SMB estimated from remote sensing data of four sets of glaciers sorted by their mean slope. Each class has the same number of glaciers but corresponds to a peculiar total glacier area. The class' glacier area of our subset and the corresponding glacier area for the entire European Alps (EA) are displayed in the inset.

Figure S9 :

 S9 FigureS9: Mean centered-reduced summer 2m-air temperature (red) and winter total precipitation in the atmospheric column (dashed olive) variability for the period 2000-2016 over the ten defined 1x1 degree tile region.

Figure S10 :

 S10 Figure S10 : Cumulative SMB of glaciers from the European Alps as a function of their maximum altitude. The dashed line represents the median cumulative SMB of 23 glaciers obtained from in situ data. Classes 1, 2, 3 and 4 illustrate the cumulative SMB estimated from remote sensing data of four sets of glaciers sorted by their maximum altitude. Each class has the same number of glaciers but corresponds to a peculiar total glacier area. The class' glacier area of our subset and the corresponding glacier area for the entire European Alps (EA) is displayed in the inset.

  Images from the MODerate resolution Imaging Spectroradiometer (MODIS) are processed to compute daily albedo maps of 30 glaciers in the French Alps over the period 2000-2015. Then we rely on the methodological framework proposed bySirguey et al. (2016) on Brewster Glacier (New Zealand), looking at the relationships between annual and seasonal SMB and the glacier-wide averaged surface albedo α. Our overall objective is to study the relationships between glacier SMB and albedo by (i) reconstructing the annual albedo cycle for 30 glaciers in the French Alps for the period 2000-2015, (ii) linking the albedo signal to the summer components of the SMB as well as to its annual values for 6 and 30 glaciers, respectively, and (iii) assessing the sensitivity of the retrieved albedo towards tuning parameters (cloud coverage threshold for images processing, reliability of detected shadows). Section 2 presents the available SMB datasets used for the comparison and describes briefly the in situ automatic weather stations (AWS) used to assess the quality of MODISretrieved albedo. The method to retrieve albedo maps is described in Sect. 7.3.Results are presented and discussed in Sects. 7.4 and 7.5. The conclusion gathers the main results of the study and provides perspectives for future works.Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remoteon 30 glaciers located in the French Alps (Fig.1). Each glacier can be classified as a mountain glacier, extending over an altitudinal range from around 1600ma.s.l. (Argentière and Mer de Glace glaciers) to 4028ma.s.l. (Blanc Glacier), and located between the coordinates 44°51" to 46° N and 6°09" to 7°08" E. The cumulative glacial coverage considered in the present study is 136km 2 , i.e. half of the glacier surface area covered by 593 inventoried glaciers over the French Alps for the period2006 -2009 (Gardent (Gardent

7. 2 . 4

 24 In situ albedo measurements Albedo measurements acquired punctually using an AWS on Saint-Sorlin Glacier have been used to evaluate the MODISretrieved albedo. In situ albedo measurements were available for three periods in the ablation zone (July-August 2006; June-August 2008; June-September 2009) and for one period in the accumulation zone (June-September 2008). Albedo data from these AWS have been calculated as the ratio of the reflected to incident shortwave radiation (0.3-2.8µm) using two Kipp and Zonen pyranometers. With a potential tilt of the instrument with respect to surface melting and the intrinsic

  the glacier-wide annual SMB. The relationship between min a

Figure 2 .

 2 Figure 2. Schematic of a typical albedo cycle over one summer, displaying parameters which have been linked to annual and summer (between 1 May and 30 September in the Northern Hemisphere) SMB.int s is retrieved using Eq. (1). The summer minimum value of albedo is represented by



  during the entire ablation season (1 May to 30 September) computed as follows and illustrated in Fig.2:

Figure 3

 3 Figure 3 illustrates the comparison between the retrieved and measured albedos at the AWS locations for various OZA classes. One can note minor differences between the data plotted in Fig. 3 and those presented in Dumont et al. (2012, Fig. 2). These differences are related to changes in the MODImLab algorithm and different computation of the in situ albedo, integrated over a 2h period in the current study.

  ) the MODImLab albedo retrieval algorithm. (Under-estimation of the Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data 169 incoming radiation computed in the MODImLab algorithm would lead to overestimated retrieved albedo values, and in addition the atmospheric corrections used to compute the incident radiation could be hypothesized as source of error -e.g. modelled transmittance through a simplified computed atmosphere; see Sirguey et al. (2009) for further description), or (2) the AWS albedo measurements. Indeed, view angles of AWS pyranometers (170 • ) could influence the retrieved albedo by monitoring out-of-glacier features (e.g. moraines, rock walls,), resulting in underestimated albedo values. However, it is worth noting that most of the points are within the combined uncertainty of both sensors and these differences in albedo retrieved from MODIS and the AWS are thus hard to interpret.

Figure 3 .

 3 Figure 3. MODIS albedo and AWS albedo data for different OZA classes on Saint-Sorlin Glacier. Years indicated in the caption correspond to the year of acquisition while subscripts express the AWS location in the accumulation or ablation areas. The mean discrepancy between MODIS and AWS albedo per OZA is quantified by the RSS

7. 4 . 3

 43 Albedo and glacier-wide SMB min a  and annual SMB The summer minimum average albedo for each year and each glacier has been linearly correlated with the glacier-wide annual SMB. Figure 6 illustrates the relationship between min a  and ba for Blanc Glacier. Error bars show the dispersion of the SMB dataset for each year, and from the glacier intrinsic variability of the albedo signal on the day of

  debris or dusts, monitoring such high albedo values averaged at the glacier scale is unrealistic. As removing 2001 from the time series does not increase the number of glaciers for which the correlation is significant, 2001 has been conserved in the time series. However, this observation reveals a limitation of the albedo method by underestimating the annual SMB value for years with very positive annual SMB.

Figure 5 .

 5 Figure 5. Albedo cycle for Argentière Glacier as a function of the SZA. Each point corresponds to glacier-wide averaged albedo for each available image. The 16 years are displayed. Colour scale gives indication on the date of the used image. The thick

Figure 6 .

 6 Figure 6. Annual SMB as a function of the MODIS-retrieved summer minimum glacierwide average albedo for Blanc Glacier. Error bars show the dispersion of the available annual SMB data and the quadratic sum of the systematic errors made on each albedo measurement. The thin dashed grey line illustrates the line of best fit, along with regression coefficients and significance.int s  and summer SMB Studying the integral of the albedo signal during the ablation season can provide insight into the intensity of the ablation season and thus into the

Figure 7 .

 7 Figure 7. Summer SMB bs expressed as a function of the integrated albedo over the entire ablation season for Saint-Sorlin Glacier. Error bars result from the uncertainties related to the glaciological method (measurements and interpolation at the glacier scale of the punctual measurements, ±0.20mw.e. in total), and on the quadratic sum of the systematic errors made on each albedo measurement. The thin dashed grey line represents the linear regression showing the best correlation between the two variables, together with correlation coefficients.

Figure 8 .

 8 Figure 8. The r 2 for the six seasonally surveyed glaciers for the albedo summer integral vs. summer SMB relationship against the cloud threshold above which images have been discarded during the summer season. For the computation, 100 thresholds have been tested between 0 and 100%. The inner histogram illustrates the number of considered images per summer and averaged on the six glaciers.

Figure 9 .

 9 Figure 9. Impact of the ratio of corrected pixels toward the difference between noncorrected and glacier-wide albedo for Argentière Glacier. Each point corresponds to one acquisition and the 16 years are therefore displayed on this graph. Colour scale gives some indication of the date of the acquired image. Grey shaded areas correspond to ratios of corrected pixels for which  non-cor - has low statistical robustness (refer to the main text). Thin grey lines represent 1σ SD of α, averaged by classes of 5% corrected pixels. The inner graph illustrates the amount of corrected pixel, as function of the selected month.

  has been computed similarly toSirguey et al. (2016) by integrating the winter albedo signal, only when exceeding a certain threshold T  , as described byaw int = a(t)ò if a(t) is found between 10.01 and 04According toSirguey et al. (2016), the use of T  allows for detection of all snowfall events on the glacier by monitoring abrupt rises of  . One of the main conclusions of the above study was the ability of the computed int w  to monitor the frequency of snowfall events, themselves proxy of the accumulation of snow on the glacier, known to be one of the main components of the winter SMB. T  has been chosen to maximize the correlation between the retrieved cumulative winter albedo int w  and the winter SMB. Threshold values have been computed independently for each of the six seasonally monitored glaciers. To evaluate the impact of this threshold, int w has also been computed without threshold over winter months (equivalent to T  =0).

8

  CONCLUSIONS ET PERSPECTIVESLe but du travail que j'ai réalisé pendant ces trois années a été d'explorer la possible automatisation de méthodes permettant d'estimer le SMB de glaciers de montagne afin d'être capable de caractériser leurs variations annuelles et saisonnières pour des glaciers individuels à une échelle régionale. Ce travail est aussi allé au-delà des lignes de codes et m'a permis d'investiguer les causes des variations de masse observées, et de comparer ces résultats aux études précédentes, tout en évaluant la plus-value des jeux de données établis.Je peux désormais apporter des éléments de réponses aux questions identifiées Sect. 2.4.3.Question 1 : Une automatisation de la détection de la ligne de neige à partir d'images satellites optique est-elle envisageable ? Dans la Sect. 5, nous avons développé un algorithme se basant sur l'analyse d'images satellites optiques de résolution décamétrique et issues de plusieurs satellites (Landsat 5, 7 et 8, ASTER et Sentinel 2). L'algorithme bénéficie de modules provenant d'études antérieures permettant, par exemple, d'estimer les ombres propres et portées, et de corriger les images des effets de la topographie(Sirguey, 2009). Cet algorithme constitue une chaine de traitement utilisable sous Python. L'automatisation de cette tache a vraisemblablement constitué le verrou le plus conséquent de mon travail. Il a été assez simple d'établir, avec confiance, l'altitude de la ligne de neige sur quelques glaciers tests. En revanche, développer une méthode générique multi-capteurs, adaptée à l'observation de glaciers ayant une topographie de surface parfois complexe s'est révélé être un vrai challenge, qui explique sans doute pourquoi aucun algorithme n'avait jusqu'aujourd'hui été développé pour une application régionale (excepté[START_REF] Rastner | On the Automated Mapping of Snow Cover on Glaciers and Calculation of Snow Line Altitudes from Multi-Temporal Landsat Data[END_REF] paru en Juin 2019).

  de chacune des images traitées, nécessaire afin d'exclure les nuages nondétectés, la présence de neige fraiche perturbant la détection, ou d'autres artefacts. Cette limitation n'a pas empêché son utilisation avec succès sur 239 Conclusions et perspectives 191 glaciers dans les Alpes européennes et 82 glaciers dans les Andes tropicales.Cet algorithme est aujourd'hui en « concurrence » avec l'algorithme publié très récemment par Philipp Rastner[START_REF] Rastner | On the Automated Mapping of Snow Cover on Glaciers and Calculation of Snow Line Altitudes from Multi-Temporal Landsat Data[END_REF]. Une autre étude est actuellement en cours de publication et est issue d'un travail d'Adina Racoviteanu, avec qui nous avions travaillé sur le massif du Hunza (hautes montagnes d'Asie). Il serait très intéressant de comparer les performances relatives de chacun de ces algorithmes sur une région commune, afin de s'enrichir respectivement des stratégies utilisées, mais aussi de proposer à des utilisateurs une évaluation des qualités et limitations de chacune des méthodes. Dans le cadre de la thèse de Gabriela Collao, nous avions produit des cartes d'albédo à partir d'images MODIS et de la chaine de traitement MODImLab sur le champ de glace Nord de Patagonie. Ces cartes ont été utilisées pour valider des estimations d'albédo faites avec le modèle atmosphérique régional (MAR, Collao-Barrios, 2018). Notre second projet était d'estimer la ligne de neige et le SMB du champ de glace Nord en utilisant la méthode ELA et de le comparer aux estimations faites avec MAR. Cette seconde partie n'a pu être réalisée, faute de temps, mais constituerait un exercice très intéressant, dans le but de pouvoir estimer l'altitude de la ligne de neige sur de grands champs de glace ou calottes à partir de l'analyse de carte d'albédos issues de MODIS. Finalement, l'algorithme que nous avons développé a permis de déduire le premier jeu de données, à notre connaissance, de SMB annuels à l'échelle régionale, issus d'images satellites optiques. Suite à des présentations faites au centre national d'études spatiales (CNES), le produit ligne de neige a été défini comme centre d'expertise scientifique (CES) dans le cadre the THEIA, le pôle national de données sur les surfaces continentales, soulignant la pertinence de cette variable pour la communauté scientifique. Une amélioration notable pouvant être portée à cet algorithme serait son optimisation, le développement d'une interface « user-friendly » ainsi que sa mise en ligne sur une plateforme partagée. Cette amélioration a été suggérée au sein du laboratoire en collaboration avec un ingénieur de recherche mais n'a pas eu de suite jusqu'à aujourd'hui. Cette amélioration permettrait à un utilisateur non-spécialiste d'utiliser l'algorithme et d'y apporter des modifications.

  bénéficiant du cadre méthodologique présenté parDumont et al. (2012) etSirguey et al. (2016), nous avons pu établir des cartes d'albédo issues du capteur MODIS, sur 30 glaciers des Alpes françaises pour la période2000 - 2016 (Sect. 7) (Sect. 7).Cette étude nous a permis de confirmer le lien entre l'albédo moyen minimum à la fin de l'été et le SMB pour 27 des 30 glaciers étudiés, validant ce paramètre comme un marqueur robuste du SMB annuel à partir d'images satellites optiques. L'albédo a ensuite été utilisé après la publication de notre étude sur d'autres régions, comme sur le plateau tibétain[START_REF] Zhang | Annual Glacier-Wide Mass Balance (2000-2016) of the Interior Tibetan Plateau Reconstructed from MODIS Albedo Products[END_REF].L'albédo intégré sur la saison estivale a pu quant à lui être relié au bilan d'été pour les six glaciers étudiés, et permet ainsi de caractériser l'ablation durant l'été.Contrairement à l'étudede Sirguey et al. (2016), nous n'avons pu confirmer, sur les six glaciers étudiés, l'aptitude de l'albédo moyen hivernal à caractériser l'accumulation à la surface du glacier. Ainsi, l'albédo moyenné sur toute la surface du glacier est un excellent marqueur du bilan estival et annuel sans pour autant être capable de l'estimer quantitativement. Des approches se basant sur une relation linéaire entre l'albédo moyen minimum estival et le bilan de masse ont émergé[START_REF] Zhang | Annual Glacier-Wide Mass Balance (2000-2016) of the Interior Tibetan Plateau Reconstructed from MODIS Albedo Products[END_REF] mais sont uniquement applicable sur le glacier étudié. En effet, notre étude révèle des coefficients pouvant varier d'un facteur six d'un glacier à l'autre et des estimations régionales annuelles et estivales de SMB ne sont aujourd'hui pas envisageables en utilisant notre méthodologie. Une nouvelle méthode suggérée par Pascal Sirguey dans la Sect. 4.2.2 pourrait constituer unesolution pour s'affranchir de données in situ mais repose sur une hypothèse forte de l'albédo de la glace et de la neige, particulièrement sur des glaciers de montagne où la topographie peut induire des erreurs dans l'estimation de l'albédo de surface(discuté Sect. 7.5.3).Le lancement des satellites Sentinel 2A&B a permis de considérablement améliorer la couverture temporelle sous nos latitudes avec un temps de revisite inférieur à cinq jours cumulé à l'archive Landsat. Une future direction de recherche pourrait se baser sur l'analyse de cartes d'albédo et/ou de réflectances pour en déduire l'altitude de la ligne de neige transitoire, grâce à un algorithme de détection automatisé (i.e. notre étude,[START_REF] Rastner | On the Automated Mapping of Snow Cover on Glaciers and Calculation of Snow Line Altitudes from Multi-Temporal Landsat Data[END_REF].La variation de cette altitude pourrait ensuite être utilisée pour estimer une variation de masse entre chaque image, et ainsi obtenir une information saisonnière voir subsaisonnière du SMB (e.g.Pelto, 2011). Une telle méthode nécessiterait d'adapter la valeur du gradient de bilan de masse grâce à des campagnes/données in situ répétées permettant d'estimer la dépendance altitudinale du gradient de bilan de masse. Ce projet de recherche a été déposé dans le cadre d'une candidature pour un poste d'ATER qui n'a malheureusement pas été acceptée du fait de l'éloignement des thématiques de recherche ouvertes pour l'année 2019/2020 (Géologie / Géographie et aménagement du territoire). Question 3 : Les méthodes développées sur les Alpes françaises et européennes peuvent-elles être exportées dans d'autres régions pouvant avoir des contextes climatiques différents (e.g. saisonnalité) et où peu de données de validations existent ? Cette problématique de transposabilité de méthodes développées sur des régions géographiques « rassurantes » (proximité géographique des glaciers étudiés, disponibilité de longues séries temporelles de données in situ…) est fondamentale. Malheureusement, les régions les mieux connus, en terme de SMB par exemple, sont souvent celles où les enjeux sont les plus faibles (e.g. ressource en eau, contribution à la hausse du niveau des mers pour les glaciers des Alpes européennes). Notre travail dans les Alpes européennes a permis, profitant des nombreuses mesures in situ de développer et de valider une approche fiable et adaptée permettant l'estimation du SMB pluriannuel et annuel à l'échelle régionale, à partir d'images satellites optiques. En plus d'avoir pu estimer un bilan moyen annuel sur la période 2000-2016 des glaciers des Alpes Européennes (-0.74±0.20 m w.e. a-1), nous avons pu observer une certaine hétérogénéité spatiale et temporelle des séries temporelles de SMB liée aux contextes climatiques locaux (i.e. analyse de données de réanalyses ERA5) et morpho-topographique de chaque glacier. Le travail effectué par Léna, dans la continuité des travaux de Rabatel et al. (2012), Loarte et al. (2015) et Soruco et al. (2009), nous a permis de confronter cette méthode développée dans les Alpes sur une région à fort enjeu, les Andes tropicales (cordillères Blanche et Royale). Cette étude (présentée en Annexe 1) nous a permis de valider notre approche sur un massif où le contexte climatique est très différent des Alpes européennes (accumulation et ablation possible toute l'année, dépendance très forte au mode de variabilité climatique principal de l'Océan Pacifique : ENSO). Les résultats obtenus ont permis de quantifier un SMB annuel sur 82 glaciers dans ces deux massifs. Cette étude nous a finalement montré qu'une analyse à l'échelle régionale ne prive pas l'opportunité d'observer des curiosités à l'échelle d'un glacier, comme l'avancée du Glacier Safuna, ayant connu une crue dans les années 2000, probablement suite un déstabilisation d'une partie du glacier entrainant son avancée, et dont la littérature scientifique ne fait aujourd'hui, pas état à notre connaissance. Question 4 : L'estimation du SMB à l'échelle régionale apporte-t-il réellement une plus-value sur la compréhension des causes à l'origine des variations observées ? Les précédentes estimations issues de données in situ sont-elles représentatives de la majorité des glaciers de la région étudiée ? Grâce au développement d'un algorithme de détection semi-automatique de la ligne de neige, nous avons pu appliquer la méthode ELA pour quantifier le bilan annuel de 239 glaciers dans les Alpes européennes (Sect. 6) et 82 glaciers dans les Andes tropicales (Annexe 1). Nous avons pu identifier deux facteurs impactant les bilans calculés. Le premier est le forçage climatique local qui impacte le bilan de masse glaciaire. Dans les Alpes européennes, nous avons utilisé des données de réanalyses climatiques (i.e. ERA5) dans le but d'estimer les variables climatiques impactant le plus les bilans calculés. Cela nous a permis de discerner deux régions où les glaciers ont une variabilité interannuelle des bilans de masse contrastée. Dans la première région, la variabilité semble majoritairement impactée par les précipitations hivernales (Nord-Est de la chaîne), alors que dans la seconde, elle semble assujettie aux flux d'humidité provenant de l'Ouest et aux conditions météorologiques estivales (versant Ouest de la chaîne).Dans les Andes tropicales, des résultats préliminaires confirment (e.g.[START_REF] Vuille | Climate change and tropical Andean glaciers: Past, present and future[END_REF] que la variabilité interannuelle du SMB semble impactée par les événements marqués de l'ENSO. Une étude plus poussée serait nécessaire, en combinant les données de SMB avec des données climatiques Conclusions et perspectives 195 issues de réanalyses (i.e. ERA5) ou de modèles physiques (i.e. WRF, RACMO, MAR) dans le but de mieux appréhender les variabilités spatiales et temporelles.



  SMB ne permettent en revanche pas d'observer un nombre important de glaciers du fait des contraintes humaines, logistiques et financières. Les études que nous avons réalisées (présentées Sect. 4, 5, 6 et 7) ont, en revanche, permis de valider l'utilisation d'images satellitaires pour quantifier le SMB estival et annuel, à l'échelle régionale. Nos résultats (Sect. 6) montrent que les SMB estimés sont en accord avec les SMB issus de données in situ, uniquement pour une partie des glaciers observés par télédétection. En effet, les glaciers « idéaux » pour des campagnes in situ[START_REF] Kaser | A manual for monitoring the mass balance of mountain glaciers[END_REF] n'étant pas représentatifs de la majorité des glaciers, notamment du fait de leurs caractères topographiques facilitant leur accès (peu pentus, parfois à relativement basse altitude), certains biais peuvent apparaitre. Dans les Alpes Européennes, les glaciers ayant une pente moyenne supérieure à 20° ne sont pas représentés par les données in situ alors que leurs SMB cumulés est significativement moins négatif. De fait, les études estimant le SMB régional à partir d'extrapolation d'observation in situ doivent nécessairement prendre en compte une dépendance au contexte morpho-topographique des glaciers étudiés (Huss, 2012) ; Perspectives personnelles Ces trois années de thèse ont été pour moi une expérience très enrichissante autant sur le plan humain que scientifique. J'ai appris ce que c'était que de se concentrer sur une question (ici scientifique) pendant trois ans, au-delà de l'ascenseur émotionnel, dans lequel chaque doctorant monte lors de la signature de son contrat doctoral. J'ai aussi eu l'opportunité d'enseigner au cours de ces trois années, autant à l'université qu'avec un public scolaire et tout public (listes des enseignements et présentation d'une réflexion pédagogique en Annexe 2 et 3). Ces expériences ont été très enrichissantes et j'ai énormément apprécié les contacts que j'ai eus avec le public. Mes interventions étant très orientées « climat », je me suis rendu compte des attentes des auditoires envers la compréhension des mécanismes climatiques à l'origine du changement actuel, mais aussi du travail de sensibilisation encore nécessaire pour favoriser une prise de conscience citoyenne, indispensable, à mon sens, pour favoriser des actions rapides d'atténuations et d'adaptations. C'est pourquoi je souhaite, au terme de cette thèse, me consacrer à plein temps à des activités de médiation scientifique autour du climat, auprès du grand public, d'enseignants, d'élus locaux, ainsi que d'un public n'ayant normalement pas accès à des interventions sur le climat. Au vu de l'effectif restreint et de leur intérêt pour l'environnement de montagne, c'était presque gagné. Et pourtant, les élèves ont été assez dubitatifs, et je l'avoue peu intéressés. Les deux expériences étaient introduites par des feuilles type TD, l'une assez scolaire sur la fonte de la neige, et une autre sous la forme d'une chasse au trésor sur un glacier. Selon moi, cela n'a pas fonctionné car : Le thème sur l'écoulement des glaciers (exemple d'une confluence glaciaire), malgré les efforts de mise en place, ne me semble pas fondamentale mais plutôt anecdotique à leur niveau.  A la manière avec laquelle elle était présentée, les étudiants n'ont pas identifié la notion essentielle que je souhaitais introduire. Lors de ce même déplacement, je suis aussi intervenu sur une conférence grand public d'1h30 sur le réchauffement climatique et son impact sur les glaciers. Une 60aine d'adultes étaient présents. Cette conférence s'est bien déroulée (auditoire attentif, questions et discussions après la présentation). Assez satisfait, ce n'est que quelques semaines/mois plus tard que je me posais les questions suivantes :

Figure

  Figure A2.16 Expérience "Albédo"

  "patauge pédagogique" ainsi que mes interventions et réflexions autour du débat scientifique ont été très intéressantes. Cela a grandement participé à mon souhait de poursuivre à plein temps une activité autour de la diffusion des savoirs, en salle et en montagne. De par le temps impartit et mon activité de recherche, je n'ai pas pu passer autant de temps que je l'aurais souhaité sur le développant d'interventions et la mise en place de ces débats. Je souhaiterais ainsi perfectionner mes interventions grand public et les focaliser sur un ou deux processus à décrire pour que les gens en aient compris le sens profond, ou à défaut, aient identifié l'obstacle épistémologique et l'intérêt du processus. J'ai fait d'autres essais de débat non discutés dans ce portfolio (écoulement du glacier, importance de la couverture détritique, inversion de température en météorologie) car j'ai moins préparé ces débats et/ou je les ai réalisés peu de fois. Je souhaite donc développer d'autres débats, toujours en lien avec le thème "Climat et Glaciers". Lors des interventions où le débat a vraiment "pris", j'ai été fasciné par l'émulation et la volonté commune de trouver une solution au problème posé, en utilisant au maximum toutes les compétences que les élèves ont réellement acquises (sens profond). Je trouve cette méthode d'une efficacité redoutable lorsque le débat est bien posé, sur un théorème/conjecture/loi dont la compréhension est difficile à appréhender (obstacle épistémologique), et qu'il est primordial de le comprendre pour la discipline enseignée.

  .

	Les premières décennies du XIXème siècle marquent le début de la
	glaciologie	"moderne".	Les scientifiques passent	d'une	approche
	principalement "contemplative", basée sur de rares expériences et mesures, à
	une approche quantitative. Même si De Saussure propose dès 1787 des
	estimations d'épaisseur de glace sur la Mer de Glace, des mesures de
	température et d'humidité de l'air ainsi que de la pression atmosphérique

jusqu'au sommet du Mont-Blanc, les premières campagnes de mesures glaciologiques datent de 1827. Par exemple, les premières mesures de vitesse d'écoulement en surface sont effectuées par Franz Joseph Hugi sur l'Unteraargletscher dans les Alpes Bernoises, en Suisse

  . Utilisée notamment par le service d'observation des glaciers français (Les GLACIers, un Observatoire

	du	CLIMat,	GLACIOCLIM)	cette	approche	peut	être	appliquée
	rétrospectivement sur des séries de bilans de masse existantes (e.g. Basantes-

  diminue drastiquement devant les variations d'épaisseurs observées. De plus, il n'existe aujourd'hui pas de satellites capables d'avoir une couverture et une répétitivité suffisante pour produire des MNTs de surface chaque année, de résolution métrique et à l'échelle globale. La méthode géodésique est en revanche fiable pour quantifier les SMB moyens décennaux de glaciers de

	montagne sur des régions englacées entières, où peu de mesures existaient
	(Hautes Montagnes d'Asie : e.g., Brun et al., 2017; Kääb et al., 2015; Andes :
	Dussaillant et al., 2019; Western North America : Menounos et al., 2019;
	plusieurs régions englacées : Zemp et al., 2019). Ces estimations ont ainsi
	permis de quantifier la contribution des glaciers de ces régions à la hausse du
	niveau des mers. Cette méthode est aussi utilisée en complément d'autres
	approches capable d'estimer le SMB annuel, dans le but de contraindre les
	résultats obtenus pour éviter le cumul d'erreurs ou de biais (méthode
	glaciologique,

  , cet indicateur est peu utilisé et des méthodes se basant sur l'albédo de surface ou la position de la ELA sont aujourd'hui préférées.

	Méthodes
	période d'ablation. Son albédo varie entre 0.2 et 0.64 en fonction de son La multiplication des capteurs satellitaires optiques permettant le calcul
	degré d'impureté (Cuffey and Paterson, 2010). de l'albédo de surface (e.g.
		LANDSAT, Sentinel, AVHRR,
	4. Les débris ou couverture détritique. Ils sont caractérisés par une MODIS) a permis d'observer
	absence de tri granulométrique, caractéristique des dépôts glaciaires presque quotidiennement le
	tels que les moraines. Les débris sont pour la plupart d'origines changement d'état de surface
	minérales et proviennent des écroulements rocheux en zone des glaciers, favorisant
	d'accumulation, des produits de l'érosion glaciaire des moraines latérales l'émergence de méthodes
	dans la zone d'ablation, de dépôts éoliens ou plus rarement de l'érosion utilisant l'albédo comme
	sous-glaciaire. Ils peuvent aussi être d'origine biologiques (Anesio and marqueur du SMB. Dès 2000,
	Laybourn-Parry, 2012). Ces débris peuvent être incorporés à la neige, Greuell et Knap utilisaient
	au névé et à la glace, modifiant leur albédo de par leurs couleurs l'albédo pour détecter
	généralement sombres. Certains glaciers peuvent être partiellement ou l'altitude de la "slush-line" au
	intégralement couvert de débris ce qui peut avoir un fort impact sur leur Groenland (ligne délimitant la
	bilan de masse (Brun et al., 2019).	partie recouverte de neige de
		la partie où la fonte peut faire
		apparaître de l'eau liquide par
		endroit, similaire à la ligne de
		neige sur les glaciers de
		montagne) à partir d'image
		AVHRR. Suite à des difficultés
		rencontrées pour délimiter
		avec certitude la ligne de neige
		sur la calotte du Vatnajökull en
		Islande, De Ruyter De Wildt et
	compaction. Le névé peut être apparent à la surface du glacier si la neige al. (2002) modifient
	accumulée par-dessus pendant l'hiver a complètement fondu. Ce l'approche de Greuell et Knap
	phénomène est de plus en plus commun, notamment du fait de bilans en calculant un albédo moyen
	de masse de plus en plus négatifs. Son albédo peut varier entre 0.30 et  sur tout le glacier. Cet albédo moyen est caractéristique de la fraction du 0.55 (Cuffey and Paterson, 2010) en fonction du degré de glacier couverte de neige. Il est ensuite possible de suivre les variations de  métamorphisme et de la quantité d'impuretés contenue dans le névé. à la surface du glacier au cours de la saison d'ablation en comparant des cartes
	d'albédo, comme celle illustrée Fig. 10 sur le Glacier Brewster, en Nouvelle-3. La glace est le stade final d'évolution de la neige/névé. En théorie Zélande (issue de Sirguey et al. 2016). Cet index, représentatif de l'AAR, est invisible pendant la période d'accumulation car masquée par le manteau intimement relié aux processus de surface. Ainsi, une année à bilan de masse neigeux, la glace est visible sur toute la zone d'ablation à la fin de la très négatif mènera à une réduction de la zone d'accumulation du glacier
	27

L'albédo de surface et le bilan de masse de surface L'albédo de surface du glacier est une caractéristique essentielle du bilan d'énergie de surface car il gouverne la quantité d'énergie de rayonnement absorbée par la surface. On distingue généralement à la surface du glacier quatre états de surface possibles :

1. La neige. Elle recouvre la majorité du glacier pendant la période d'accumulation et seulement la zone d'accumulation en fin de période d'ablation. En fonction du métamorphisme et de la quantité d'impuretés ou de débris qu'elle contient, la neige peut avoir un albédo compris entre 0.4 et 0.85

(Cuffey and Paterson, 2010)

.

2. Le névé correspond à de la neige ayant subsisté au moins une année sans fondre totalement, et ayant subi du métamorphisme et de la

moyenne des ELAs mesurées in situ sur la période est exprimée par le point "ELA mean 1967-88". Figure extraite de Hagen et al., (1991) et modifiée par

				. Cette
		méthode se base sur l'estimation de
		la différence entre la ELA de l'année
		et la ELA du glacier si son bilan était
		équilibré (appelé ELA0 et représenté
		Fig. 13 par l'intersection de la droite
		avec l'axe des ordonnées). Afin de
		déterminer le SMB, cette différence
		est multipliée par ce que Braithwaite
		appela le gradient de bilan effectif.
		Initialement	considéré	comme
		constant sur une région entière,
		Braithwaite suggéra que le gradient
		de bilan de masse pouvait varier
	Figure 13 : Lien entre la ELA et le bilan de masse annuel de surface pour le	d'une région à l'autre et était relié au gradient de bilan de masse de surface,
	glacier de Broggerbreen, de 1967 à	exprimé en m.w.e. m -1 . Ce gradient
	1988. La	exprime, le taux de variation du bilan
		de masse de surface en fonction de
		l'altitude au sein du glacier. Ce
		gradient de bilan est souvent calculé
	à partir de données in situ. Si l'on compare par exemple le gradient de bilan
	de masse calculé sur des glaciers tropicaux dans les Andes (1.2 m 100 m -1
	d'après Rabatel et al., 2012) et celui calculé sur des glaciers des Alpes
	françaises (0.78 m 100 m -1 d'après Rabatel et al., 2005), le bilan de masse de
	surface semble varier plus fortement avec l'altitude dans les Andes que dans
	les Alpes françaises. Ce gradient est évidemment influencé par des facteurs
	locaux comme la variabilité de l'accumulation de neige due à la topographie
	de surface, ou la quantité de débris à la surface, même si les hypothèses de
	gradients uniformes et constants à l'échelle d'une région mènent à des
	résultats très bien corrélés avec les données in situ (Benn and Lehmkuhl, 2000;

Tableau 1 :

 1 Tableau synthétisant les capacités des méthodes présentées Sect. 2.3 à déterminer le SMB à plusieurs échelles spatio-temporelles. La méthode glaciologique est encadrée et considérée comme référence. D/I spécifie si la méthode est directe ou indirecte (détails Sect. 2.3). « Méthode calibrée » renseigne sur la dépendance de la méthode à utiliser des données in situ et ainsi, son habilité à être utilisée sur des régions/glaciers où peu ou pas de données in situ existent. * indique qu'une caractérisation de la densité de surface est nécessaire pour déduire le SMB. ⁺ indique que cette méthode requiert l'estimation du gradient de SMB, obtenu à partir de données in situ. Les croix indiquent que des études ont permis une estimation du SMB à l'échelle spatio-temporelle considérée alors que les zones grisées reflètent l'absence ou l'utilisation anecdotique de la méthode à cette échelle.

2.4.

2 Méthodes permettant le calcul du SMB saisonnier L'estimation du SMB saisonnier consiste à estimer le SMB au cours de la période d'accumulation (hiver) et/ou au cours de la période d'ablation (été). Cette estimation permet de mieux contraindre les contributions des glaciers

  Annual and seasonal glacier-wide surface mass balance quantified from changes in glacier surface state: a review on existing methods using optical satellite imagery Antoine Rabatel, Pascal Sirguey, Vanessa Drolon, Philippe Maisongrande, Yves Arnaud, Etienne Berthier, Lucas Davaze, Jean-Pierre Dedieu and Marie Dumont. Remote Sensing 9, 507 (2017).

	densité et d'une très forte accumulation due au contexte climatique de la
	région étudiée réduisant ainsi l'impact des incertitudes évoquées Sect. 2.3.
	Cette étude reste anecdotique et une application à plus grande échelle et dans
	une région à plus faible accumulation n'est aujourd'hui pas envisageable. En
	revanche, la méthode basée sur les cartes d'enneigement a été appliquée avec
	succès par Drolon et al. (2016) sur 55 glaciers dans les Alpes européennes,
	permettant l'estimation du SMB hivernal, sans pour autant parvenir à estimer
	de manière satisfaisante le SMB estival (détails Sect. 2.3).
	Utilisant l'évolution de l'albédo à la surface du glacier, Sirguey et al.
	(2016) ont proposé une méthodologie adaptée à l'estimation du SMB estival
	et hivernal, avec une étude focalisée sur le Glacier Brewster (Nouvelle
	Zélande). Malgré sa dépendance à des données in situ, la méthode albédo
	semble être adaptée à une application régionale et les verrous liés au
	traitement massif d'images pour produire des cartes d'albédo ont été levés par
	Sirguey (2009), Sirguey et al. (2016) et Dumont et al. (2011, 2012) grâce au
	développement d'algorithmes de traitement et de corrections d'images.
	2.4.3 Méthodes privilégiées et objectifs de thèse
	L'objectif principal de mon travail de recherche effectué pendant cette thèse a
	été de développer et d'appliquer des méthodes permettant d'estimer le SMB
	de glaciers de montagne à une échelle régionale, annuelle et saisonnière, à
	partir d'images satellites. Cet objectif découle d'un manque de données
	d'observation à l'échelle régionale, limitant notre compréhension de l'impact
	du dérèglement climatique au sein de régions englacées où peu ou pas de
	données existent. Améliorer la couverture spatiale de la quantification du SMB
	au cycle hydrologique dans des bassins versants fortement englacés, ainsi que
	de mieux comprendre les raisons de la variabilité interannuelle observée sur
	des glaciers de montagne.
	Parmi les méthodes décrites précédemment, peu permettent
	l'estimation du SMB saisonnier. Les récentes améliorations en termes de
	résolutions spatiale, radiométrique et temporelle ont permis de rares
	applications de la méthode géodésique à l'échelle saisonnière, notamment en
	bénéficiant d'images Pléiades et Worldview2. Une étude, menée par Belart et
	al. (2017), a conduit à la quantification du SMB hivernal du champ de glace
	Drangajökull (Islande) pendant l'hiver 2014-2015, mais n'a été rendu possible
	que par la disponibilité d'images très haute résolution, de mesures in situ de

de glaciers de montagne permettrait, de plus, de corriger certains biais de non représentativité issus de SMB mesurés in situ par la méthode glaciologique

(Gardner et al., 2013)

. La méthode géodésique permet de pallier ce verrou de nonreprésentativité spatiale grâce à sa capacité à estimer un changement de masse cumulé sur une période d'environ 10 ans. Elle ne permet en revanche pas d'analyser la variabilité annuelle, essentielle pour caractériser l'impact de certaines conditions climatiques et météorologiques sur le SMB. Bénéficiant d'une synergie entre l'amélioration d'outils informatiques permettant l'utilisation de la méthode géodésique à une échelle globale (i.e. projet ASTERix) et l'augmentation du nombre de capteurs optiques de résolutions décamétriques adaptés à l'observation des glaciers de montagne (e.g., Sentinel 2, LANDSAT, ASTER), ) et la mise à disposition libre des données satellites nous avons choisi d'utiliser la méthode ELA dont les principaux verrous résidaient dans l'automatisation de la détection de la ligne de neige, la disponibilité d'un nombre suffisant d'images au cours de l'été et la disponibilité de SMBs moyens pluriannuels issus de la méthode géodésique. Dans un premier temps, nous nous sommes donc concentrés sur l'automatisation de la détection de la ligne de neige à partir d'images satellites optiques. En découle ainsi une première question de recherche : Question 1 : Une automatisation de la détection de la ligne de neige à partir d'images satellites optique est-elle envisageable ? Poursuivant cette volonté de perfectionner notre compréhension des processus climatiques et météorologiques à l'origine des changements de masse glaciaire observés à l'échelle saisonnière, une deuxième méthode a été privilégiée. Ainsi, une deuxième partie de mon travail de recherche s'est concentrée sur le développement et la validation de la méthode utilisant l'albédo de surface comme marqueur des processus de surface. Comme illustré précédemment, cette méthode permet d'estimer le SMB à un pas de temps saisonnier et bénéficie du développement d'algorithmes de traitement d'image, rendant l'estimation du SMB à une échelle régionale envisageable, sans pour autant avoir été effectuée : Question 2 : l'analyse de cartes d'albédo de surface permet-elle de quantifier sur les glaciers de montagne les processus d'accumulation et d'ablation à l'échelle régionale ? Le développement des deux méthodes choisies s'est concentré sur des glaciers de montagne localisés dans les Alpes françaises. Ce choix provient de l'existence de données météorologiques ainsi que de SMB annuels et saisonniers sur six glaciers français, obtenues par la méthode glaciologique et maintenus par le réseau d'observation GLACIOCLIM. De plus, des lignes de neige ont été manuellement délimitées sur 44 glaciers dans les Alpes occidentales dans le cadre du travail de Rabatel et al. (2013a), et ont abouti à l'estimation du SMB annuels de 30 de ces glaciers (Rabatel et al., 2016). Ces données ont ainsi pu être utilisées pour valider les lignes de neige estimée automatiquement ainsi que les SMB obtenues. Question 3 : Les méthodes développées sur les Alpes françaises et européennes peuvent-elles être exportées dans d'autres régions pouvant avoir des contextes climatiques différents (e.g. saisonnalité) et où peu de données de validations existent ? L'automatisation de certaines méthodes, couplée avec l'augmentation du nombre de capteurs satellitaires opérationnels ouvre aujourd'hui la voie à des estimations de variations de masse glaciaire sur des régions englacées entières. Au-delà du challenge technique, scientifique et technologique, une autre question est apparue au cours de ma thèse : Question 4 : L'estimation du SMB à l'échelle régionale apporte-t-il réellement une plus-value sur la compréhension des causes à l'origine des variations observées ? Les précédentes estimations issues de données in situ sont-elles représentatives de la majorité des glaciers de la région étudiée ? 3 RESUME DE MES RECHERCHES Les efforts de la communauté glaciologique pour quantifier le changement de masse des glaciers de montagne, notamment à partir d'images satellites, ont permis l'élaboration de méthodes novatrices (détaillé Sect. 2). En revanche, peu d'applications ont été menées pour estimer le SMB annuel et saisonnier à l'échelle régionale. Ma première motivation scientifique a été d'adapter et de développer ces méthodes dans le but de les rendre opérationnelles à ces échelles spatio-temporelles. Cette démarche vise, comme évoqué précédemment, à mieux comprendre l'impact du changement climatique dans des régions où parfois peu de variables climatiques sont mesurées, mais aussi à établir des jeux de données valides et robustes et de les rendre disponibles à la communauté scientifique pour des utilisations diverses (glaciologiques, hydrologiques, climatiques…). Ces motivations ainsi que les principaux résultats obtenus lors de mes trois années de thèse sont regroupés dans cette Section 3. Le travail de recherche effectué a mené à l'écriture de trois articles scientifiques ainsi qu'à la participation à un quatrième article de synthèse, tous ayant comme ligne directrice, le développement et l'application de méthodes permettant la quantification du SMB annuel et saisonnier à l'échelle régionale, à partir d'images satellites optiques. Pour chaque étude, un résumé ainsi que le ou les articles associés (en anglais, du fait de leur publication dans des journaux à visée internationale) est présenté. Cette Section se divise ainsi en trois sous-parties : 1. Synthèse des méthodes indirectes permettant la quantification du SMB annuel et saisonnier des glaciers de montagne à partir d'images satellites optiques (1 article) 2. Développement d'un algorithme permettant l'identification automatique de la SLA sur des glaciers de montagne, à partir de différents capteurs satellitaires optiques (1 article) 3. Estimation du SMB annuel et saisonnier à partir d'images satellites optiques : application sur les Alpes européennes et les Andes tropicales (3 articles) 3.1 SYNTHESE DES METHODES INDIRECTES PERMETTANT LA QUANTIFICATION DU SMB ANNUEL ET SAISONNIER DES GLACIERS DE MONTAGNE A PARTIR D'IMAGES SATELLITES OPTIQUES Article 1 : Comme illustré dans la Sect. 2 de ce manuscrit, il existe différentes méthodes permettant d'estimer le SMB annuel et saisonnier de glaciers de montagne à partir d'images satellites. En revanche, peu d'articles de synthèse existent dans la littérature (e.g. Bamber and Rivera, 2007; Cuffey and Paterson, 2010; Rees, 2005), et sont parfois trop anciens pour présenter des méthodes développées récemment. Dans cet article de synthèse, trois méthodes permettant l'estimation du SMB annuel et/ou saisonnier à partir d'images satellites optiques sont présentées. Pour chaque méthode présentée, le contexte historique est rappelé ainsi que son principe et sa mise en place. Un exemple d'application de chaque méthode illustre ensuite leurs potentiels et permet de soulever certaines limitations. Puis, une section discute des limites communes identifiées, comme la présence de nuages non-détectés, véritable verrou de la télédétection spatiale optique, ou encore la pertinence d'appliquer ces méthodes sur des glaciers couverts de débris. Les limites spécifiques à chaque méthode sont aussi énoncées, notamment en lien avec la résolution spatiale des différents capteurs utilisés (e.g. taille de glacier minimale requise) où les caractéristiques des glaciers étudiés (e.g. altitude maximale, régime de précipitations). Finalement, une dernière partie évoque les perspectives d'améliorations. J'ai participé au développement de deux de ces trois méthodes, ainsi qu'à la relecture et à la modification du manuscrit. Ma contribution à la méthode Résumé de mes recherches 47 albédo a fait l'objet d'une publication indépendamment de cet article de synthèse (Sect. 7, article 4). De par l'état d'avancement de mes travaux au moment de la publication (2017), ma contribution à la méthode ELA n'est pas prise en compte dans cette étude mais la motivation de mes recherches sur cette méthode découle des perspectives évoquées, notamment quant à la nécessité d'automatiser cette méthode et de l'appliquer à grande échelle.

  le lien entre le SMB et la pente moyenne, les altitudes médiane et maximale du glacier, expliquant respectivement 19, 24 et 12% de la variance observée. Notre étude suggère ainsi que les glaciers fortement pentus et hauts en altitude présentent des pertes de masse moins négatives. Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data Lucas Davaze, Antoine Rabatel, Yves Arnaud, Pascal Sirguey, Delphine Six, Anne Letreguilly, and Marie Dumont. The Cryosphere 12, 271-286 (2018). 2L'albédo est représentatif de l'état de surface du glacier et détermine la quantité de rayonnement solaire absorbée par la surface. Grâce à l'apparition de capteurs optiques multi-spectraux ayant un temps de revisite quasi-Annual surface mass balances quantification for 82 glaciers in the tropical Andes over a 18-yr period from optical remote sensing data Léna Gurriaran, Lucas Davaze, Inès Dussaillant, Antoine Rabatel, Etienne Berthier, Yves Arnaud 3

	Résumé de mes recherches
	Article 4 :
	53

Dans l'échantillon de SMB issu de données in situ (23 glaciers), aucun glacier n'a une pente moyenne supérieure à 21.6°. Notre étude révèle en revanche une forte dépendance de la perte de masse cumulée à la pente moyenne du glacier, où les glaciers ayant des pentes supérieures à 20.5° (correspondant à 45% de la surface englacées des Alpes européennes) subissent une perte de masse ~25% moindre que les glaciers moins pentus, sur la période 2000-2016. Nous suggérons ainsi que les mesures in situ des 23 glaciers utilisées dans cette étude ne sont pas représentatives des bilans calculés à l'échelle des Alpes européennes. En conclusion, les estimations du bilan de masse régionale basées sur une extrapolation de ces mesures doivent nécessairement prendre en compte les caractères morpho-topographiques des glaciers lors de l'extrapolation.

Nous avons ensuite analysé la variabilité spatio-temporelle des SMBs annuels. Une analyse statistique nous a permis d'identifier deux principales régions, ayant des variabilités interannuelles similaires en termes de SMB et de forçages climatiques (données issues de réanalyses ERA5). La première région correspond à des glaciers situés au Nord-Est de la chaîne, semblant être fortement impactés par la quantité de précipitations (eau précipitable sur toute la colonne atmosphérique expliquant de 31 à 37% de la variance de SMB observée). Cette observation confirme une précédente étude

[START_REF] Abermann | Climatic controls of glacier distribution and glacier changes in Austria[END_REF]

suggérant que les glaciers dans cette région existent à de telles altitudes du fait d'importantes précipitations. La deuxième région correspond à la partie occidentale de la chaîne, où les variabilités observées de SMB semblent davantage expliquées par l'advection d'humidité provenant de l'Ouest, ainsi qu'à des conditions météorologiques estivales, dominées par des flux allant vers le sud (45 à 53% de variance expliquée). Au-delà des résultats, cette étude propose un cadre méthodologique, contraint sur une région où de nombreuses données de validation existent, permettant d'estimer le bilan pluri-décennal et annuel à l'échelle régionale pour des glaciers individuels, à partir d'images satellites optiques. Nous espérons que cette étude ouvrira la porte à d'autres applications sur des régions où les enjeux liés aux glaciers de montagne, est plus importante que dans les Alpes Européennes (e.g., hautes montagnes d'Asie, Andes).

journalier (e.g. MODIS), des méthodes ont été développées afin de déterminer le SMB annuel et saisonnier de glaciers de montagne à partir de l'albédo de surface glaciaire

(détails Sect. 2.3 et 4.2.2.)

. Les récentes avancées en termes de puissance de calcul et de développement d'algorithmes permettant le traitement d'un grand nombre d'images permettent aujourd'hui d'appliquer ces méthodes à l'échelle régionale. En utilisant la méthode présentée par

Sirguey et al. (2016)

, des cartes d'albédo ont été produites à l'échelle de 30 glaciers dans les Alpes françaises entre 2000 et 2016. L'analyse de ces cartes a permis de calculer un albédo moyen à l'échelle du glacier, marqueur des processus de surface

(détails Sect. 2.3 et 4.2.2.)

. Grâce à des données de SMB de validation, nous avons pu confirmer la capacité de l'albédo moyen minimum à être un indicateur du SMB annuel pour 27 des 30 glaciers considérés, validant ainsi l'utilisation de cette méthode à une échelle régionale.

A l'échelle saisonnière et pour six glaciers pour lesquels des données de SMB saisonniers existent, nos résultats montrent que le suivi de l'évolution temporelle de l'albédo moyen sur la période estivale est un indicateur du SMB estival et une application régionale est aujourd'hui réalisable. A contrario, l'analyse de l'évolution temporelle de l'albédo moyen pendant l'hiver permet de quantifier la fréquence des précipitations mais pas la quantité. Ainsi, notre étude sur plusieurs glaciers révèle que la méthodologie développée par Article 5 :

2

  National School of Surveying, Univ. Otago, Dunedin, New Zealand; pascal.sirguey@otago.ac.nz

									et
	océanographie	spatiales	(LEGOS),	F-31000	Toulouse,	France;
	vanessa.drolon@gmail.com;		philippe.maisongrande@cnes.fr;
	etienne.berthier@legos.obs-mip.fr
	4	Météo-France, CNRS, Centre national de recherches météorologiques / Centre
	d'études	de	la	neige	(CNRM/CEN),	F-38000	Grenoble,	France;
	marie.dumont@meteo.fr

3

Univ. Toulouse, CNES, CNRS, IRD, Laboratoire d'études en géophysiques

Table S1 .

 S1 Sensor specifications. λNIR and λSWIR stands for wavelength of near and shortwave infrared radiations, respectively. *Res stands for the spatial resolution of the image for the specified band. **METI stands for the Japan Ministry of Economy, Trade and Industry.

	Sensor	λNIR [μm]	Res* [m]	λSWIR [μm]	Res* [m]	Tile size km] [km *	Revisit [days] time	Time span	Data provider	Processing level	Failure
	Landsat 5	0.76-0.90	30	1.55-1.75	30	170*185 16	1984-2013	NASA	L1T
											Scan line
	Landsat 7	0.77-0.90	30	1.55-1.75	30	170*185 16	1999-… NASA	L1T	corrector failure
											May 2003
	Landsat 8	0.85-0.88	30	1.57-1.65	30	170*185 16	2013-… NASA	L1T
											SWIR
	ASTER	0.76-0.86	15	1.60-1.70	30	60*60	16	2000-April 2008	NASA/ METI**	L1T	detectors failure Apr.
											2008
	Sentinel 2A&B	0.76-0.91	10	1.54-1.68	20	100*100 5	2015-… ESA	L1C

Table 1 :

 1 Glacier datasets used in the present study, together with their respective cumulative area and the corresponding fraction of the RGI 6.0.

	Glacier dataset	Glacier count [Nb]	Area [km2]	Fraction of the RGI 6.0 [%]
	RGI 6.0	3,927	2,092	100
	Geodetic	361	1,563	75
	SLA-method	239	1,215	58
	6.3 DATA & METHODS

Table 2 :

 2 Results from multi-linear LASSO regression analysis in between median annual SMB and nine climatic variables from ERA5. r 2 b3 displays the correlation coefficient with the computed inter-annual SMB considering the most explicative climatic variables (three at maximum), and r 2 b, the correlation coefficient of the most explicative variable. Scores displayed in bold are significant (risk<5%) according to a Student t-test. Most explicative variables are classified function of their relative significativity, with red the most explicative, salmon, the second and blue the third.As described in Sect 6.4.2, computed annual SMB with the SLA-method underestimates extremely high and low SMB values (i.e. 2001, the only year with positive SMB from the studied period and 2003, a very low SMB year in the European Alps). As our studying period has both extremes in the first four years, the cumulative SMB is at the beginning poorly reproduced by remotelysensed estimates but is then able to match in situ monitored SMB for the rest of the period. The glacier mean slope of the in situ monitored glaciers falls into class 1 and 2 (slopes gentler than 20.5°). The median of the 23 in situ SMB time series is in agreement with the remote-sensing SMB for glaciers having a similar mean slope (inset Tablein Figure

	Region-wide annual glacier surface mass balance for the European Alps from
	2000 to 2016

  for the correction of shadows. It is first computed at 125m resolution, providing Booleantype products of self and cast shadows per pixel. Results are then averaged and aggregated to 250m resolution, producing a sub-pixel fraction of shadow (further detailed inSirguey et al.,

	7.3.3 Surface albedo and glacier-wide mass balance relationship
	Basis of the method
	2009). Finally, MODIS data processed with MODImLab provide, among others,
	near-daily maps of white-sky albedo at 250m resolution together with cloud
	masks and cast and projected shadows.
	Albedo maps have been processed for 5068 images for the Ecrins range,
	4973 for Mont Blanc and 5082 for Vanoise over the period 2000-2015. Only
	images acquired between 09:50 and 11:10UTC (+2h in summer for local time
	conversion) were selected to get minimum SZA and limit projected shadows
	of surrounding reliefs.
	7.3.2 Glacier masks
	Following Dumont et al. (2012) and Brun et al. (2015), we manually created
	raster masks of the 30 glaciers, based on the glaciers' outlines from the 1985-
	1987 (Rabatel et al., 2013) and high spatial resolution (6m) SPOT-6 images
	from 2014. All debris-covered areas, together with mixed pixels (rocksnow/ice)
	have been removed to capture only the snow/ice albedo signal. The resulting
	number of pixels per glacier is listed in Table 1.

Table 2 .

 2 Filtering the images from OZA values.

  ). Conversely to

	Monitoring glacier albedo as a proxy to derive summer and annual surface
	mass balances from optical remote-sensing data
	method but not remotely sensed. To conclude,
		min a	, the	int s  is slightly more
	robust to the presence of undetected clouds as its value does not rely on a
	single image. The lowest correlation has been found for Talèfre Glacier. The
	latter accounts for a relatively large debris-covered tongue that has been
	excluded when delineating the glacier mask. Consequently, the low correlation
	could be partly explained by this missing area, considered in the glaciological
	175

int s

Table 3

 3

	gathers all the

  . La variable morpho-topographique expliquant le plus la variance des bilans cumulés observés semble varier d'une région à l'autre : pente moyenne, altitudes médiane et maximale dans les Alpes européennes ; altitudes minimale et maximale en Cordillère Blanche ; altitude maximale en Cordillère Royale). En règle générale, les glaciers situés plus haut en altitude, et les plus pentus (dans les Alpes européennes) subissent des pertes de masse moins importantes. Finalement, l'impact des variables morpho-topographiques et climatiques est difficilement transposable d'une région à l'autre et des études ciblées sur des régions englacées sont désormais nécessaires, pour estimer le SMB pluriannuel (i.e. projet ASTERIX) et annuel (i.e. méthode ELA). Les mesures in situ d'accumulation et d'ablation constituent la mesure de référence pour estimer le bilan de masse de surface, notamment grâce à la méthode glaciologique. Ces mesures permettent aussi d'améliorer notre connaissance des processus de surface, notamment en lien avec le bilan d'énergie de surface lorsqu'elles sont couplées à l'installation d'une station météorologique (e.g.

Ce travail a mené à une présentation à l'AGM, à l'AGU en

ainsi qu'à l'IGE lors de la venue des thématiciens du CNES en Janvier 2018.

Grâce au développement de l'algorithme permettant d'estimer de manière semi-automatique l'altitude de la SLA, nous avons pu envisager une application de la méthode ELA à une échelle régionale. Poursuivant le travail effectué sur 44 glaciers des Alpes occidentales par Rabatel et al. (2016), les SLAs de 239 glaciers ont été estimées à l'échelle des Alpes européennes. Bénéficiant d'une estimation du SMB pluri-décennal issu de la méthode géodésique, et produite par Romain Hugonnet et Etienne Berthier dans le cadre du projet ASTERix, ces SLAs ont pu être utilisées pour en déduire le SMB annuel de ces glaciers sur la période 2000-2016. Cette approche étant novatrice à une échelle régionale, nous avons tout d'abord validé les SMB annuels reconstruits grâce à des SMB issus de données in situ de 23 glaciers, répartis sur les Alpes européennes. Les estimations in situ et télédétectées sont significativement corrélé (r 2 =0.33, n=282) et la variabilité interannuelle médiane reconstruite est capable de reproduire celle estimée à partir de données in situ (r 2 =0.45, n=17). Les différences observées proviennent principalement des lacunes de notre algorithme de détection à estimer les SLAs pour les années exceptionnellement déficitaires, ou à bilan positif, du fait de l'emplacement de la SLA, haute dans la zone d'accumulation, ou proche du front, complexifiant son identification (discuté en détail dans l'article 3, Sect. 5.4). Une autre source de différence entre les deux jeux de données provient du choix des glaciers pour le suivi in situ et de l'extrapolation

Ce travail a mené à une présentation à l'EGU en 2017(Davaze et al., 2017conf) ainsi qu'à l'IGE lors de la venue des thématiciens du CNES en janvier 2017.

Sirguey et al. (2016) ne peut pas être utilisée pour déterminer le SMB hivernal sur les glaciers des Alpes françaises. Finalement, une analyse détaillée des masques de nuage produits par la chaîne de traitement utilisée confirme la robustesse de l'algorithme utilisé et le faible d'impact du seuil du masque des nuages sur nos résultats. Une analyse similaire révèle en revanche des biais systématiques issus des corrections topographiques et des améliorations sont requises si les cartes d'albédo sont utilisées hors d'une comparaison interannuelle.

Ce travail a mené à une présentation faite par Antoine Rabatel à l'IUGG en 2019[START_REF] Gurriaran | Annual glacier-wide mass balances of 82 glaciers in the tropical Andes over a 18-yr period from optical remote sensing data[END_REF].

termes de transposabilité dans une autre région, du fait des nombreux « crashtest » qu'elle a effectué. Dans une région où peu de données existent et où les glaciers ont une importance capitale dans les cycles hydrologiques locaux, ce jeu de données de SMB annuel, lorsqu'il sera publié, apportera une vraie plusvalue à la communauté scientifique internationale et locale, notamment intéressée par l'utilisation de l'algorithme de détection de la ligne de neige pour des inventaires nationaux (collaboration avec Katy Medina du National Glacier and Mountain Ecosystem Research Institute (INAIGEM)). Le rapport de de fin de stage de Léna tel qu'elle l'a écrit est présenté dans l'annexe 1.Figure15: SMB pluriannuels issus de la méthode géodésique sur la période 2000-2018 en fonction de six variables morpho-topographiques. Les coefficients de corrélation sont indiqués pour chaque comparaison. La ligne rouge représente la régression linéaire robuste calculée à partir de la méthode M-estimateur[START_REF] Susanti | [END_REF].

Acknowledgments: The authors are grateful to the spatial agencies: CNES, ESA, NASA, VITO, for providing free access (or with a reduced cost) to the satellites images. The French Service National d'Observation GLACIOCLIM (http://glacioclim.osug.fr/) and the World Glacier Monitoring Service (http://wgms.ch/) are acknowledged for the collection and distribution of in situ glaciological data used to calibrate/validate the remote-sensing based method. Authors from the University of Grenoble Alpes acknowledge the support of LabEx OSUG@2020 (Investissements d'avenir -ANR10 LABX56) and of EquipEx GEOSUD (Investissements d'Avenir ANR-10-EQPX-20). P. Sirguey was supported by a University of Otago Research Grant (ORG12-10864401RFW) and support of the Grenoble Institute of Technology (Grenoble-INP) and the University Grenoble-Alpes (UGA). The MODIS Level-1B data were processed by the MODIS Adaptive Processing System (MODAPS) and the Goddard Distributed Active Archive Center (DAAC), and are archived and distributed by the Goddard DAAC. V. Drolon and P. Maisongrande thank Matthias Huss for providing the seasonal SMB of 48 European Alpine glaciers. They are also grateful to VITO and Belspo for the SPOT/VGT satellite images distribution and for their scientific contribution to carry out this study. V. Drolon PhD was supported by the CNES and VITO/CLS. Finally, we acknowledge Aubrey Miller for proofreading the manuscript, as well as Dr. Frank Paul (Editor) and the three anonymous reviewers for their constructive comments used to improve the paper.

Acknowledgments

This study was conducted in the context of the French glacier observatory GLACIOCLIM (https://glacioclim.osug.fr) and the Labex OSUG@2020 (Investissements d'Avenir -ANR10 LABX56). We are grateful to the satellite data providers: Copernicus/EU/ESA for Sentinel-2, USGS for Landsat, and JPL/NASA and METI for ASTER. The authors also acknowledge Fabien Maussion for processing the glacier central flow lines, Etienne Berthier and Delphine Six for fruitful discussions on the method and the manuscript, Juliette Blanchet for advises on the statistical approach developed for this algorithm, Maxim Lamare for discussions on topographic correction models and Fanny Brun for glaciological discussions. Acknowledgements. This study was conducted within the Service National d'Observation GLACIOCLIM and supported by the French National Research Institute for Sustainable Development (IRD) and the French National Research Agency (ANR) through ANR-13-SENV-0005-04-PRESHINE. Equipex GEOSUD (Investissements d'avenir -ANR-10-EQPX-20) is acknowledged for providing the 2014 SPOT-6 images. The MODIS Level-1B data were processed by the MODIS Adaptive Processing System (MODAPS) and the Goddard Distributed Active Archive Center (DAAC) and are archived and distributed by the Goddard DAAC. In situ mass balance data for the Glacier Blanc were kindly provided by the Parc National des Ecrins. The authors acknowledge the contribution the Labex OSUG@2020 (Investissements d'avenir -ANR10 LABX56). Pascal Sirguey thanks the University of Grenoble Alpes and Grenoble-INP for the 6-month "invited professor grant" obtained in 2015-2016.

Data availability. Processing of the albedo images has been performed using the opensource MODImLab algorithm. This algorithm can be accessed by contacting its administrator, P. Sirguey. MODIS images have been downloaded using the permanent ftp server: n5eil01u.ecs.nsidc.org maintained by the National Aeronautics and Space Administration (NASA) Distributed Active Archive Center (DAAC).

Author Contributions: Antoine Rabatel coordinated the writing, Pascal Sirguey, Vanessa Drolon and Philippe Maisongrande contributed to the writing of the manuscript. Yves Arnaud, Lucas Davaze, Marie Dumont, Etienne Berthier and Jean-Pierre Dedieu have contributed to the development of the methods presented in this review paper and revised the manuscript.

Conflicts of Interest: The authors declare no conflict of interest.

obstruction by the surrounding topography from the input DEM at 30 m resolution. The used algorithm is fully described by Sirguey (2009) and was inspired by Dozier and others (1981) and Dozier and Frew (1990). Self shadows are merged to the cast shadow product derived from the illumination angle. An example of detected shadows on a Landsat image is illustrated Fig. 3D. Quality assessment of the produced shadow masks has been encompassed in Davaze and others (2018) and shaded areas have then been masked. Fig. S4. Example of different band composites for the Blanc Glacier in the French Alps on a Landsat 8 image from the 8 Sept. 2016. (A) corresponds to a RGB composite from the band 6, 5 and 3 (green, NIR, SWIR) as recommended by Rabatel and others (2012), (B) to a band ratio 4/6 (Red/SWIR) as recommended by Winsvold and others (2016), (C) to a band 5/6 (NIR/SWIR), and (D) to an enhanced band ratio 5²/6. The red glacier outline corresponding to 1985 is illustrative.

A semi-automated method to estimate the glacier end-of-summer snowline altitude from optical satellite images 109 Visual verification of eligible images Only relying on the NDSI index to detect high altitude and snow clouds is challenging, especially using the rather low NDSI threshold we chose. To mitigate the impact of undetected clouds, a visual verification of the images is currently performed. It also allows to discard images when the glacier is completely covered with fresh snow that generates erroneous detection of SLA.

For convenience, a RGB composite of the green, NIR and SWIR band as recommended by Rabatel andothers (2012, 2016) is produced for each image (as displayed Fig. S3E). At this stage, glaciers with particularly patchy snowlines are discarded to avoid detection errors.

Definition of a snowline detection footprint

Because of outline uncertainties, compounded with georectification errors, contrasting targets at the edge of glacier areas can lead to mapping errors.

Therefore, only a central buffer along the glacier flow line, far away enough from the margin, is considered for snowline detection (example illustrated Fig. S3F). This footprint extends from 50 m above the glacier terminus to 100 m below its top. The width of this footprint is located along the main glacier flow line and is defined as 40% of the mean width of the glacier and is computed as follow:

Where WB is the width of the buffer, A the glacier area, and L its maximum length.

Unsupervised classification

An unsupervised K-Means classification is performed on the resulting processed band ratio within the snowline detection footprint, following the algorithm proposed by [START_REF] Lloyd | Least squares quantization in PCM[END_REF], with two classes (ice = 1, snow = 2).

CONCLUSION

We quantified a spatially resolved estimate of individual annual SMB time series for 239 glaciers in the European Alps, using both geodetic and snowline altitude estimates for the period 2000-2016. Benefiting from in situ SMB data from 23 glaciers, we could validate our approach, based on the SLA-method to quantify the SMB of individual glaciers at regional scale. The observed discrepancies with in situ estimates have been assigned to both the uncertainties in the SLA using a semi-automatic algorithm (described in the supplementary material) and to glacier-wide estimates from in situ surface mass balances.

Geodetic estimates of long term mass changes were used to analyse the impact of glacier morpho-topographic features on the mass changes. Three morpho-topographic variables: mean glacier slope, median and maximum altitude, appear to explain 19, 24 and 12%, respectively, of the observed variance of the long term mass changes. An analysis of the cumulative SMB from in situ and remote-sensing estimates, as a function of the glacier slopes, suggests that steeper glaciers experienced less mass losses. Our analysis suggests that SMB time series quantified from the glaciers with in situ measurements are not representative of the regional behaviour of glaciers in the European Alps. Indeed, these glaciers belong to the "gentle slope glaciers" (< 20 °), and half of the glaciers in the European Alps, with higher slope experimented as lower mass loss. Nonetheless, in situ measurements still provide essential data for the comprehension of surface processes and glacier response to climate forcing at fine scales.

Analysis of the temporal SMB variability in combination with climate forcing from ERA5 reanalysis data, allowed identifying two main sub-regions within the European Alps. In the first one at the north-eastern part of the European Alps, the temporal SMB variability seems to be mainly driven by the precipitable water of the entire atmospheric column, suggesting that, in agreement with previous studies (e.g. (GLAMOS, 2018). Data from glaciers from Italy and Austria were downloaded on the WGMS database (WGMS, 2017). We kindly acknowledge all the data providers for their work on the respective glaciers. J'ai ainsi décidé de l'expérimenter sur les expériences que j'avais initiées avec des élèves.

Réflexion sur le sens profond et la préparation de mes enseignements

Le sens profond caractérise le fait d'avoir compris dans toute sa complexité un concept, un théorème ou une conjecture. Si l'on a compris le sens profond d'une loi en physique par exemple, on ne peut plus l'oublier, ou se tromper dans les termes de cette loi. J'ai découvert que je n'avais pas le sens profond de certains concepts que j'enseignais, notamment en physique, alors que j'ai eu, au long de ma scolarité des résultats en physique honorables. Je trouve la notion de sens profond étroitement liée avec celle d'obstacle épistémologique, décrivant quelque chose de complexe à comprendre pour l'apprenant en mobilisant uniquement les concepts qu'il maitrise déjà et dont il a acquis le sens profond. Je me suis donc questionné sur ces obstacles épistémologiques.

Où sont-ils ? Quel concept est difficile à appréhender pour mon auditoire ? Je me suis rendu compte, notamment en discutant avec le directeur de l'unité d'enseignement dans laquelle j'ai enseigné ma première année de doctorat, et profitant de son expérience, que certains concepts étaient "difficiles". A savoir,  Les élèves butent sur sa compréhension chaque année Annexe 3 : "Patauge pédagogique " 241  Les élèves ont l'air de le comprendre, mais parviennent rarement à réutiliser et s'approprier cette notion par la suite Lors de la préparation du cours, j'ai ressenti deux façons d'identifier ces concepts difficiles/obstacles épistémologiques. La première est que je me suis rendu compte parfois de ne pas avoir réellement compris ce que j'enseignais, ce que j'ai vécu comme une pure trahison envers les élèves, lorsque je m'en rendais compte (e.g. "Prenez cela comme acquis" ou je vais présenter cela rapidement pour ne pas avoir de questions qui me piégeraient). La deuxième est que j'avais au début une démarche où je me basais sur mes propres expériences d'apprenant dans mes préparations de cours. Si un concept m'avait paru clair et limpide lorsque j'étais élève, alors je ne prévoyais pas de temps pour l'expliquer en détail. En revanche, je concentrais mon enseignement sur des points qui m'avait paru difficiles.

Au cours de nombreuses discussions avec d'autres doctorants/enseignants, il est commun de préparer les cours de cette façon (e.g. je ne comprends pas que les élèves ne comprennent pas les fonctions (en informatique) alors que pour moi cela a été très simple). Après quelques réflexions, je suis convaincu qu'en tant qu'enseignant nous ne sommes pas comme tous les élèves. En Le fait de changer l'énoncé n'a pas forcément détourné les élèves de ce problème de glace "bleue". En réalité ceci n'est pas forcément un problème car cela me permet d'introduire les couleurs, même si cela intervient un peu tôt pour des élèves de 5ème/4ème. Lors de mes cinq essais avec ce nouvel énoncé, le débat a été plus orienté sur des problématiques glaciologiques que des problèmes de chien et de jardin, ce dont je me réjouis. En revanche, beaucoup d'élèves se préoccupent plus du problème la qualité de l'eau destinée au chien (e.g. proposition de faire fondre de la neige ou de la neige sale en fonction des saletés qui peuvent altérer la potabilité de l'eau...). En effet, j'ai par la suite fait attention à employer des termes plus "scolaire" et plus rigoureux.