# Nonlinear sensitivity of glacier mass balance to future climate change unveiled by deep learning

Jordi Bolibar, Antoine Rabatel, Isabelle Gouttevin, Harry Zekollari, Clovis Galiez

► To cite this version:

Jordi Bolibar, Antoine Rabatel, Isabelle Gouttevin, Harry Zekollari, Clovis Galiez. Nonlinear sensitivity of glacier mass balance to future climate change unveiled by deep learning. Nature Communications, 2022, 13 (1), pp.409. 10.1038/s41467-022-28033-0. insu-03668394

HAL Id: insu-03668394

https://insu.hal.science/insu-03668394

Submitted on 14 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ARTICLE

https://doi.org/10.1038/s41467-022-28033-0

OPEN

<figure><img src="image_3.png" /></figure>

Check for updates

# Nonlinear sensitivity of glacier mass balance to future climate change unveiled by deep learning

Jordi Bolibar ID <fnref n="1" />,²,<fnref n="3" />✉, Antoine Rabatel ID <fnref n="1" />, Isabelle Gouttevin ⁴, Harry Zekollari ⁵,⁶ & Clovis Galiez ⁷

Glaciers and ice caps are experiencing strong mass losses worldwide, challenging water availability, hydropower generation, and ecosystems. Here, we perform the first-ever glacier evolution projections based on deep learning by modelling the 21st century glacier evolution in the French Alps. By the end of the century, we predict a glacier volume loss between 75 and 88%. Deep learning captures a nonlinear response of glaciers to air temperature and precipitation, improving the representation of extreme mass balance rates compared to linear statistical and temperature-index models. Our results confirm an over-sensitivity of temperature-index models, often used by large-scale studies, to future warming. We argue that such models can be suitable for steep mountain glaciers. However, glacier projections under low-emission scenarios and the behaviour of flatter glaciers and ice caps are likely to be biased by mass balance models with linear sensitivities, introducing long-term biases in sea-level rise and water resources projections.

<footnote n="1">¹Univ. Grenoble Alpes, CNRS, IRD, G-INP, Institut des Géosciences de l'Environnement, Grenoble, France. ²INRAE, UR RiverLy, Lyon-Villeurbanne, France.</footnote>

<footnote n="3">³Institute for Marine and Atmospheric research Utrecht, Utrecht University, Utrecht, Netherlands. ⁴Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM, Centre d'Études de la Neige, Grenoble, France. ⁵Department of Geoscience and Remote Sensing, Delft University of Technology, Delft, Netherlands. ⁶Laboratoire de Glaciologie, Université Libre de Bruxelles, Brussels, Belgium. ⁷Univ. Grenoble Alpes, CNRS, G-INP, Laboratoire Jean Kuntzmann, Grenoble, France. ✉email: j.bolibar@uu.nl</footnote>
Glaciers are experiencing important changes throughout the world as a consequence of anthropogenic climate change¹. Despite marked differences among regions, the generalized retreat of glaciers is expected to have major environmental and social impacts²,³. Water resources provided by glaciers sustain around 10% of the world's population living near mountains and the contiguous plains⁴, depending on them for agriculture, hydropower generation⁵, industry or domestic use. Several aquatic and terrestrial ecosystems depend on these water resources as well, which ensure a base runoff during the warmest or driest months of the year⁶. Predicting future glacier evolution is of paramount importance in order to correctly anticipate and mitigate the resulting environmental and social impacts. During the last decade, various global glacier evolution models have been used to provide estimates on the future sea-level contribution from glaciers⁷,⁸. All these glacier models, independently from their approach, need to resolve the two main processes that determine glacier evolution: (1) glacier mass balance, as the difference between the mass gained via accumulation (e.g. snowfall, avalanches and refreezing) and the mass lost via different processes of ablation (e.g. melt and sublimation of ice, firn and snow; or calving)⁹; and (2) ice flow dynamics, characterized by the downward movement of ice due to the effects of gravity in the form of deformation of ice and basal sliding. Simulating these processes at a large geographical scale is challenging, with models requiring several parametrizations and simplifications to operate. Recent efforts have been made to improve the representation of ice flow dynamics in these models, replacing empirical parametrizations with simplified physical models⁹,¹⁰. Nonetheless, to represent the glacier mass balance, the vast majority of large-scale glacier evolution models relies on temperature-index models. This type of model uses a calibrated linear relationship between positive degree-days (PDDs) and the melt of ice or snow¹¹. The main reason for their success comes from their suitability to large-scale studies with a low density of observations, in some cases displaying an even better performance than more complex models¹². However, both the climate and glacier systems are known to react non-linearly, even to pre-processed forcings like PDDs¹³, implying that these models can only offer a linearized approximation of climate-glacier relationships. Deep artificial neural networks (ANNs) are nonlinear models that offer an alternative approach to these classic methods. However, the use of ANNs remains largely unexplored in glaciology for regression problems, with only a few studies using shallow ANNs for predicting the ice thickness¹⁴ or mass balance¹³ of a single glacier.

The French Alps, located in the westernmost part of the European Alps, experience some of the strongest glacier retreat in the world¹⁵⁻¹⁷. Long-term historical interactions between French society and glaciers have developed a dependency of society on them for water resources, agriculture, tourism¹⁸—particularly the ski business¹⁹—and hydropower generation. This rapid glacier retreat is already having an environmental impact on natural hazards²⁰, mountain ecosystems²¹, and biodiversity⁶. Without these cold water resources during the hottest months of the year, many aquatic and terrestrial ecosystems will be impacted due to changes in runoff, water temperature or habitat humidity⁶,²¹,²². An accurate prediction of future glacier evolution will be crucial to successfully adapt socioeconomic models and preserve biodiversity. Glaciers in the European Alps have been monitored for several decades, resulting in the longest observational series in the world²³,²⁴. In this study, we investigate the future evolution of glaciers in the French Alps and their nonlinear response to multiple climate scenarios. We perform, to the best of our knowledge, the first-ever deep learning (i.e. deep artificial neural networks) glacier evolution projections by modelling the regional evolution of French alpine glaciers through the 21st century. For

this, a newly-developed state-of-the-art modelling framework based on a deep learning mass balance component and glacier-specific parametrizations of glacier changes is used. The climatic forcing comes from high-resolution climate ensemble projections from 29 combinations of global climate models (GCMs) and regional climate models (RCMs) adjusted for mountain regions for three Representative Concentration Pathway (RCP) scenarios: 2.6, 4.5, and 8.5²⁵. The 29 RCP-GCM-RCM combinations available, hereafter named climate members, are representative of future climate trajectories with different concentration levels of greenhouse gases (Table S1). With this study, we provide new predictions of glacier evolution in a highly populated mountain region, while investigating the role of nonlinearities in the response of glaciers to multiple future climate forcings.

## Results

**Glacier evolution through the 21st century.** Our projections show a strong glacier mass loss for all 29 climate members, with average ice volume losses by the end of the century of 75%, 80%, and 88% compared to 2015 under RCP 2.6 (±9%, *n* = 3), RCP 4.5 (−17% +11%, *n* = 13) and RCP 8.5 (−15% +11%, *n* = 13), respectively (Fig. 1 and S1). Differences in projected glacier changes become more pronounced from the second half of the century, when about half of the initial 2015 ice volume has already been lost independent of the considered scenario. Annual glacier-wide mass balance (MB) is estimated to remain stable at around −1.2 m.w.e. *a*⁻¹ throughout the whole century under RCP 4.5, with glacier retreat to higher elevations (positive effect on MB) compensating for the warmer climate (negative effect on MB). Conversely, for RCP 8.5, annual glacier-wide MB are estimated to become increasingly negative by the second half of the century, with average MB almost twice as negative as today's average values (Fig. 1a). MB rates only begin to approach equilibrium towards the end of the century under RCP 2.6, for which glaciers could potentially stabilize with the climate in the first decades of the 22nd century depending on their response time (Fig. S1a). An analysis of the climate signal at the glaciers' mean altitude throughout the century reveals that air temperature, particularly in summer, is expected to be the main driver of glacier mass change in the region (Fig. 1). Interestingly, future warmer temperatures do not affect annual snowfall rates on glaciers as a result of both higher precipitation rates in the EURO-CORDEX ensemble (Fig. 1g-i)²⁶ and glaciers shrinking to higher elevations where precipitation rates are higher as a result of orographic precipitation enhancement²⁷. The increase in glacier altitude also causes the solid to liquid precipitation ratio to remain relatively constant. Therefore, solid precipitation is projected to remain almost constant at the evolving glaciers' mean altitude independently from the future climate scenarios, while air temperature is projected to drive future glacier-wide mass changes (Fig. 1d, g). These results are in agreement with the main known drivers of glacier mass change in the French Alps²⁸. Overall, the evolving glaciers are expected to undergo rather stable climate conditions under RCP 4.5, but increasingly higher temperatures and rainfall under RCP 8.5 (Fig. 1). These differences in the received climate signal are explained by the retreat of glaciers to higher altitudes, which keep up with the warming climate in RCP 4.5 but are outpaced by it under RCP 8.5.

The glacier ice volume in the French Alps at the beginning of the 21st century is unevenly distributed, with the Mont-Blanc massif accounting for about 60% of the total ice volume in the year 2015 (7.06 out of 11.64 km³, Fig. 2a). The vast majority of glaciers in the French Alps are very small glaciers (<0.01 km²), that are mainly remnants from the Little Ice Age, with a strong imbalance with the current climate¹⁵. Our projections highlight the almost complete disappearance of all glaciers outside the Mont-Blanc and Pelvoux
<figure><img src="image_4.png" /><figcaption>**Fig. 1** 21st century glacier-wide MB, geometry and climate signal evolution of French Alpine glaciers. **a** Glacier-wide annual MB, **b** Ice volume, **c** Glacier area. The cumulative positive degree days (CPDD), snowfall and rainfall **d–I**, are at the glaciers' annually evolving centroids. As such, these values reflect both the climatic forcing and the changing glacier geometry. Summer climate is computed between April 1st and September 30th and winter climate between October 1st and March 31st. Thin lines represent each of the 29 individual member runs, while the thick lines represent the average for a given RCP.</figcaption></figure>

(Ecrins region) massifs under RCP 4.5 (Fig. 2) and RCP 8.5 by the end of the century. By 2100, under RCP 4.5, these two high-altitude massifs are predicted to retain on average 26% and 13% of their 2015 volume, respectively, with most of the ice concentrated in a few larger glaciers (>1 km², Fig. S4). Glacier landscapes are expected to see important changes throughout the French Alps, with the average glacier altitude becoming 300 m (RCP 4.5) and 400 m (RCP 8.5) higher than nowadays (Fig. 2a and S3).

**Nonlinear climate-glacier interactions.** Glacier surface mass changes are commonly modelled by relying on empirical linear relationships between PDDs and snow, firn or ice melt$^{8–10,29}$.

Since the climate and glacier systems are known to be nonlinear$^{13}$, we investigate the benefits of using a model treating, among others, PDDs in a nonlinear way in order to simulate annual glacier-wide MB at a regional scale. We compare model runs using a nonlinear deep learning MB model (the reference approach in our study) against a simplified linear machine learning MB model based on the Lasso$^{30}$, i.e. regularized multi-linear regression. Both MB models were trained with exactly the same data, and all other glacier model parameters were unchanged in order to allow isolating the effects of the nonlinearities in the MB. A comparison between the two MB models shows that a nonlinear response to climate forcings is captured by the deep learning MB model, allowing a better representation of
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Projected glacier and climate evolution per glacierized massif between 2015 and 2100 under RCP 4.5. **a** Projected mean glacier altitude evolution between 2015 and 2100. Massifs without glaciers by 2100 are marked with a cross, **b** Glacier ice volume distribution per massif, with its remaining fraction by 2100 (with respect to 2015), **c** Annual glacier-wide MB per massif, **d** Annual snowfall per massif, **e** Annual cumulative positive degree-days (CPDD) per massif. All values correspond to ensemble means under RCP 4.5. Years in white in c-e indicate the disappearance of all glaciers in a given massif.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Nonlinear and linear response to climatic forcing of MB models.</figcaption></figure>

Nonlinear deep learning response and linear Lasso response to a
Cumulative positive degree days (CPDD) anomalies, b winter snowfall, and
c summer snowfall. All climate anomalies are computed with respect to the
1967-2015 mean values. Envelopes indicate ±σ based on results for all 660
glaciers in the French Alps for the 1967-2015 period. The vertical blue and
red lines indicate the distribution of extreme (top 5%) values for all 21st
century projected climate scenarios, with the mean value in the center and
±1σ indicated by dashed lines.

glacier mass changes, including significantly reduced biases for
extreme values (see Methods). A sensitivity analysis of both MB
models revealed nonlinear relationships between PDDs, snowfall
(in winter and summer) and glacier-wide MB, which the linear
model was only able to approximate (r² = 0.41 for the Lasso vs.
r² = 0.76 for deep learning in cross-validation³¹; Fig. 3).
Regarding air temperature forcings, the linear Lasso MB model
was found to be slightly under-sensitive to extreme positive
cumulative PDD (CPDD) and over-sensitive to extreme negative
CPDDs. Conversely, the linear MB model appears to be over-
sensitive to extreme positive and negative snowfall anomalies.
This behaviour is particularly clear for summer snowfall, for
which the differences are the largest (Fig. 3). In summary, the
linear approximations used by the Lasso manage to correctly fit
the main cluster of average values but perform poorly for extreme
values³¹. This has the strongest impact under RCP 2.6, where
positive MB rates are more frequent (Fig. 4), as the linear model
tends to over-estimate positive MB rates both from air tem-
perature and snowfall (Fig. 3). When using the linear MB model
(Lasso), glaciers are close to reaching an equilibrium with the
climate in the last decades of the century, which is not the case for

the nonlinear MB model (deep learning). Under warmer condi-
tions (RCP 8.5), the differences between the linear and nonlinear
MB model become smaller, as the topographical feedback from
glacier retreat compensates for an important fraction of the losses
induced by the late century warmer climate (Fig. S8 and Fig. 4).
This behaviour is expected for mountain glaciers, as they are
capable of retreating to higher altitudes, thus producing a positive
impact on their glacier-wide MB (Fig. 5). Alternatively, flatter
glaciers (i.e. ice caps) that are found in other glacierized regions
such as the Arctic, where the largest volumes of glacier ice (other
than the ice sheets) are stored³², cannot retreat to higher eleva-
tions. This means that these flatter ice bodies, under a warming
climate, will be subject to higher temperatures than their steeper
counterparts. In order to investigate the implications of these
results for flat glaciers, we performed additional synthetic
experiments in order to reproduce this lack of topographical
feedback (Fig. 4). Our results point out that this lack of topo-
graphical feedback leads to an increased frequency of extreme
negative MB rates and to more pronounced differences between
the nonlinear and linear MB models (Figs. 4e and 5). Our pre-
vious work³¹ has shown that linear MB models can be correctly
calibrated for data around the mean temperature and precipita-
tion values used during training, giving similar results and per-
formance to deep learning. However, as shown in our previous
work and confirmed here, the accuracy of linear models drasti-
cally drops as soon as the input climate data diverges from the
mean cluster of values used for training. A similar behaviour is
observed when comparing temperature-index models to more
complex models (e.g. energy balance), with differences increasing
when the conditions considerably differ from the calibration
period³³. The lower fraction of variance explained by linear
models is present under all climate scenarios. Both models agree
around the average values seen during training (i.e. -0.78 m.w.e.
a⁻¹), but when conditions deviate from this mean training
data centroid, the Lasso can only linearly approximate the
extremes based on the linear trend set on the main cluster of
average values (Fig. 4). Overall, this results in linear MB models
overestimating both extreme positive (Fig. 4a, b) and negative
(Fig. 4e) MB rates.

We further assessed the effect of MB nonlinearities by
comparing our simulated glacier changes with those obtained
from other glacier evolution studies from the literature, which
rely on temperature-index models for MB modelling. Previous
studies on 21st century large-scale glacier evolution projections
have covered the French Alps⁷,⁸. Here, we compare our results
with those from a recent study that focused on the European
Alps¹⁰. In that study, a temperature-index model with a separate
degree-day factor (DDF) for snow and ice is used, resulting in
piecewise linear functions able to partially reproduce nonlinear
MB dynamics. Both the Lasso and the temperature-index
MB model rely on linear relationships between PDDs, solid
precipitation and MB. Therefore, their sensitivities to the
projected 21st century increase in PDDs are linear. Despite the
differences in the two modelling approaches (Table S2), both
regional glacier volume projections present relatively similar
results by the end of the century, with volume differences ranging
between 14% for RCP 2.6 to less than 2% for RCP 4.5 (Fig. S7).
Nonetheless, a close inspection of the annual glacier-wide MB
rates from both models reveals similar patterns to those found
when comparing deep learning and Lasso approaches (Figs. S5
and S6). Despite the existence of slightly different trends during
the first half of the century, both the Lasso and the temperature-
index model react similarly under RCP 4.5 and 8.5 during the
second half of the century, compared to the deep learning model.
The two models with linear MB responses to PDDs and
accumulation simulate more positive MB rates under RCP 2.6,
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Effect of deep learning nonlinearities on glacier mass balance projections. Average cumulative MB projections of French Alpine glaciers with a nonlinear deep learning vs. a linear Lasso model for 29 climate scenarios; **a** with topographical feedback (allowing for glacier retreat) and **e** without topographical feedback (synthetic experiment with constant mean glacier altitude). The projections without glacier geometry adjustment explore the behaviour of glaciers which cannot retreat to higher elevations (i.e. ice cap-like behaviour). **b, c, d** and **f, g, h** annual glacier-wide MB probability distribution functions for all *n* scenarios in each RCP. Vertical axes are different for the two analyses. Photographs taken by Simo Räsänen (Bossons glacier, European Alps, CC BY-SA 3.0) and Doug Hardy (Quelccaya ice cap, Andes, CC BY-SA 4.0).</figcaption></figure>

highlighting their over-sensitivity to negative air temperature anomalies and positive snowfall anomalies (Fig. S6).

Discussion

Studies have warned about the use of temperature-index models for snow and ice projections under climate change for

decades³⁴⁻³⁶. Temperature-index models are known to be over-sensitive to temperature changes, mainly due to important differences in the processes contributing to future warming. This oversensitivity directly results from the fact that temperature-index models rely on linear relationships between PDDs and melt and that these models are calibrated with past MB and climate data. In the past, shortwave radiation represented a more
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Main implications of nonlinear mass balance response to future warming and topographical feedback for mountain glacier and ice cap projections. Glacier topography is a crucial driver of future glacier projections and is expected to play an important role in determining the magnitude that nonlinearities will have on the mass balance signal: ice caps and large flatter glaciers are expected to be more influenced by these nonlinear sensitivities than steep mountain glaciers in a warming climate. Graphics inspired by Hock and Huss⁴⁰.</figcaption></figure>

important fraction in the glacier surface energy budget than the energy fluxes directly related to air temperature (e.g. longwave radiation budget, turbulent fluxes), in comparison with a future warmer climate. Indeed, the projected 21st century warming will lead to increasing incoming longwave radiation and turbulent fluxes, with no marked future trends in the evolution of short-wave radiation³⁷. This will reduce the importance of shortwave radiation for future ablation rates, and it is expected to result in a reduction in values of degree-day factors (DDFs) and therefore a significant change in melt sensitivity to air temperature variations³⁶. This behaviour has already been observed for the European Alps, with a reduction in DDFs for snow during the ablation season of -7% per decade³⁴. Changes in DDFs with respect to air temperature also strongly depend on albedo, with ice presenting a substantially more nonlinear response than snow. Ice melt sensitivity to PDDs strongly decreases with increasing summer temperatures, whereas snow melt sensitivity changes at a smaller rate³⁴. These trends explored with energy balance models from the literature correspond to the behaviour captured by our deep learning MB model, with a clearly less sensitive response of glacier-wide MB to extreme climate forcings, particularly in summer (Fig. 3c). This reduced sensitivity is captured through the response to summer snowfall anomalies, since the sensitivity to positive CPDD anomalies is quite similar for the linear and nonlinear models, as it encompasses both the accumulation and ablation seasons (Fig. 3a).

At this point, it is important to clarify the different ways of treating PDDs in the Lasso and the temperature-index MB models analysed in this study in order to justify analogies. The temperature-index model includes up to three different DDFs, for ice, firn and snow, resulting in three parameters. Alternatively, the Lasso model used here includes 13 DDFs: one for the annual CPDDs and 12 for each month of the hydrological year. Due to the statistical nature of the Lasso model, the response to snowfall anomalies is also highly influenced by variations in PDDs (Fig. 3c). This results in a higher complexity of the Lasso compared to a temperature-index model. Nonetheless, since they are

both linear, their calibrated parameters establishing the sensitivity of melt and glacier-wide MB to temperature variations remain constant over time. This is not the case for the nonlinear deep learning MB model, which captures the nonlinear response of melt and MB to increasing air temperatures, thus reducing the MB sensitivity to extreme positive and negative air temperature and summer snowfall anomalies (Fig. 3). Moreover, these differences between nonlinear and linear models appear to come from an over-sensitivity of linear models to increasing ablation season air temperatures, when ice is exposed in a large fraction of glaciers. Interestingly, this matches the nonlinear, less sensitive response to summer snowfall in the ablation season of our deep learning model (Fig. 3c), which is directly linked to summer air temperatures and has a strong influence on surface albedo. Conversely, during the accumulation season, glaciers are mostly covered by snow, with a much higher albedo and a reduced role of shortwave radiation in the MB that will persist even under climate change. This translates into a more linear response to air temperature changes compared to the ablation season (Fig. 3b). Despite their limitations, temperature-index models, owing to their simplicity and parsimonious data requirements, have been widely used for large-scale glacier projections⁷,⁸. Here, with our newly presented approach, we were able for the first time to quantify the effect that stationary parameters in temperature-index mass balance models have on transient glacier evolution. Our results serve as a strong reminder that the outcomes of existing large-scale glacier simulations should be interpreted with care, and that newly available techniques (such as the nonlinear deep learning approach presented here) and observations (e.g. on various mass balance and radiation components) are opening the door for updated and better constrained projections.

Our results also highlight the important role played by glacier geometry adjustment under changing climatic conditions, which is typical of mountain glaciers³⁸. Our analysis suggests that due to this positive impact on the MB signal, only relevant differences are observed between nonlinear and linear MB models for the lowest emission climate scenarios (Fig. 4a). However, many
glacierized regions in the world present different topographical setups, with flatter glaciers, commonly referred to as ice caps, covering the underlying terrain³⁹. Such ice caps cannot retreat to higher elevations in a warming climate, which inhibits this positive impact on MB⁴⁰ (Fig. 5). In fact, in many cases the surface lowering into warmer air causes this impact on the MB to be negative, further enhancing extreme negative mass balance rates. In order to investigate the effects of MB nonlinearities on flatter glaciers, we conducted a synthetic experiment using the French Alps dataset. We ran glacier evolution projections for both the deep learning and Lasso MB models, but we kept the glacier geometry constant, thus preserving the glacier centroid where the climate data is computed constant through time. With this setup, we reproduced the ice cap-like behaviour with a lack of topographical adjustment to higher elevations. The effect of glaciers shrinking to smaller extents is not captured by these synthetic experiments, but this effect is less important for flat glaciers that are dominated by thinning (Fig. 5). Additionally, glacier surface area was found to be a minor predictor in our MB models³¹. These synthetic experiments suggest that, for equal climatic conditions, flatter glaciers and ice caps will experience substantially more negative MB rates than steeper mountain glaciers. This translates into more frequent extreme negative MB rates, and therefore greater differences due to nonlinearities for the vast majority of future climate scenarios (Fig. 4e).

These different behaviours and resulting biases can potentially induce important consequences in long-term glacier evolution projections. On the one hand, MB nonlinearities for mountain glaciers appear to be only relevant for climate scenarios with a reduction in greenhouse gases emissions (Fig. 4a). The linear Lasso MB model suggests a stabilization of glacier evolution, reaching neutral MB rates by the end of the century. This behaviour is not observed with the nonlinear model, hinting at a positive bias of linear MB models under RCP 2.6. When comparing our deep learning simulations with those from the Lasso, we found average cumulative MB differences of up to 17% by the end of the century (Fig. S5b). For intermediate and pessimistic climate scenarios, no significant differences were found (Fig. S5c-f), except for the largest glaciers (e.g. Mer de Glace, 29 km² in 2015), which did show important differences under RCP 8.5 (up to 75%), due to their longer response time. This suggests that linear MB models are adequate tools for simulating MB of mountain glaciers with important topographical adjustment, with the only exception being the most optimistic climate scenarios and glaciers with long response times. On the other hand, for flatter glaciers large differences between deep learning and Lasso are obtained for almost all climate scenarios (Fig. 4e). Since these flatter glaciers are more likely to go through extreme negative MB rates, nonlinear responses to future warming play a more important role, producing cumulative MB differences of up to 20% by the end of the century (Fig. S5h, j, l). Therefore, linear MB models present more limitations for projections of ice caps, showing a tendency to negative MB biases. As previously mentioned, here these differences are computed at regional level for a wide variety of glaciers. Differences for individual glaciers can be much more pronounced, as large and flat glaciers will have topoclimatic configurations that produce more extreme MB rates than small and steep glaciers with a short response time.

The main uncertainties in future glacier estimates stem from future climate projections and levels of greenhouse gas emissions (differences between RCPs, GCMs, and RCMs), whose relative importance progressively increases throughout the 21st century. With a secondary role, glacier model uncertainty decreases over time, but it represents the greatest source of uncertainty until the middle of the century⁸. Taking into account that for several regions in the world about half of the glacierized volume will be

lost during this first half of the 21st century, glacier models play a major role in the correct assessment of future glacier evolution. The two recent iterations of the Glacier Model Intercomparison Project (GlacierMIP⁷,⁸) have proved a remarkable effort to aggregate, compare and understand global glacier evolution estimates and their associated uncertainties. Despite the existence of a wide variety of different approaches to simulate glacier dynamics, all glacier models in GlacierMIP rely on MB models with linear relationships between PDDs and melt, and precipitation and accumulation. Some of these models use a single DDF, while others have separate DDFs for snow and ice, producing a piecewise function composed of two linear sub-functions that can partially account for nonlinear MB dynamics depending on the snowpack. As we have previously shown, these models present a very similar behaviour to the linear statistical MB model from this study (Fig. 4 vs. S5). Interestingly, our analysis indicates that more complex models using separate DDFs for ice, firn and snow might introduce stronger biases than more simple models using a single DDF. DDFs are known to vary much less with increasing temperatures for intermediate values of albedo (i.e. 0.5) than lower values typical from ice³⁴. Consequently, a simple MB model with a single DDF (e.g. the Open Global Glacier Model - OGGM⁹) is likely to be less affected by an over-sensitivity to future warming than a more complex model with dedicated DDFs for ice, snow, and firn. This creates an interesting dilemma, with more complex temperature-index MB models generally outperforming simpler models for more climatically homogeneous past periods but introducing important biases for future projections under climate change. Our results suggest that, except for the lowest emissions climate scenarios and for large glaciers with long response times, MB models with linear relationships for PDDs and precipitation are suitable for mountain glaciers with a marked topographical feedback. On the one hand, this improves our confidence in long-term MB projections for steep glaciers made by most GlacierMIP models for intermediate and high emissions climate scenarios. On the other hand, ice caps present a different response to future warming, with our results suggesting a negative MB bias by models using linear PDD and accumulation relationships. Ice caps in the Canadian Arctic, the Russian Arctic, Svalbard, and parts of the periphery of Greenland are major reservoirs of ice, as well as some of the biggest expected contributors to sea level rise outside the two polar ice sheets⁷. On top of that, they happen to be among the glacierized regions with the largest projected uncertainties⁸. Together with recent findings by another study⁴¹ highlighting the increased uncertainties in ice thickness distribution estimates of ice caps compared to mountain glaciers, our results raise further awareness on the important uncertainties in glacier projections for ice caps. These conclusions drawn from these synthetic experiments could have large implications given the important sea-level contribution from ice cap-like ice bodies⁸. However, to further investigate these findings, experiments designed more towards ice caps, and including crucial mechanisms such as ice-ocean interactions and thermodynamics, should be used for this purpose.

In this study, we demonstrated the advantages of using deep learning to model glacier MB at regional scales, both in terms of variance and bias. Nonetheless, a better understanding of the underlying processes guiding these nonlinear behaviours at large geographical scales is needed. The machine learning models used in this study are useful to highlight and quantify how nonlinearities in MB affect climate-glacier interactions, but are limited in terms of process understanding. At present, using complex surface energy balance models for large-scale glacier projections is not feasible yet, mainly due to the lack of input data. Therefore, an alternative nonlinear parameterization for the reduction in MB sensitivity under increasing air temperatures would be useful.
This is particularly important for the ablation season and for ice DDFs, which need to accommodate the progressively decreasing role that shortwave radiation will play in the future glacier surface energy budget under warmer conditions. New methods bridging the gap between domain-specific equations and machine learning are starting to arise<sup>42</sup>, which will play a crucial role in further investigating the physical processes driving these nonlinear climate-glacier interactions.

By unravelling nonlinear relationships between climate and glacier MB, we have demonstrated the limitations of linear statistical MB models to represent extreme MB rates in long-term projections. Our analyses suggest that these limitations can also be translated to temperature-index MB models, as they share linear relationships between PDDs and melt, as well as precipitation and accumulation. By performing glacier projections both with mountain glaciers in the French Alps and a synthetic experiment reproducing ice cap-like behaviour, we argue that the limitations identified here for linear models will also have implications for many other glacierized regions in the world. Uncertainties of existing projections of future glacier evolution are particularly large for the second half of the 21st century due to a large uncertainty on future climatic conditions. Our results indicate that these uncertainties might be even larger than we previously thought, as linear MB models are introducing additional biases under the extreme climatic conditions of the late 21st and 22nd centuries. Through synthetic experiments, we showed that the associated uncertainties are likely to be even more pronounced for ice caps, which host the largest reserves of ice outside the two main ice sheets<sup>32</sup>. This implies that current global glacier mass loss projections are too low for the lowest emissions climate scenarios and too high for the highest emissions ones, which has direct consequences for related sea-level rise and water resources projections.

## Methods

**Glacier mass balance modelling.** Glacier-wide MB is simulated annually for individual glaciers using deep learning (i.e. a deep artificial neural network) or the Lasso (regularized multilinear regression)<sup>30</sup>. This modelling approach was described in detail in a previous publication dedicated to the methods, where the ALPine Parameterized Glacier Model (ALPGM<sup>43</sup>) was presented<sup>31</sup>. ALPGM uses a feed-forward fully connected multilayer perceptron, with an architecture (40-20-10-5-1) with Leaky-ReLu<sup>44</sup> activation functions and a single linear function at the output. A He uniform initialization<sup>45</sup> was used for the network parameters. The smallest best performing architecture was used, in order to find a good balance between predictive power, speed, and extrapolation outside the training data. The training was performed with an RMSprop optimizer, batch normalization<sup>46</sup>, and we used both dropout and Gaussian noise in order to regularize it. The Lasso<sup>30</sup>, used for the linear mass balance model, is a linear regression analysis method which shrinks model parameters, thus performing both variable selection and regularization. A dataset of 32 glaciers with direct annual glacier-wide MB observations and remote sensing estimates was used to train the models. For these 32 glaciers, a total of 1048 annual glacier-wide MB values are available, covering the 1967–2015 period with gaps. In order to simulate annual glacier-wide MB values, (a) topographical and (b) climate data for those glaciers and years were compiled for each of the 1048 glacier-year values. (a) Topographical predictors were computed based on the glaciers' annually updated digital elevation model (DEM). These predictors are composed of: the mean glacier altitude, maximum glacier altitude, slope of the lowermost 20% altitudinal range of the glacier, glacier surface area, latitude, longitude and aspect. (b) Climate predictors are based on climatic anomalies computed at the glaciers' mean altitude with respect to the 1967–2015 reference period mean values. Models were trained using the SAFRAN reanalysis dataset<sup>47</sup>, including observations of mountain regions in France for the 1958–2015 period. This reanalysis is specifically designed to represent meteorological conditions over complex mountain terrain, being divided by mountain massif, aspect and elevation bands of 300 m. Winter climate data are computed between October 1 and March 31, and summer data between April 1 and September 30. Climate predictors consist of: the annual CPDD, winter snowfall, summer snowfall, monthly mean temperature and monthly snowfall. This creates a total of 34 input predictors for each year (7 topographical, 3 seasonal climate, and 24 monthly climate predictors).

In order to avoid overfitting, MB models were thoroughly cross-validated using all data for the 1967–2015 period in order to ensure a correct out-of-sample performance. Three different types of cross validation were performed: a Leave-One-Glacier-Out (LOGO), a Leave-One-Year-Out (LOYO) and a Leave-Some-Years-and-

Glaciers-Out (LSYGO). Each one of these cross-validations served to evaluate the model performance for the spatial, temporal and both dimensions, respectively. When working with spatiotemporal data, it is imperative to respect spatial and temporal data structures during cross-validation in order to correctly assess an accurate model performance<sup>48</sup>. With this cross-validation we determined a deep learning MB model spatiotemporal (LSYGO) RMSE of 0.59 m.w.e. $a^{-1}$ and an $r^2$ of 0.69, explaining 69% of the total MB variance. Alternatively, the Lasso MB model displayed an RMSE of 0.85 m.w.e. $a^{-1}$ and an $r^2$ of 0.35<sup>31</sup>. Simulations for projections in this study were made by generating an ensemble of 60 cross-validated models based on LSYGO. Each one of these models was created by training a deep learning model with the full dataset except all data from a random glacier and year, and evaluating the performance on these hidden values. This ensures that the model is capable of reproducing MB rates for unseen glaciers and years. Simulations were then performed by averaging the outputs of each one of the 60 ensemble members. This approach is known as a cross-validation ensemble<sup>49</sup>. Future projections of glacier-wide MB evolution were performed using climate projections from ADAMONT<sup>25</sup>. This dataset applies a statistical adjustment specific to French mountain regions based on the SAFRAN dataset, to EURO-CORDEX<sup>26</sup> GCM-RCM-RCP members, covering a total of 29 different future climate scenarios for the 2005–2100 period. This adjustment represents a major improvement over most climate data used to force regional and global glacier models. The high spatial resolution enables a detailed representation of mountain weather patterns, which are often undermined by a coarser resolution climate datasets.

**Glacier geometry evolution.** A well-established parametrization based on empirical functions<sup>50</sup> was used in order to redistribute the annually simulated glacier-wide mass changes over each glacier. This parametrization reproduces in an empirical manner the changes in glacier geometry due to the combined effects of ice dynamics and MB. As for the MB modelling approach, a detailed explanation on this method can be found in a previous dedicated paper on the methods<sup>31</sup>. In our model, we specifically computed this parameterized function for each individual glacier larger than 0.5 km<sup>2</sup>, representing 80% of the total glacierized area in 2015, using two DEMs covering the whole French Alps: a photogrammetric one in 1979 and a SPOT-5 one in 2011. We previously demonstrated that this period is long enough to represent the secular trend of glacier dynamics in the region<sup>31</sup>. Both DEMs were resampled and aligned at a common spatial resolution of 25 m. For each glacier, an individual parameterized function was computed representing the differences in glacier surface elevation with respect to the glacier's altitude within the 1979–2011 period. This method has the advantage of including glacier-specific dynamics in the model, encompassing a wide range of different glacier behaviours. Using this function, the glacier-specific ice thickness and the DEM are updated every year, adjusting the 3D geometry of each glacier. This enables the recalculation of every topographical predictor used for the MB model, thus updating the mean glacier altitude at which climate data for each glacier are retrieved. This annual geometry adjustment accounts for the effects of glacier retreat on the climate signal received by glaciers. Glaciers smaller than 0.5 km<sup>2</sup> often display a high climate imbalance, with their equilibrium line being higher than the glacier's maximum altitude. Such glaciers are often remnants of the Little Ice Age, and mainly lose mass via non-dynamic downwasting<sup>51</sup>. For such cases, we assumed that ice dynamics no longer play an important role, and the mass changes were applied equally throughout the glacier.

The performance of this parametrization was validated in a previous study, indicating a correct agreement with observations<sup>31</sup>. The dataset of initial glacier ice thickness, available for the year 2003, determines the starting point of our simulations. We performed a validation simulation for the 2003–2015 period by running our model through this period and comparing the simulated glacier surface area of each of the 32 glaciers with MB to observations from the 2015 glacier inventory<sup>16,52</sup>. Then, we ran multiple simulations for this same period by altering the initial ice thickness by ±30% and the glacier geometry update parametrizations by ±10%, according to the estimated uncertainties of each of the two methods<sup>31</sup>. These results revealed that the main uncertainties on glacier simulations arise from the initial ice thickness used to initialize the model. This is well in agreement with the known uncertainties of glacier evolution models, with glacier ice thickness being the second largest uncertainty after the future GCM-RCM-RCP climate members used to force the model<sup>29</sup>. Glacier ice thickness observations are available for four different glaciers in the regions, which were compared to the estimates used in this model. Ice thickness accuracy varied significantly, with an overall correct representation of the ice distribution but with local biases reaching up to 100%. The ice thickness data for two of the largest glaciers in the French Alps were modified in order to improve data quality. Ice thickness data for Argentière glacier (12.27 km<sup>2</sup> in 2015) was taken from a combination of field observations (seismic, ground-penetrating radar or hot-water drilling<sup>53</sup>) and simulations<sup>32</sup>. The estimated ice thickness for Mer de Glace (28.87 km<sup>2</sup> in 2015) was increased by 25% in order to correct the bias with respect to field observations<sup>31</sup>. Since these two glaciers are expected to be some of the few large glaciers that will survive the 21st century climate, an accurate representation of their initial ice thickness has an important effect on the estimates of remaining ice.

**Model comparison and extraction of nonlinearities.** The nonlinearities present in the simulated annual glacier-wide MB values were assessed by running two different glacier simulations with two different MB models. The advantage of this
method is that by only changing the MB model, we can keep the rest of the model
components (glacier dynamics and climate forcing) and parameters the same in
order to have a controlled environment for our experiment. Therefore, we were
capable of isolating the different behaviours of the nonlinear deep learning model
and a linear machine learning model based on the Lasso³⁰. Both machine learning
MB models were trained with exactly the same data coming from the 1048 annual
glacier-wide MB values, and both were cross-validated using LSYGO. In order to
investigate the effects of MB nonlinearities on ice caps, we performed the same type
of comparison between simulations, but the glacier geometry update module
described in the "Glacier geometry evolution" section was deactivated. This syn-
thetic setup allowed us to reproduce the climatic conditions to be undergone by
most ice caps, with their mean surface altitude hardly evolving through time. This
removes the topographical feedback typical from mountain glaciers, and repro-
duces the more extreme climate conditions that ice caps are likely to endure
through the 21st century⁴⁰. This synthetic experiment is an approximation of what
might occur in other glacierized regions with ice caps. In many aspects, it might be
too optimistic, as many ice caps will have a negative impact on MB through
thinning, bringing their mean surface elevation to lower altitudes, thus further
warming their perceived climate. This means that these differences linked to MB
nonlinearities observed in this experiment could be even greater for such ice caps.
However, the impact of different climate configurations, such as a more continental
and drier climate or a more oceanic and humid climate, would certainly have an
impact on the results, albeit a much less important one than the lack of topo-
graphical feedback explored here. Our synthetic experiment does not account for
glacier surface area shrinking either, which might have an impact on the glacier-
wide MB signal. Nevertheless, we previously demonstrated that glacier surface area
is not an important predictor of MB changes in our models²⁹, and ice caps evolve
mostly through thinning and not shrinking (Fig. 5).

Additionally, the specific responses of the deep learning and Lasso MB models to air temperature and snowfall were extracted by performing a model sensitivity analysis. Since the neural network used here virtually behaves like a black box, an alternative way is needed to understand the model's behaviour. In order to do so, we applied a deterministic sampling process as a sensitivity analysis to both the deep learning and the Lasso MB models. For that, a dataset of input predictors covering all the glaciers in the French Alps for the 1967–2015 period was generated from a past MB reconstruction study¹⁵. Multiple copies of this dataset were created, and for each individual copy a single predictor (i.e. CPDD, winter snowfall or summer snowfall) was modified for all glaciers and years. This allows us to assess the MB models responses at a regional scale to changes in individual predictors (Fig. 3). Regarding air temperature, a specific CPDD anomaly ranging from –1500 PDD to +1500 PDD in steps of 100 PDD was prescribed to all glaciers for each dataset copy. Since both MB models also include monthly temperature data as predictors, this CPDD anomaly was distributed evenly between the ablation season (April 1–September 30), following the expected increase in mostly summer temperatures instead of winter temperatures in the future (Fig. 1). Tests were performed distributing the CPDD anomalies equally among all months of the year with very similar results. The same was done with winter snowfall anomalies, ranging between –1500 mm and +1500 mm in steps of 100 mm, and summer snowfall anomalies, ranging between –1000 mm and +1000 mm in steps of 100 mm. The anomaly in snowfall was evenly distributed for every month in the accumulation (October 1–April 31) and ablation seasons, respectively. This experiment enabled the exploration of the response to specific climate forcings of a wide range of glaciers of different topographical characteristics in a wide range of different climatic setups, determined by all meteorological conditions from the years 1967–2015 (Fig. 3).

Alternatively, the comparisons against an independent large-scale glacier evolution model were less straightforward to achieve. GloGEMflow¹⁰ is a state-of-the-art global glacier evolution model used in a wide range of studies, including the second phase of GlacierMIP⁷,⁸. Several differences are present between ALPGM, the model used in this study, and GloGEMflow (Table S2), which hinder a direct meaningful comparison between both. In order to overcome these differences, some adaptations were performed to the GloGEMflow output, accompanied with some hypotheses to ensure a realistic comparison. The first main difference is related to the climate data used to force the models. GloGEMflow relies on EURO-CORDEX ensembles²⁶, whereas ALPGM uses ADAMONT²⁵, an adjusted version of EURO-CORDEX specifically designed for mountain regions. ADAMONT provides climate data at 300 m altitudinal bands and different slope aspects, thus having a significantly higher spatial resolution than the 0.11° from EURO-CORDEX. This implies that specific climatic differences between massifs can be better captured by ALPGM than GloGEMflow. Nonetheless, since the main GCM-RCM climate signal is the same, the main large-scale long-term trends are quite similar. We reduced these differences by running simulations with GloGEMflow using exactly the same 29 climate members used by ALPGM in this study (Table S1). The initial glacier ice thickness data for the year 2003 also differs slightly between both models. The original ice thickness estimates of the methods used by both models are different¹⁰,³², and for ALPGM we performed some additional modifications to the two largest glaciers in the French Alps (see ‘Glacier geometry evolution’ for details). Despite these differences, the average altitude difference of the glaciers between both models is never higher than 50 m (Fig. S10). Since in ALPGM the climate forcing of glaciers is extracted at the mean glacier altitude, we do not expect these altitude differences to drive important MB

differences between models. Another source of discrepancy between both models
comes from the different MB data used to calibrate or train the MB models.
GloGEMflow has been previously applied in a study over the whole European Alps,
and its temperature-index model was mainly calibrated with MB data from the
Swiss Alps. Swiss glaciers have displayed less negative MB rates than French
glaciers during the last decades, thus likely introducing a bias in simulations
specific to the French Alps. In order to improve the comparability between both
models, a MB bias correction was applied to GloGEMflow's simulated MB, based
on the average annual MB difference between both models for the 2003–2015
period (–0.4 m.w.e. a⁻¹) over the French Alps. Finally, there are differences as well
in the glacier dynamics of both models, with ALPGM using a glacier-specific
parameterized approach and GloGEMflow explicitly reproducing the ice flow
dynamics. Nonetheless, these differences have been shown to be rather small,
having a lower impact on results than climate forcings or the initial glacier ice
thickness¹⁰.

Data availability

The model output data generated in this study have been deposited in netCDF and CSV
format in a Zenodo repository under accession code Creative Commons Attribution 4.0
International. https://zenodo.org/record/5549758.

Code availability

The source code of the glacier model can be freely accessed in the following repository:
https://github.com/JordiBolivar/ALPGM.

Received: 7 April 2021; Accepted: 3 January 2022;
Published online: 20 January 2022

**References**

1. IPCC. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (2018).

2. Huss, M. et al. Toward mountains without permanent snow and ice: mountains without permanent snow and ice. *Earths Future* **5**, 418–435 (2017).

3. Zemp, M. et al. Global glacier mass changes and their contributions to sea-level rise from 1961 to 2016. *Nature* **568**, 382–386 (2019).

4. Immerzeel, W. W. et al. Importance and vulnerability of the world's water towers. *Nature* **577**, 364–369 (2020).

5. Farinotti, D., Round, V., Huss, M., Compagno, L. & Zekollari, H. Large hydropower and water-storage potential in future glacier-free basins. *Nature* **575**, 341–344 (2019).

6. Caavy-Fraunié, S. & Dangles, O. A global synthesis of biodiversity responses to glacier retreat. *Nat. Ecol. Evol.* **3**, 1675–1685 (2019).

7. Hock, R. et al. GlacierMIP – A model intercomparison of global-scale glacier mass-balance models and projections. *J. Glaciol.* **65**, 453–467 (2019).

8. Marzeion, B. et al. Partitioning the uncertainty of ensemble projections of global glacier mass change. *Earths Future* https://doi.org/10.1029/2019EF001470 (2020).

9. Maussion, F. et al. The Open Global Glacier Model (OGGM) v1.1. *Geosci. Model Dev.* **12**, 909–931 (2019).

10. Zekollari, H., Huss, M. & Farinotti, D. Modelling the future evolution of glaciers in the European Alps under the EURO-CORDEX RCM ensemble. *The Cryosphere* **13**, 1125–1146 (2019).

11. Hock, R. Temperature index melt modelling in mountain areas. *J. Hydrol.* **282**, 104–115 (2003).

12. Réveillet, M. et al. Relative performance of empirical and physical models in assessing the seasonal and annual glacier surface mass balance of Saint-Sorlin Glacier (French Alps). *The Cryosphere* **12**, 1367–1386 (2018).

13. Steiner, D., Walter, A. & Zumbühl, H. J. The application of a non-linear backpropagation neural network to study the mass balance of Grosse Aletschgletscher, Switzerland. *J. Glaciol.* **51**, 313–323 (2005).

14. Clarke, G. K. C., Berthier, E., Schoof, C. G. & Jarosch, A. H. Neural networks applied to estimating subglacial topography and glacier volume. *J. Clim.* **22**, 2146–2160 (2009).

15. Boljar, J., Rabatel, A., Gouttevin, I. & Galiez, C. A deep learning reconstruction of mass balance series for all glaciers in the French Alps: 1967–2015. *Earth Syst. Sci. Data* **12**, 1973–1983 (2020).

16. Paul, F. et al. Glacier shrinkage in the Alps continues unabated as revealed by a new glacier inventory from Sentinel-2. *Earth Syst. Sci. Data* **12**, 1805–1821 (2020).

17. Hugonnet, R. et al. A globally complete, spatially and temporally resolved estimate of glacier mass change: 2000 to 2019. https://meetingorganizer.copernicus.org/EGU2020/EGU2020-20908.html (2020) https://doi.org/10.5194/egusphere-egu2020-20908.
18. Salim, E., Ravanel, L., Deline, P. & Gauchon, C. A review of melting ice adaptation strategies in the glacier tourism context. *Scand. J. Hosp. Tour.* **21**, 229–246 (2021).

19. Spandre, P. et al. Winter tourism under climate change in the Pyrenees and the French Alps: relevance of snowmaking as a technical adaptation. *The Cryosphere* **13**, 1325–1347 (2019).

20. Magnin, F., Haeberli, W., Linsbauer, A., Deline, P. & Ravanel, L. Estimating glacier-bed overdeepenings as possible sites of future lakes in the de-glaciating Mont Blanc massif (Western European Alps). *Geomorphology* **350**, 106913 (2020).

21. Carlson, B. Z. et al. Monitoring the Seasonal hydrology of alpine wetlands in response to snow cover dynamics and summer climate: a novel approach with sentinel-2. *Remote Sens.* **12**, 1959 (2020).

22. Robinson, C. T., Thompson, C. & Freestone, M. Ecosystem development of streams lengthened by rapid glacial recession. *Fundam. Appl. Limnol. Arch. Für Hydrobiol.* **185**, 235–246 (2014).

23. Vincent, C. et al. Common climatic signal from glaciers in the European Alps over the last 50 years: Common Climatic Signal in the Alps. *Geophys. Res. Lett.* **44**, 1376–1383 (2017).

24. GLAMOS. Swiss Glacier Mass Balance (release 2019). (2019) https://doi.org/10.18750/MASSBALANCE.2019.R2019.

25. Verfaillie, D., Déqué, M., Morin, S. & Lafaysse, M. The method ADAMONT v1.0 for statistical adjustment of climate projections applicable to energy balance land surface models. *Geosci. Model Dev.* **10**, 4257–4283 (2017).

26. Smiatek, G., Kunstmann, H. & Senatore, A. EURO-CORDEX regional climate model analysis for the Greater Alpine Region: performance and expected future change: climate change in the gar area. *J. Geophys. Res. Atmospheres* **121**, 7710–7728 (2016).

27. Roe, G. H. Orographic precipitation. *Annu. Rev. Earth Planet. Sci.* **33**, 645–671 (2005).

28. Six, D. & Vincent, C. Sensitivity of mass balance and equilibrium-line altitude to climate change in the French Alps. *J. Glaciol.* **60**, 867–878 (2014).

29. Huss, M. & Hock, R. A new model for global glacier change and sea-level rise. *Front. Earth Sci.* **3** (2015).

30. Tibshirani, R. Regression Shrinkage and Selection via the Lasso. *J. R. Stat. Soc. Ser. B Methodol.* **58**, 267–288 (1996).

31. Bolibar, J. et al. Deep learning applied to glacier evolution modelling. *The Cryosphere* **14**, 565–584 (2020).

32. Farinotti, D. et al. A consensus estimate for the ice thickness distribution of all glaciers on Earth. *Nat. Geosci.* **12**, 168–173 (2019).

33. Gabbi, J., Carenzo, M., Pellicciotti, F., Bauder, A. & Funk, M. A comparison of empirical and physically based glacier surface melt models for long-term simulations of glacier response. *J. Glaciol.* **60**, 1140–1154 (2014).

34. Braithwaite, R. J. Positive degree-day factors for ablation on the Greenland ice sheet studied by energy-balance modelling. *J. Glaciol.* **41**, 153–160 (1995).

35. Pellicciotti, F. et al. An enhanced temperature-index glacier melt model including the shortwave radiation balance: development and testing for Haut Glacier d'Arolla, Switzerland. *J. Glaciol.* **51**, 573–587 (2005).

36. Huss, M., Funk, M. & Ohmura, A. Strong Alpine glacier melt in the 1940s due to enhanced solar radiation. *Geophys. Res. Lett.* **36**, L23501 (2009).

37. Bartók, B. et al. Projected changes in surface solar radiation in CMIP5 global climate models and in EURO-CORDEX regional climate models for Europe. *Clim. Dyn.* **49**, 2665–2683 (2017).

38. Zekollari, H., Huss, M. & Farinotti, D. On the Imbalance and Response Time of Glaciers in the European Alps. *Geophys. Res. Lett.* **47** (2020).

39. Consortium, R. G. I. Randolph Glacier Inventory 6.0 (2017) https://doi.org/10.7265/N5-RGI-60.

40. Hock, R. & Huss, M. Glaciers and climate change. In *Climate Change* 157–176 (Elsevier, 2021). https://doi.org/10.1016/B978-0-12-821575-3.00009-8.

41. Millan, R., Mouginot, J., Rabatel, A., & Morlighem, M. Ice velocity and thickness of the world's glaciers. *Nature Geosciences*, https://doi.org/10.1038/s41561-021-00885-z (2022).

42. Rackauckas, C. et al. Universal Differential Equations for Scientific Machine Learning. *ArXiv200104385 Cs Math Q-Bio Stat* (2020).

43. Bolibar, J. ALPGM (ALpine Parameterized Glacier Model) v1.1. (Zenodo, 2020). https://doi.org/10.5281/zenodo.3609136.

44. Xu, B., Wang, N., Chen, T. & Li, M. Empirical Evaluation of Rectified Activations in Convolutional Network. *CoRR abs/1505.00853* (2015).

45. He, K., Zhang, X., Ren, S. & Sun, J. Delving Deep into Rectifiers: Surpassing Human-Level Performance on ImageNet Classification. *2015 IEEE Int. Conf. Comput. Vis. ICCV* (2015) https://doi.org/10.1109/iccv.2015.123.

46. Ioffe, S. & Szegedy, C. *Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift* (2015).

47. Durand, Y. et al. Reanalysis of 47 Years of Climate in the French Alps (1958–2005): Climatology and Trends for Snow Cover. *J. Appl. Meteorol. Climatol.* **48**, 2487–2512 (2009).

48. Roberts, D. R. et al. Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure. *Ecography* **40**, 913–929 (2017).

49. Hastie, T., Tibshirani, R. & Friedman, J. The Elements of Statistical Learning. (Springer, New York, 2009).

50. Huss, M., Jouvet, G., Farinotti, D. & Bauder, A. Future high-mountain hydrology: a new parameterization of glacier retreat. *Hydrol. Earth Syst. Sci.* **14**, 815–829 (2010).

51. Paul, F., Kääb, A., Maisch, M., Kellenberger, T. & Haeberli, W. Rapid disintegration of Alpine glaciers observed with satellite data: disintegration of alpine glaciers. *Geophys. Res. Lett.* **31**, n/a–n/a (2004).

52. Gardent, M., Rabatel, A., Dedieu, J.-P. & Deline, P. Multitemporal glacier inventory of the French Alps from the late 1960s to the late 2000s. *Glob. Planet. Change* **120**, 24–37 (2014).

53. Rabatel, A., Sanchez, O., Vincent, C. & Six, D. Estimation of glacier thickness from surface mass balance and ice flow velocities: a case study on Argentière Glacier, France. *Front. Earth Sci.* **6** (2018).

## Acknowledgements

We acknowledge the more than 50 years of glaciological monitoring performed by the GLACIOCLIM French National Observatory (https://glacioclim.osug.fr), which provided essential observations for our modelling study. This work was funded by the Labex OSUG@2020 (Investissements d'avenir, ANR10 LABX56) and the Auvergne-Rhône-Alpes region through the BERGER project. J.B. was supported by a NWO VIDI grant 016.Vidi.171.063. H.Z. acknowledges the funding received from a EU Horizon 2020 Marie Skłodowska-Curie Individual Fellowship (grant no. 799904) and from the Fonds de la Recherche Scientifique - FNRS (postdoctoral grant - chargé de recherches).

## Author contributions

J.B. developed the main glacier model, performed the simulations, analysed the results, and wrote the paper. A.R. provided glacier mass balance data and performed the glaciological analyses. I.G. contributed to the climate analyses. H.Z. performed simulations with another glacier model, provided results for comparison, and contributed to the glaciological analyses. C.G. contributed to the extraction of nonlinear mass balance responses and to the statistical analysis. All authors provided inputs to the paper and helped to write it.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-28033-0.

**Correspondence** and requests for materials should be addressed to Jordi Bolibar.

**Peer review information** *Nature Communications* thanks Mohd Anul Haq, Lauren Vargo, and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022