1950–2100 climate trends in avalanche activity in Haute-Maurienne, French Alps

Avalanche activity in alpine regions is sensitive to climate change. However, without consistent historical data, it is challenging to estimate past trends in avalanche activity and assess future avalanche scenarios from climate projections. To tackle this challenge, we use avalanche observations and simulated snowpack conditions to train a machine learning gradient-boosting regression model, which predicts the number of avalanches per day. We focus on a small alpine domain with high-quality data: the Haute-Maurienne valley in the French Alps, where avalanche paths span elevations from approximately 1800 to 2700 m a.s.l. First, we demonstrate that accounting for the uncertainties in avalanche occurrence dates and using only the most recent period (2006–2023) with homogeneous observations during the training step is essential for achieving consistent results. We then use this machine learning model to reconstruct the past avalanche activity (1958–2023) from reanalysed meteorological and snow data, and to project future avalanche activity (1950–2100) from a downscaled ensemble of snow-climate simulations. We evaluate climatic trends in avalanche activity using three indicators: the number of avalanches per winter season, the number of avalanches per month, and the annual maximum number of avalanches in one week, which quantifies the largest avalanche cycles. Based on reanalysed snow-climate simulations, the model estimates that avalanche activity decreased in the past: the mean number of avalanches per year declined by 4.8 % per decade between 1958 and 2023, with a stronger decrease in spring avalanche activity, and the 30-year return level associated with large avalanche cycles decreased at a slower rate of 1.4 % per decade. In the future, avalanche activity is also expected to decrease. For the emission scenarios RCP4.5 and RCP8.5, the annual number of avalanches is expected to decrease by 4.6 % per decade and 8.8 % per decade, respectively, mainly due to a reduction in spring avalanche activity. Large avalanche cycles, quantified by the 30-year return level, are also expected to decrease in intensity but at slower rates: 2.2 % per decade for RCP4.5 and 5.5 % per decade for RCP8.5. This study quantifies the impact of climate change on avalanche activity in an exemplary alpine valley. It demonstrates that combining statistical learning with climate simulations can help produce reference scenarios for mitigation strategies in high mountain environments.

1 Introduction

Snow avalanches are rapid mass movements of snow down a steep slope (Schweizer et al., 2003). They represent a significant natural hazard in snow-covered mountainous regions, threatening human lives, settlements, and infrastructure. Studying avalanche risk is essential in populated areas of mountainous regions, such as the European Alps, due to the economic and human stakes involved (Fuchs et al., 2007; Cappabianca et al., 2008; Eckert et al., 2012; Techel et al., 2016). In these areas, mitigation strategies are mandatory and mainly rely on heavy protective structures or restrictive land-use policies, which are designed to remain effective over several decades (Höller, 2007). Therefore, implementing effective risk mitigation strategies requires anticipating potential non-stationarity in risk components, among which is the possible impact of climate change on avalanches (Fuchs et al., 2013; Eckert and Giacona, 2023).

The frequency, magnitude, and spatial distribution of avalanches are closely related to snowpack conditions (Schweizer et al., 2003), i.e. the quantity and quality of the snow on the ground. Their evolution with climate change is the primary driver of climatic trends in avalanche activity and is thus briefly described here for the European Alps (e.g. Dumont et al., 2025). Durand et al. (2009) reported a marked decline in snow depth in the French Alps between 1958 and 2005, with the most significant reductions observed at the lowest elevations and at the end of the winter season. Matiu et al. (2021) quantified this shortening of the snow season by roughly one month since 1950, below 2000 m elevation in the European Alps. The reduction in snow on the ground is primarily attributed to rising air temperatures, which cause precipitation to change from a solid to a liquid state (Serquet et al., 2011). Meanwhile, mean precipitation amounts are almost constant (Bozzoli et al., 2024; Masson and Frei, 2015). Therefore, the impact of climate change is strongly dependent on elevation: at relatively high altitude sites such as Weissfluhjoch (located at 2536 m in the Swiss Alps), mean winter snow depth has remained nearly constant between 1960 and 2023 because the warming in winter has not been sufficient to shift precipitation from snow to rain (Dumont et al., 2025). Future trends were computed from climate simulations downscaled on the complex topography of alpine regions (e.g., Marty et al., 2017; Verfaillie et al., 2018; Kotlarski et al., 2023). Marty et al. (2017) demonstrated that the resulting snow cover changes may be roughly equivalent to an elevation shift of 500 to 1000 m by the end of the century depending on the emission scenarios. Verfaillie et al. (2018) showed that the impact of emission scenarios becomes significant for the second half of the 21st century. While the mean winter snow depth is expected to decrease, climate change may also impact other snowpack properties. The snowpack will be wetter, and rain-on-snow events are expected to increase (Beniston and Stoffel, 2016). Extreme snowfall events will also be impacted: by the end of the century, they are expected to decrease in the Pyrenees (Bonsoms et al., 2025), as in the French Alps under 2700 m, but the 100-year return level of daily snowfall above 2700–3000 m is expected to stay constant or to increase slightly (Le Roux et al., 2023). The decrease in mean winter snow depth might also be accompanied by shifts in the physical properties of the snowpack, as the snow metamorphism is dependent on the temperature gradient (Schneebeli and Sokratov, 2004). Because snow microstructure has a primary influence on snowpack stability (Schweizer et al., 2003), the combined effects of changing snow depth, snow layer properties, and resulting snow microstructure make the overall impact on snowpack stability difficult to assess.

One direct way to estimate past climate trends in avalanche activity is to assess it in a consistent, representative spatial domain using a fixed sampling strategy over several decades. Whereas satellite monitoring is a promising tool (Karas et al., 2022) to acquire data in a systematic way, this new technology lacks temporal coverage to capture anthropogenic climate trends above the decadal internal climate variability (Hafner et al., 2021). Given that systematic measurements are rare, multi-source historical archives (e.g., administrative documents, newspapers) (Giacona et al., 2017, 2022), or proxy records, such as tree rings (Corona et al., 2012), seismic signals (Heck et al., 2018) or geomorphic markers (Johnson and Smith, 2010) are often used (Eckert et al., 2024). For instance, Giacona et al. (2021) used historical records and statistical techniques to examine the impact of long-term past climate change on avalanche activity in the low-elevation Vosges Mountains, France. The study mainly revealed a transition from the late Little Ice Age to the early twentieth century, with a drastic reduction in the annual number of avalanches. Peitzsch et al. (2021) used tree rings to quantify a decline of about 14 % in avalanche activity in the Rocky Mountains between 1950 and 2017. In a higher mountain environment, Ballesteros-Cánovas et al. (2018) demonstrated an increase in avalanche activity in the western Indian Himalayas, using dendrochronology in avalanche paths located between 2600 and 4200 m elevation, but the study was limited to a few slopes. In the French Alps, the Enquête Permanente sur les Avalanches (EPA) reports avalanche activity observed in numerous and specified avalanche paths since the early twentieth century (Bourova et al., 2016). It is thus likely to vary over time with human factors and formal task instructions, creating biases in trends and thereby precluding inference of climatic control. At the scale of entire mountain ranges, biases associated with individual observers may be negligible because large-scale trends may still emerge despite local observational artefacts (Eckert et al., 2010d). However, such biases can become critical when the observation record is based on only a few observers. The temporal resolution of such observations is high – a few days – compared to indirect observations, but still too coarse to capture the correlation to snow conditions evolving within hours. Based on this EPA dataset for the period 1950–2009, Eckert et al. (2013) reported a peak in avalanche activity in the 1980s in the French Alps, followed by a 19 % decrease in the mean number of avalanches per winter through 2009. Due to unknown sampling variations over time, it is, however, almost impossible to evaluate to which extent this overall trend is altered over smaller domains.

Capturing future trends from direct or indirect avalanche observations is by definition impossible. Projecting avalanche activity necessarily requires linking climate simulations forced by emission scenarios to snow conditions and finally to avalanche activity. The main idea is to build a relationship between snow conditions and avalanche activity using either physical knowledge of the involved processes or machine learning techniques on a reference period, where the targeted avalanche activity and predictive snow conditions are known with sufficient accuracy. This impact model can then be applied to snow simulations forced by downscaled climate simulations. This strategy also works on the past, with a forcing based on reanalysed meteorological data. We employ this strategy in the paper and briefly review here the studies which used similar methods to capture climate-driven trends in avalanche activity. Reuter et al. (2025) developed a physically-based model to relate simulated snow conditions to avalanche problems. Applied to the S2M reanalysis (Vernay et al., 2022b) on the French Alps over 1958–2020, they showed a decreasing frequency of persistent weak layer problems, an increasing frequency of new snow problems and that the onset of wet-snow activity now occurs about three weeks earlier than in 1958. This study provides clues about the evolution of the main drivers of avalanche activity but lacks quantification of the severity of the detected avalanche problems. Castebrunet et al. (2014) pioneered the projection of future avalanche activity. Instead of modeling individual avalanche occurrences, they used a Composite Index of avalanche activity derived from the EPA database and from the MEPRA hazard index (Durand et al., 1999). This synthetic indicator, which cannot be associated to a direct observable avalanche activity index, was computed at seasonal and annual time scales over large regions of the French Alps. They developed statistical regression models relating this avalanche activity index to snowpack and meteorological conditions simulated by the SAFRAN – Crocus – MEPRA chain (Durand et al., 1999). Based on the global trajectories of the CMIP3 framework (Nakicenović and Swart, 2000), they showed a general decrease in avalanche activity for the whole domain and observed inconsistent evolutions at high altitude. Ortner et al. (2025) investigated future changes in avalanche risks and the associated monetary loss in the Central Swiss Alps. They built a hazard map based on the avalanche flow RAMMS forced by the projected 3 d snow depth accumulation and temperature, which is then combined with the vulnerability of the exposed infrastructures. The results show a marked decline in the average annual impact of avalanches by the end of the century. However, it should be noted that this integrated approach, from climate simulations to monetary costs, remains relatively simple at the avalanche formation step. Indeed, the study assumes that any three-day snowfall can, without distinction, trigger an avalanche and that factors known to contribute to avalanche formation, such as wind, pre-existing snowpack, and snow wetting, are omitted. Mayer et al. (2024) used a snow cover model driven by downscaled climate simulations to compute future alterations in dry- and wet-snow avalanche occurrences throughout the 21st century across seven sites in the Swiss Alps. Their machine learning model classifies each day as either a non-avalanche day, a dry-snow avalanche day or a wet-snow avalanche day. The results show a decrease in the overall annual number of avalanche days until 2100, by 20 %–40 % relative to the RCP8.5 scenario, and approximately 10 % relative to the RCP4.5 scenario. To date, Mayer et al. (2024) have provided one of the most comprehensive studies of projected changes in avalanche activity under a warming climate. However, the binary classification of avalanche and non-avalanche days does not provide information about the intensity of the avalanche cycle, which is of great importance for mitigation strategies. This projected decrease in the number of avalanche days can be partly attributed to the shortening of the winter season (Kruyt et al., 2023). However, additional processes may also contribute to changes in snowpack stability, including increased air temperature and associated complex modifications of snow stratigraphy, as well as other mechanical effects influencing snow evolution.

In broad terms, current studies suggest that warming temperatures are generally associated with a decrease in average avalanche activity. This decreasing trend is exacerbated where the snowpack vanishes, i.e., at lower elevations or at the end of the snow season (Eckert et al., 2024). However, at higher elevations and for high-impact avalanche cycles, the picture is less clear and lacks support by quantitative findings. Long and consistent observational series are rare, and indirect sources such as historical archives or tree rings provide trends at annual/seasonal resolution only. Human-based inventories introduce observer biases and uncertainties that may blur true climatic signals. Model-based studies, whether physical, statistical, or machine learning, generally simplify key processes and often provide broad patterns without quantifying the actual intensity of avalanche cycles. There are thus many gaps in the current understanding of the impact of climate change on avalanche activity.

In this study, we focus on a populated valley in the French Alps, the Haute-Maurienne, and quantify changes in avalanche activity between 1950 and 2100. This valley is a model territory at the European scale: frequent avalanches affect human activities, it is a relatively high elevation site, and avalanches have been systematically reported in recent decades. We characterize the changes in avalanche activity not only by the average number of avalanches per year but also by the seasonality of avalanche occurrence and the return period of intense avalanche cycles. To capture the link between snow conditions and avalanches, we train a machine learning model on a human-based avalanche inventory and reanalysed weather and snow simulations. The training step accounts for potential biases and uncertainties in the inventory, and the set of predictors is comprehensive enough to capture the complex and multivariate daily relationship between snow and avalanches. To reconstruct the past or project the future of avalanche activity, we use the best meteorological dataset available to date for the considered spatio-temporal domain as input to our machine learning model.

2 Material and methods

2.1 Overall approach

We developed a machine learning model that predicts the daily number of avalanches per aspect sector from simulated weather and snowpack conditions, after training on a human-based avalanche inventory. Similar to Sielenou et al. (2021); Mayer et al. (2022); Viallon-Galinier et al. (2023); Hendrick et al. (2023), we relate daily meteorological and snow conditions to daily avalanche activity. We split the area under study into four aspect sectors: North, East, South, and West-oriented sectors, as snow-related variables depend heavily on incoming solar radiation and sometimes on wind direction.

In this section, we first present the target variable (i.e., the variable we want to predict): the daily number of avalanches in the Haute-Maurienne valley, as reported by the EPA inventory. Next, we describe the model input features (i.e., predictors of the target variable). These features are either derived from reanalysed weather and snow data or from snow-climate simulations. Then, we detail the machine learning model itself. In particular, we describe how we account for uncertainties and biases in the avalanche observations. Finally, we define three indicators derived from the daily prediction that we use to evaluate the model and quantify avalanche activity.

2.2 Target: daily number of avalanches per aspect sector

The Enquête Permanente sur les Avalanches (EPA) is a French survey that reports all avalanches that meet a specified threshold on designated avalanche paths (Bourova et al., 2016; Eckert et al., 2010a). Approximately 3900 avalanche paths are monitored in France, mainly by operators of the forestry service (Eckert et al., 2013). Each path has an associated flow threshold, and any avalanche that exceeds this threshold is expected to be included in the inventory. Each reported avalanche is at least described by the associated path and the time window during which it occurred. Each path is characterized by the average aspect and altitude of the typical release zone. One of the main advantages of this dataset is its presumed temporal homogeneity, as the flow thresholds have remained unchanged over the years. However, the EPA has been rigorously revised after the Montroc avalanche in 1999 and, in 2006, the observation protocol was revised, including the addition of new descriptors of the avalanches (Bourova et al., 2016). A second advantage of the EPA is that avalanches are reported almost daily. The occurrence of an avalanche is controlled by the snow conditions at the time it occurred, and is only very weakly controlled, for instance, by the snow conditions averaged over the winter season (Schweizer et al., 2003). To learn the relationship between snow conditions and avalanche occurrence, it is thus important to rely on data of high temporal resolution of the order one day. However, the EPA is not perfect in this aspect: for each avalanche event recorded, the observer identifies a time window defined by 2 d [d₁,d₂] during which the avalanche probably occurred. This period may correspond either to the interval between two successive observations or to a narrower time frame based on the observer's evaluation of meteorological and snow conditions. The length of this interval is rarely 1 d but generally spans 1 to 5 d due to avalanche danger, visibility conditions, holidays, etc. Lastly, it should be noted that the EPA dataset does not inventory all avalanches but only those in pre-designated paths. This survey was primarily developed for risk management in areas near the valley floor, so avalanches in unpopulated areas or those threatening snow recreationists are far from all included.

This study focuses on the Haute-Maurienne valley, located in the French Alps. This valley has long been a critical area for avalanche observation due to its geography, with relatively high mountains and steep avalanche paths (e.g. Kern et al., 2020; Zgheib et al., 2020; Viallon-Galinier et al., 2023; Barkat et al., 2024). We focus on the upper valley, which is characterized by more active avalanche paths. The study area comprises the three district municipalities of Bessans, Bonneval-sur-Arc and Lanslevillard (Fig. 1). It is characterized by about a hundred monitored EPA avalanche paths and thousands of avalanches have been recorded between 2006 and 2023. The typical flow threshold is at an elevation around 1800 m (slightly above the valley bottom) while the avalanche paths extend to maximum elevations of up to 3300 m, with an average upper limit of about 2700 m. In this subset, the observations show very few avalanches in the 1950s and significantly higher avalanche counts in recent decades (Fig. 2), which is inconsistent with larger-scale trends (Eckert et al., 2013). This behavior may be specific to the considered domain. The number of days of uncertainty associated with each avalanche release date is variable: the exact day a given avalanche occurred is known only in about 18 % of the cases and 96 % of the avalanche occurrence days are known with less than five days of uncertainty (Fig. 3).

Figure 1 EPA paths in the selected area of the Haute-Maurienne valley. © OpenStreetMap contributors 2025. Distributed under the Open Data Commons Open Database License (ODbL) v1.0. © Mapbox. All rights reserved. © IGN France. All rights reserved. © ONF France. All rights reserved.

Figure 2 Raw number of avalanches reported in the EPA dataset from 1950 to 2023 in the studied area. Winter seasons are denoted by the year in which they end.

Figure 3 Distribution of the number of days of uncertainty concerning the triggering date from winter season 2007 to 2023 in the studied area.

In practice, we filtered the raw avalanche events recorded in the EPA to overcome its limitations. First, we excluded the data prior to the last renovation of the observation protocol in 2006 to evaluate the machine learning model, considering only the period spanning from winter seasons 2006/2007 to 2022/2023. Second, we assumed that an event reported in the time window [d₁, d₂] occurred between 18:00 UTC on day d₁−1 and d₂ at 18:00 UTC. We also chose to exclude events with more than 5 d of uncertainty from this analysis, which correspond to 3.6 % of the listed events. In other words, we arbitrarily assumed that these very uncertain but sporadic events never happened. Third, we considered only the paths that have been monitored during the entire period 2006–2023. This final filtered dataset consists of n = 89 EPA avalanche paths and m = 2080 recorded avalanches between the winter seasons 2006/2007 and 2022/2023. Winter seasons will subsequently be denoted by the year in which they end; for example, winter 2020/2021 will be denoted as 2021.

2.3 Features: weather and snowpack

The input features for the machine learning model are generated from a series of models. An atmospheric model provides input data to a snowpack model, which simulates the evolution of snowpack stratigraphy. From this, various predictors related to avalanche formation are calculated.

2.3.1 Weather forcings

We used two types of weather forcings: reanalysed meteorological data and climate simulations. These forcings define the spatial resolution of the model chain. It is produced over elementary areas designed to represent the main drivers of spatial variability in mountains: the so-called “massifs” (here Haute-Maurienne) for which meteorological and snow conditions are estimated at discrete elevation levels spaced every 300 m (here 1800, 2100, 2400 and 2700 m), aspect classes (here N, E, S, W) and slope angles (here, 40°). The weather forcings are generated at hourly resolution.

To study past weather and snow conditions, we used the SAFRAN meteorological reanalysis, which provides meteorological data for the period 1958–2023 (Vernay et al., 2022b). This reanalysis combines information from atmospheric models (the ERA-40 reanalysis from the winter seasons 1958/1959 to 2001/2002, the model ARPEGE since winter season 2002/2003) and in situ meteorological observations.

To study climate trends, we used the ADAMONT data (Verfaillie et al., 2017), which provides climate runs based on different pairs of global climate models (GCMs, CMIP5) and regional climate models (RCMs) in Europe (EURO-CORDEX initiative) (Jacob et al., 2014). The ADAMONT data is based on quantile mapping during the historical period (from winter seasons 1951 to 2005) (Themeßl et al., 2011; Maurer and Pierce, 2014), which is a statistical bias-correction method that adjusts the distribution of climate model outputs to match that of a reference dataset, which is here the SAFRAN reanalysis. The climate projections (winter seasons 2006–2100) are available for three Representative Concentration Pathways (RCPs) (Van Vuuren et al., 2011): RCP8.5, RCP4.5 and RCP2.6. A limitation of quantile-mapping approaches is that they rely on the assumption that the statistical relationship between modeled and observed variables estimated during the reference period remains valid under future climate conditions. This assumption may be particularly problematic for extreme values, whose characteristics can evolve in a changing climate and may fall outside the range represented in the calibration dataset. The ADAMONT method partly addresses this issue through a specific treatment of extreme values designed to improve the correction of distribution tails.

Similar to Bonsoms et al. (2025), we investigated only two greenhouse gas emission scenarios: RCP4.5 and RCP8.5, excluding RCP2.6 because it is available only for a few pairs of GCM/RCM models. Overall, we used 18 GCM/RCM pairs that were adapted locally using ADAMONT quantile mapping. These pairs provide weather data over the historical period and for the two emission scenarios (see details in Table 1).

Table 1 Selected GCM/RCM model pairs and their associated winter seasons in which they are defined. These pairs provide RCP4.5 and RCP8.5 scenarios.

2.3.2 Snowpack simulations and derived avalanche formation proxies

The meteorological data are used as input for the Earth continental surface model SURFEX. The soil model is ISBA (Decharme et al., 2011). The snowpack model is the detailed model Crocus (Vionnet et al., 2012). Crocus simulates the temporal evolution of the snowpack stratigraphy, represented with up to 50 layers, continuously throughout the year. We used a 10-year spin-up to obtain realistic initial conditions at the start of the simulation period. In this study, Crocus outputs were produced at a 6 h temporal resolution. The resulting dataset includes numerous snow-related variables (numerous variables for each layer) across four elevation bands and four slope aspects, four times a day, leading to a high-dimensional feature space that must be reduced to limit model complexity and avoid overfitting (Ying, 2019).

The release of an avalanche can be formally classified into five avalanche problems (European Avalanche Warning Services (EAWS), 2022): new snow, wind slab, persistent weak layers, wet snow and gliding snow. The machine learning model must predict all avalanches, regardless of the avalanche problem. Therefore, it is essential to construct predictors that reflect this diversity.

Most predictors were selected to characterize snowpack conditions in the release area and therefore target avalanche release processes, such as wind slab formation or persistent weak-layer instability. Accordingly, variables directly related to release conditions (e.g., wind speed and direction, strength-to-stress ratio) were provided only at 2700 m, which approximately corresponds to the mean elevation of the avalanche release zones. In contrast, a subset of variables was extracted at 1800, 2100, 2400, and 2700 m because conditions along the avalanche path may influence avalanche flow and runout altitude. This distinction is particularly relevant because the dataset only includes avalanches that reach a predefined observation threshold. Although these factors are expected to play a secondary role compared with release-area conditions, they may affect whether a released avalanche is recorded in the dataset. For example, sufficient snow depth or available transportable snow at lower elevations may be required along the path for an avalanche to reach the threshold. Solid and liquid precipitation were also provided at all elevations because precipitation can affect both avalanche release and avalanche flow. Wet-snow-related variables were similarly provided at all elevations because snowpack wetting generally progresses from lower to higher elevations. Accounting for this elevation-dependent evolution allows a more complete characterization of wet-snow conditions throughout the avalanche path. We retained variables with a plausible physical link to avalanche activity rather than applying a strict pre-selection procedure. Gradient boosting methods are generally robust to the presence of weakly informative predictors, allowing physically relevant variables to be included without substantially degrading predictive performance. This approach preserves the physical interpretability of the model while allowing the algorithm to identify the most informative predictors among the provided features.

For each day, for a given aspect sector, we computed 103 features (Table 2) as follows:

• New snow (28 predictors): the primary driver of avalanche events is the amount of new snow. We characterized it by the depth of new snow accumulated over the past 24, 72, and 120 h at all the studied elevations. To quantify the intensity of solid precipitation, the mean and maximum snowfall rates over the last 24 and 48 h are also included as predictors.

• Wind slab (12 predictors): the occurrence and magnitude of snow drift are mainly driven by the wind speed and the characteristics of the snow surface (Mott et al., 2010). Wind speed must exceed a critical threshold velocity to initiate wind-blowing snow. The mean and maximum wind speed over the last 24 and 48 h are thus used as features. Wind direction is also used as a predictor to account for preferential deposition and erosion (with four dummy variables for North, West, South, and East components). In addition, snow cohesion affects its susceptibility to wind transport. The surface penetration depth of the Rammsonde was therefore included as a proxy for snow erodibility.

• Persistent weak layers (6 predictors): these problems are related to the presence of a weak layer in the pre-existing snowpack. The strength-stress ratio S_n (Roch, 1966) compares the shear strength of one layer to the shear stress due to the overlying snowpack. This ratio is used to detect slab avalanche-prone situations (e.g. McClung, 1981; Reuter et al., 2022; Viallon-Galinier et al., 2023). We used the minimum values of S_n across defined depth intervals: 0.3 m bins from 0.2 m (similar to Reuter et al. (2022)) to 1.7 m, and one additional bin for depths > 1.7 m. This predictor is computed only for the highest elevation (2700 m), as it is mainly involved in the release process.

• Wet and gliding snow (44 predictors): Mitterer and Schweizer (2013) showed that the release of wet avalanches is related to the snowpack liquid water content and the thickness of the snowpack that is isothermal and/or wet. We defined the snowpack liquid water content ratio r as the ratio of liquid water content to snow water equivalent. We used the maximum values of r on the last 24 and 72 h and its temporal evolution ($r_d-r_{d-\mathrm{24}\mathrm{h}}$, $r_d-r_{d-\mathrm{72}\mathrm{h}}$, $r_d-r_{d-\mathrm{168}\mathrm{h}}$) as predictors. We also added the liquid precipitation rate as a predictor.

• General (13 predictors): the number of avalanches is related to the number of avalanche paths in the considered domain. For instance, in Haute-Maurienne, only 14 avalanche paths face West and 36 paths face North. This number is thus added to the predictor list. Total snow depth is also an important predictor for all avalanche problems. For instance, a heavy snowfall at high elevation on a snow-free ground may not lead to a substantial avalanche flow. Similarly, the Rammsonde surface penetration depth was additionally retained as a general indicator of snow availability along the avalanche path.

Table 2 Selected predictors per aspect sector. All elevations means 1800, 2100, 2400, 2700 m. We recall that 1 d d is defined between 18:00 UTC on day d−1 and 18:00 UTC on day d.

For each days, those quantities are computed at 18:00 UTC. Note that the goal here is not to define the best predictors of avalanche formation, in general. The idea is to span different avalanche situations to achieve a sufficient predictive power with a limited number of features. In addition, the presented classification of features into avalanche problems is somewhat arbitrary, as some variables may be associated with multiple problems.

2.4 Machine learning model

2.4.1 Machine learning algorithm

Most comparable studies (Sielenou et al., 2021; Mayer et al., 2022; Viallon-Galinier et al., 2023; Hendrick et al., 2023) used the random forest algorithm. This method aggregates predictions from numerous independent decision trees. The extreme gradient boosting (XGBoost) algorithm (Chen and Guestrin, 2016) is another learning method based on decision trees. Unlike Random Forests, which build multiple independent trees in parallel, XGBoost constructs trees sequentially, where each new tree corrects the errors of the previous ones. Unlike Random Forests, XGBoost uses a loss function that is explicitly minimized through gradient boosting. This technique iteratively reduces prediction errors by fitting new trees to the residuals of previous trees. We used this model especially to benefit from a custom loss function that accounts for uncertainties in the observed release days (see Sect. 2.4.2). To predict the daily number of avalanches, we chose the regression variant of the XGBoost algorithm (implemented using the XGBoost Python package) rather than a classifier.

2.4.2 Explicit modeling of uncertainty regarding the release date

A loss function aims to compute the difference between the prediction (here the daily number of avalanches) and the observation at each step of the training process. We defined a custom loss function to account for the uncertainty in the release date of each observed avalanche. For each day, we calculate the largest time period associated with all the observed avalanches that could occur on that day. The number of avalanches predicted by the model during this period is then compared to the range of the number of avalanches (the exact number is often not known) that could have occurred during the same period. This process is illustrated in Fig. 4. In more detail, the loss function for day d is built as follows:

• We identify all the avalanches in the EPA dataset that could have occurred on day d. Each avalanche i is characterized by an uncertainty interval defined by two dates d_i,1 and d_i,2. We define ΔT = $[d_{min},d_{max}]$ as the interval spanning from the earliest d_min = ${min}_i\left(d_{i,\mathrm{1}}\right)$ to the latest d_max = ${max}_i\left(d_{i,\mathrm{2}}\right)$. If no avalanche is observed on day d, then ΔT consists only of the single day d. For instance, for d = $\mathrm{25}/\mathrm{12}$ in Fig. 4, ΔT = $[\mathrm{24}/\mathrm{12},\mathrm{29}/\mathrm{12}]$.

• We compute the range $[n_{min},n_{max}]$ of the number of observed avalanches within the time window ΔT. For example, in Fig. 4, for ΔT = $[\mathrm{24}/\mathrm{12},\mathrm{29}/\mathrm{12}]$, n_min = n_max = 7 avalanches. Similarly, for ΔT = $[\mathrm{13}/\mathrm{12},\mathrm{16}/\mathrm{12}]$, n_min = 2 and n_max = 4, depending on the triggering day of the avalanches associated with the paths 73144202 and 73040020.

• We compute the sum of the daily predicted number of avalanches over ΔT.

• We compute the signed distance between this predicted value and the interval $[n_{min},n_{max}]$.

• This distance is then distributed uniformly across all days within ΔT by dividing it by the duration of the interval (in days). The resulting quantity defines the daily residual R used as the loss.

Figure 4 Illustration of the custom loss-function used in the training process, taking the example of the avalanches recorded in the Southern aspect sectors between 10 and 29 December 2019. (a) Avalanches reported in the EPA dataset. Each continuous line corresponds to an avalanche, and each point represents a potential day for one avalanche. (b) Time-period ΔT corresponding to each date. (c) Minimum and maximum number of avalanches observed during the time-period ΔT associated to each date.

The details of the chosen XGBoost hyper-parameters are provided in Appendix A for reproducibility.

2.4.3 Derived indicators of avalanche activity

The model predicts the daily number of avalanches for each aspect sector. To investigate trends in avalanche activity, we defined three indicators that quantify the avalanche activity at different time scales over the whole domain (i.e., we sum the number of simulated avalanches for all four aspect sectors):

• the number of avalanches per winter season (September to May).

• the number of avalanches per month.

• the maximum number of avalanches that occur during 1 rolling week for each year. This indicator reflects the magnitude of the largest avalanche cycle each year, assuming that large cycles are associated with a high number of avalanches. To investigate changes in high values of this indicator, we fit non-stationary Generalized Extreme Value (GEV) distributions (Coles, 2001) on the data, and investigate the 30-year return level (i.e., the level associated with a 30-year return period). The application of this framework is detailed in Appendix B.

To assess temporal trends in these indicators and to fit the GEV distributions, we employed the Bayesian statistical framework implemented in the PyMC Python library (Abril-Pla et al., 2023). This approach enables us to incorporate prior knowledge and compute the full posterior distributions of the parameters, including their uncertainties, which proved useful for many related avalanche problems (e.g., Eckert et al., 2010c; Fischer et al., 2020).

2.4.4 Model evaluation

The model is evaluated on the mean absolute errors (MAE) of these three indicators during the winter seasons 2006/2007 to 2022/2023. As in Viallon-Galinier et al. (2023), we used a “leave-one-year-out” (LOYO) evaluation method: each winter of the studied period is, in turn, treated as a test set, with a model trained on the remaining winters.

To assess the benefits of our new method, we trained four versions of the model.

• $M_{d_{\mathrm{1}}-d_{\mathrm{2}}}^{\mathrm{2007}-\mathrm{2023}}$ was trained on the period defined between the winter seasons 2007 to 2023 and accounts for the uncertainties in the triggering date by using the custom-loss function described in Sect. 2.4.2.

• $M_{d_{\mathrm{1}}-d_{\mathrm{2}}}^{\mathrm{1961}-\mathrm{2023}}$ was trained on the period defined between the winter seasons 1961 to 2023 and accounts for the uncertainties in the triggering date.

• $M_{d_{\mathrm{2}}}^{\mathrm{2007}-\mathrm{2023}}$ was trained on the period defined between the winter seasons 2007 to 2023 and uses the default loss function of XGBoost and assumes that each avalanche occurred at the end of its period of uncertainty.

• $M_{d_{\mathrm{2}}}^{\mathrm{1961}-\mathrm{2023}}$ was trained on the period defined between the winter seasons 1961 to 2023 and uses the default loss function of XGBoost.

Using a similar approach, Viallon-Galinier et al. (2023) trained their models on the 1960–2018 period, assuming that each avalanche occurred at the end of its period of uncertainty. In the four models investigated here, we used the same features computed from the S2M reanalysis and the models were evaluated on the same winter seasons from 2006/2007 to 2022/2023.

3 Results

In the following sections, the four versions of the model are first evaluated over the winter seasons from 2006/2007 to 2022/2023. Then, the best one is used to correct the historical trend in avalanche activity between winter seasons 1959 to 2023 based on the S2M reanalysis. Finally, we compute 1950–2100 trends from climate simulations.

3.1 Model evaluation between 2007 and 2023

We first evaluate the four versions of the model between winter seasons spanning from 2007 to 2023 (Table 3). Accounting for both observation uncertainties and potential observation bias before the 2007 winter season reduces prediction errors. The annual MAE is reduced by half (from about 50 avalanches to 28 avalanches per year), the seasonal MAE by a factor of 2.5 (from about 5.5 avalanches to 2.3 per month) and the worst week MAE by half to one fourth (from 10–19 avalanches to 5 avalanches per week). Thus, the relevance of the modeling choices is confirmed, as they substantially improve the methodology and lead to the best scores on the three metrics.

Table 3 Scores for the four models. Training periods are expressed as winter seasons denoted by the year in which they end (e.g. 2007–2023 corresponds to the winter seasons spanning from 2006/2007 to 2022/2023). The mean absolute errors correspond to a number of avalanches. The values in brackets are percentages calculated from the average value of each indicator. The model in bold is the one selected for the rest of this study.

The predictions of the model $M_{d_{\mathrm{1}}-d_{\mathrm{2}}}^{\mathrm{2007}-\mathrm{2023}}$ are detailed in Fig. 5. We recall that these predictions are generated with the leave-one-year-out method (see Sect. 2.4.4). The inter-annual variability of the number of avalanches per year (between about 40 and 500) is well predicted by the machine learning model (Fig. 5a, b). The model correctly identifies specific years, such as winter season 2018 with a record-breaking number of avalanches, and 2011 with very low avalanche activity. The mean absolute annual error is 32.5 avalanches per year (26.5 % of the mean number of avalanches per winter season). The bias is small: 8.5 avalanches (7 % of the mean number of avalanches per winter), showing that the error is almost centered on zero. The model is also capable of reproducing the average number of avalanches per winter month (Fig. 5c, d). There are almost no avalanches predicted or reported before November. Then the number of avalanches progressively rises until the peak activity in March. Avalanche activity almost disappears in May. The mean absolute monthly error is approximately 2.3 avalanches (17 % of the mean number of avalanches per month during the winter season), with a small bias of 0.9 avalanches. Focusing on the maximum number of avalanches during 1 rolling week per season (Fig. 5e, f), the mean absolute error is approximately 4.5 avalanches, representing 12.5 % of the mean value. The model correctly predicts the inter-annual variability, but slightly underestimates the maximum number of avalanches per week, by around 10 % over the entire period, as shown by a negative bias of around −4.0 avalanches for the most active week of the year. All in all, the results show that the model performs well at annual and monthly and at the most active week scales. This demonstrates that the model is suitable for investigating climate-driven trends in these avalanche activity indicators.

Figure 5 Evaluation of the model $M_{d_{\mathrm{1}}-d_{\mathrm{2}}}^{\mathrm{2007}-\mathrm{2023}}$ using the LOYO method on the annual number of avalanches (a, b), the mean number of avalanches per month over the 17 winter seasons (c, d), and the maximum number of avalanches during a week per winter season (e, f). The MAE and the bias are computed using the distance to the interval associated with the uncertainties (error bars on observations).

3.2 1958–2023 reanalysis

A single XGBoost model is trained using all the winter seasons from 2007 to 2023. The model is then applied using the S2M reanalysis dataset from winter seasons 1959 to 2006 as input data. To obtain a continuous time series spanning winter seasons 1959 to 2023, these reconstructed values were combined with the predictions for the winter seasons between 2007 and 2023 obtained using the leave-one-year-out method described above. The predicted daily number of avalanches is used to derive the three previously defined indicators of avalanche activity (Fig. 6). In addition, we compute the temporal linear trends of these indicators using a Bayesian framework (Abril-Pla et al., 2023). The linear trends are systematically expressed in % per decade relative to the 1975–2005 climatology, with its median value and the 10th–90th percentile posterior credible interval.

The predicted number of avalanches per winter season shows large inter-annual variability but an average decreasing trend (Fig. 6a). For instance, years with more than 200 avalanches are rather frequent (approximately every three years on average) before the late 1980s and become rare thereafter. As expected, the trends from the raw observations are completely off from these assessments, which highlights that raw observations over a long period may be affected by sampling biases. Quantitatively, our model predicts a median decreasing linear trend of −4.8 % per decade in the annual number of avalanches over the period 1958–2023, with a 10th–90th percentile credible interval ranging from −10.0 % to +0.4 % (Appendix Fig. A1).

The decrease of the annual number of avalanches is not uniform during the winter season as shown by contrasting 30-year periods at the beginning and at the end of the dataset: 1958–1988 and 1993–2023 (Fig. 6b). For instance, the resulting trends indicate a median slope of −1.6 % (with a credible interval ranging from −7.7 % to +4.6 %) per decade in the avalanche activity in December, January, and February (DJF) and of −8.3 % (between −13.4 % and −3.5 %) per decade in March, April and May, with a peak at −14.3 % for the April avalanche activity (Appendix Fig. A1). Therefore, the decrease in avalanche activity in the past has been greater in spring than in winter.

The weekly maximum number of avalanches for each winter season is shown in Fig. 6c. As we use a Bayesian framework to fit the GEV distribution on the data, the slope distribution of the return levels is directly provided within the inference procedure. The results show a decreasing slope of the 30-year return level by −1.4 % (between −7.3 % and +4.3 %) per decade relative to the mean return level between 1975 and 2005. This downward trend is less significant than the one observed in the total number of avalanches, suggesting that the magnitude of large cycles has decreased at a lower rate than the total number of avalanches.

Figure 6 Past avalanche activity: (a) annual number of avalanches, (b) mean number of avalanches per month, (c) maximum number of avalanches per week. The model predictions are shown in orange. The raw observations are shown in blue. For panels (a) and (c), linear trends are represented with their 10, 25, 50, 75, 90th percentiles.

3.3 1950–2100 climate simulations

The same machine learning model is applied to all members of the climate simulations (18 per scenario), and the avalanche activity is quantified by the three previously defined indicators. In contrast to the single reanalysis, the linear trends derived from climate simulations for a given scenario are not only made uncertain by the inter-annual variability but also by the inter-model variability. To compute the probability distribution function of the trend slope, we arbitrarily assumed that all GCM/RCM pairs are equiprobable and evaluated the overall credibility interval by pooling posterior samples derived from individual projections (see Appendix Fig. A1).

The annual number of avalanches clearly shows a decreasing trend under advanced climate change and its consequences for snowpack characteristics (Fig. 7a). The mean number of avalanches in the historical simulations is around 25 % higher than that in the reanalysis. However, the trends are similar with a historical decrease at a median rate of −5.9 % per decade, but with large uncertainties (10th–90th percentile posterior credible interval ranging from −17.1 % to 7.6 %). The projections for the 21st century show a decrease in the annual number of avalanches at median rates of −4.6 % per decade for the RCP4.5 scenario and −8.8 % for the RCP8.5 scenario with respect to the period spanning from the winter seasons 1975 to 2005. The associated uncertainties are shown in Table 4. In other words, the results indicate a reduction of approximately 37 % (respectively 65 %) in the mean number of avalanches over the 2070–2100 period compared with the 1975–2005 baseline in the RCP4.5 scenario (resp. RCP8.5).

Figure 7 Climate simulations of avalanche activity: (a) 30-year running mean of the annual number of avalanches, (b) mean number of avalanches per month over different periods and (c) 30-year return level of the maximum number of avalanches per week.

The expected decline of avalanche activity is heterogeneous between seasons (Fig. 7b). First, the higher number of avalanches predicted by historical climate simulations compared to the reanalysis mainly arises from stronger activity in March and April, while DJF avalanche activity is broadly consistent between the two datasets. Then, the projected DJF avalanche activity appears almost unaffected by climate change, except for the RCP8.5 at the end of the century. Conversely, a marked decline in avalanche activity is expected in March, April, and May (MAM) for all scenarios. The projected decrease in the annual number of avalanches is thus primarily driven by the strong reduction in spring avalanche activity. For instance, the number of avalanches in DJF is expected to decrease only by 3.1 % per decade under RCP4.5 and by 6.9 % per decade under RCP8.5, while in MAM, the expected decreases correspond to 5.4 % per decade under RCP4.5 and 10 % per decade under RCP8.5, relative to the 1975–2005 reference period (Table 4). April was historically one of the most active months for avalanches; however, under the RCP8.5 scenario, by the end of the century, avalanche activity in April is projected to be very limited. Similarly, in May, avalanche activity is expected to be nearly absent under both scenarios, indicating a substantial shortening of the avalanche season.

Table 4 Median values and 10th–90th percentile posterior credible interval associated with the slope of the avalanche activity indicators in the future in % per decade with respect to the 1975–2005 period. The associated ensemble distributions are represented in Appendix Fig. A1.

It is projected that there will be a decrease in the amplitude of large avalanche cycles by the end of the century, with a lower rate for RCP4.5 and a higher rate for RCP8.5 (Fig. 7c). Yet, considering the credible intervals, for the near future (2020–2050) or under the RCP4.5 scenario at all temporal horizons, large avalanche events with intensities comparable to those observed during the 1975–2005 reference period remain possible with a similar probability. Indeed, the 30-year return levels show respective trends of −2.2 % per decade and −5.5 % per decade for the RCP4.5 and RCP8.5 scenarios relative to the mean return levels computed from the historical runs between winters from 1975 to 2005. In other words, the 30-year avalanche cycle is expected to be composed of approximately 15 % (resp. 35 %) fewer avalanches in 2100 compared to the 1975–2005 reference period in the RCP4.5 scenario (resp. RCP8.5).

4 Discussion

4.1 Methodology

This work uses a new methodology compared to previous studies to investigate future changes in avalanche activity by using a regressor and considering two key aspects of the training dataset: the heterogeneity of observations over time and uncertainties associated with triggering dates.

First, using a machine learning regressor algorithm enables us to progress beyond studies that used binary classifiers. To our knowledge, all studies that predict daily avalanche activity from snow stratigraphy using machine learning methods rely on classification algorithms to predict whether a day is avalanche-prone (Mayer et al., 2024; Viallon-Galinier et al., 2023; Hendrick et al., 2023; Kronholm et al., 2006), or to categorize daily avalanche activity into three classes: null, moderate, or high (Sielenou et al., 2021). These models do not provide quantified information about the magnitude of avalanche activity during the avalanche days, and are mainly used for short-term avalanche hazard forecasting and for pre-detecting situations that may require specific attention by the forecaster. Their predictions generally show high recall but low precision: most observed avalanche days are well predicted, but the predictions may contain many false alarms. Using a regressor thus allows us to predict the daily number of avalanches and to quantify climate trends in avalanche activity using different indicators, without treating all avalanche days equally. However, the number of avalanches gives no information on the size of avalanches and the potential damage: many small avalanches may have no impact at all, while a small number of large ones could cause significant damage. Future research could investigate additional target variables, such as deposit volume (Kern et al., 2020), to enhance the scope of potential applications. In addition, we may improve our regressor algorithm. Indeed, the observed daily avalanche counts follow a zero-inflated count distribution, characterized by integer-valued observations with an excess of zeros (Young et al., 2022). We did not use a specific model to handle this distribution. More complex models adapted to zero-inflated count data (Fávero et al., 2024), such as a two-stage model that combines a classifier and a Poisson regressor, might improve the methodology (e.g. Orhobor et al., 2023; Rožanec et al., 2023). In such a framework, the GEV component would need to be formulated using a discrete extreme value distribution, as proposed by Evin et al. (2021). In the present study, it relies on a continuous GEV formulation as the model predicts continuous daily numbers of avalanches. Moreover, the model tends to underestimate the most active week of the year (Fig. 5e), indicating limitations in its ability to predict extreme avalanche activity. These underestimations likely result from the combined effects of the rarity of extreme avalanche cycles in the training dataset, the smoothing introduced by the loss function that distributes avalanche occurrences across potential release days, and the spatial averaging of meteorological and snowpack conditions over the SAFRAN geometry, which can attenuate the predictor signals associated with the most intense events. Combining a similar machine learning framework with an extreme value approach could help better represent the tail behavior of the distribution and improve the credibility of predicted extremes (Evin et al., 2021).

Secondly, to address heterogeneity in past observations, we trained a machine learning model using the EPA dataset between the winter seasons 2006/2007 and 2022/2023. We have shown that training the model on a larger time period (Table 3), as done in previous studies (e.g. Viallon-Galinier et al., 2023; Sielenou et al., 2021), leads to lower scores. Besides, the predicted number of avalanches based on the reanalysis significantly deviates from the observations in the distant past (Fig. 6). These results suggest that avalanche reporting was not consistent throughout the 1958–2023 period in the Haute-Maurienne valley, with possible under-reporting during the earlier decades. While this result is specific to the studied area, it confirms that the EPA dataset, as with most avalanche records stemming from human observations or archives, should be used with caution when considering the distant past (Giacona et al., 2021). However, training the model on a smaller dataset also has its drawbacks. In particular, the approach implicitly assumes that the range of snow conditions encountered during the model application period (1950–2100) is represented within the training dataset (2006–2023). This assumption becomes more plausible as the length and diversity of the training period increase. When the model is applied to conditions outside the range of its training data, its predictions are likely to be unreliable. In other words, the model cannot handle new, unknown avalanche situations and only predicts the future probability of situations that have already occurred. This limitation is particularly critical when analyzing extreme avalanche cycles. In this regard, including the 2017/2018 winter season – characterized by exceptionally unusual snow conditions – within the training dataset is especially valuable. To overcome this limitation, physics-based avalanche models (Durand et al., 1999; Reuter et al., 2022) may help this machine learning approach handle new conditions not encountered during training.

Thirdly, we showed that accounting for the uncertainties in the release dates in the training step significantly improves the prediction scores (Table 3). This was achieved through the use of a custom loss function in the XGBoost algorithm. This approach enabled us to accurately associate daily snow conditions with the corresponding avalanche activity without making any assumptions about the exact triggering dates. Standard random forests widely used in the literature (e.g. Mayer et al., 2022; Viallon-Galinier et al., 2023; Sielenou et al., 2021) do not allow for the integration of uncertainties through a loss function.

4.2 Input data

The results are strongly influenced by the quality of the input data, specifically the avalanche observations and associated predictors. Enhancing the quality of these data sources may lead to improved outcomes. Furthermore, as this study is limited to a single location and utilizes high-quality data, additional research is necessary to determine whether these findings can be generalized to other regions.

The avalanche activity described in this paper corresponds to the one reported by the EPA, which is not representative of all the avalanches that can occur, particularly those at higher altitudes that do not reach the EPA thresholds (Fig. 1). Nevertheless, these avalanches can still pose a threat to snow recreationists. In addition, the EPA thresholds of the studied area are located within a relatively narrow elevation range, and we assumed that avalanche activity is driven by snow conditions aggregated over all the elevation bands between 1800 and 2700 m. As a result, the influence of elevation on temporal trends in avalanche activity cannot be explicitly isolated in the present analysis. This limitation is noteworthy, as elevation has been shown to play a major role in governing snowpack evolution and long-term changes in avalanche activity, particularly in the context of climate change (Eckert et al., 2024). Moreover, this study is specific to the Haute-Maurienne valley, characterized by numerous recorded avalanches. The exact same methodology is not applicable in other areas where the too small number of recorded avalanches does not yield a sufficient training data set. Working on model transferability (Wang et al., 2022) or by reasoning at alternative spatial scales (e.g. considering individual avalanche paths and accounting for their topographical specificity by adding topographic properties in the features) could further strengthen the methodology and facilitate its application to other locations.

The model is trained on weather and snow predictors derived from the S2M reanalysis at the massif scale. In the Haute-Maurienne valley, snowfall is highly influenced by specific meteorological events (easterly return flows) that can produce substantial localized accumulation near the Italian border while nearby areas receive minimal snowfall, as exemplified by the large 2008 avalanche cycle in Southeast France (Eckert et al., 2010b). Consequently, the snowpack simulations cannot fully capture this heterogeneity and may yield locally unrepresentative conditions. Nevertheless, as these issues are present in both the training and application datasets, there should still be a predictive relationship between simulated meteorological and snowpack conditions and avalanche activity. Besides, the S2M reanalysis is also subjected to biases that can vary with time, which affects reanalysed climatic trends (Vernay et al., 2022b): the number of assimilated observations varies over time and the meteorological conditions are simulated by the ERA40 model (Uppala et al., 2005) before 2002 and by the operational forecasts of the French global NWP model ARPEGE thereafter. In addition, we assumed here that avalanche hazard evolution is solely controlled by changes in snow cover. However, avalanche activity can also be impacted by changes in land use, such as land abandonment and reforestation (Mainieri et al., 2020; Zgheib et al., 2022; Moen et al., 2004). The projected reduction of avalanches may also initiate a positive feedback loop: fewer avalanches may increase the density of small-diameter trees in avalanche starting zones, reducing future avalanche magnitudes even more (Teich et al., 2012). The choice of predictors could also be further refined. For example, we only used the strength-to-stress ratio to characterize persistent weak-layer avalanche problems. It provides a very simplified description of snowpack stability. More physically based metrics, such as the critical crack length (Gaume et al., 2017), could potentially offer a more realistic representation of avalanche release processes and improve predictive skill. However, the identification, evaluation, and optimization of relevant predictors constitute a substantial research topic in their own right and are beyond the scope of the present study. The objective here was not to determine the optimal predictor set, but rather to assess the ability of a physically informed machine-learning framework to reproduce and project avalanche activity.

The climate trends computed in this study are based on the ADAMONT framework, which relies on a projection onto a massif-elevation grid of RCM simulations and on an advanced quantile mapping with the reanalysis on the historical period (Verfaillie et al., 2017). The EURO-CORDEX Regional Climate Model Ensemble simulations are generally too cold and too wet (Matiu et al., 2024; Vautard et al., 2021; Smiatek et al., 2016). The ADAMONT method is designed to reduce systematic discrepancies between climate model outputs and the reference reanalysis. However, in high-mountain environments, the correspondence between regional climate model grid points and the associated reanalysis points is challenging due to strong elevation gradients and complex topography. Furthermore, the ADAMONT bias-correction procedure is applied independently to each elevation band, without explicitly enforcing vertical coherence among corrected variables. This approach may therefore disrupt altitude-dependent physical relationships and introduce inconsistencies across elevations. These limitations are further amplified near the rain–snow transition zone, where strong non-linear responses to small temperature biases can generate significant biases in snow-related variables. Indeed, in the Haute-Maurienne valley at 2700 m, the mean winter snow depth is still around 30 % higher in ADAMONT climate simulations than in the reanalysis (not shown), which may explain why the avalanche activity computed from historical ADAMONT simulations is significantly higher compared to one computed from the reanalysis (Fig. 7). A strong assumption of quantile mapping is also the temporal stationarity of the bias, which remains questionable under pronounced climate change (Verfaillie et al., 2017). As a result, we based our analysis on climatic trends (in % per decade) rather than absolute values, which are more prone to systematic biases. Moreover, Cannon et al. (2015) analyzed the impact of bias correction via quantile mapping on extreme precipitation and concluded that quantile mapping can inflate the magnitude of relative trends in precipitation extremes with respect to the raw dataset. A specific treatment of extreme values is used in ADAMONT: for RCM values greater than the 99.5 % quantile, the quantile mapping method is not applied and a constant adjustment is applied to allow for new extremes (Verfaillie et al., 2017). This treatment of the extreme values might be better than a simple quantile mapping method, but may be insufficient to capture the correct trends in the extreme values. Future work should therefore focus on more advanced extreme-value frameworks to better constrain the upper tail and assess the robustness of the projected changes in the extreme events.

4.3 Climate trends in avalanche activity

The past trend in avalanche activity based on the S2M reanalysis (1958–2023) shows a decrease in the mean annual number of avalanches between 1958 and 2023, at a rate of approximately 5 % per decade. Based on the S2M reanalysis and the raw EPA dataset without any prior statistical pre-processing, Castebrunet et al. (2012) found no significant trend in avalanche occurrences at the scale of the Northern French Alps between 1958 and 2009. Using a more advanced statistical treatment of the EPA dataset, Eckert et al. (2013) found a decrease of 19 % in the mean number of avalanches per winter between 1980 and 2009 in the entire French Alps, which corresponds to approximately 6 % per decade, in line with our result. We also showed that the decrease in avalanche activity is more pronounced during spring, which is consistent with Reuter et al. (2025) who reported an advancement of about 3 weeks of the onset date of wet-snow activity from 1958 to 2020. Moreover, we found that the magnitude of large avalanche cycles, quantified by the maximum number of avalanches during 1 week, decreased by around 1.4 % per decade, at a lower rate compared to the annual number of avalanches. Peitzsch et al. (2021) used dendrochronology analysis in the Rocky Mountains and demonstrated that the probability of large avalanche occurrence decreased by approximately 2 % per decade from 1950 to 2017. Although this metric is not directly comparable to our results, since we analyze avalanche counts during one week rather than the occurrence of large avalanches, both metrics address changes in the frequency of large avalanche cycles over time.

Regarding future climate simulations, our results are in agreement with the small number of existing studies that provide quantitative assessments, which consistently report a marked decrease in the mean annual avalanche occurrence. Mayer et al. (2024) reported a decrease in the number of avalanche days between the periods 1990–2020 and 2069–2099, with the magnitude of the decrease depending on station elevation. At the station with the lowest elevation (1824 m, similar to the Haute-Maurienne valley floor), the number of avalanche days declined by approximately 30 % under the RCP4.5 scenario and by about 60 % under the RCP8.5 scenario. At higher altitudes, between 2300 and 2800 m, this decline is less pronounced, reaching around 10 % in the RCP4.5 scenario and 20 % to 30 % in the RCP8.5 scenario. Mayer et al. (2024) also showed that seasonality in dry-snow avalanche activity (characterized by peaks in January and February) remains mostly unchanged while the wet snow avalanches are expected to appear earlier in the season and decrease in number in April and May. Taking into account all avalanche problems, it leads to a significant reduction in spring avalanches. Similarly, in the Rocky Mountains, the wet avalanches are expected to start earlier than historical averages (around 40 to 45 d earlier in the high emission scenario) (Lazar and Williams, 2008). This is consistent with our findings which indicate a drop in the number of avalanches in April and May, especially in the RCP8.5 scenario. More generally, the observed decrease in the annual number of avalanches appears to be largely associated with the shortening of the avalanche season, as discussed in Kruyt et al. (2023). We also found that the return levels of large avalanche cycles are expected to decrease, but at a slower rate than the mean avalanche activity. This is consistent with the results of Ortner et al. (2025), which show that despite climate change, some locations continue to be at high risk of avalanches. Furthermore, Bonsoms et al. (2025) demonstrated that the mean daily snowfall in the Pyrenees at 2500–3000 m will decline significantly and emphasize that extreme snowfall events will change only slightly. Similarly, Le Roux et al. (2023) demonstrated that the expected decrease in the 100-year return level of daily snowfall is slower than the expected decrease in the mean annual maximum value. This illustrates that a decrease in mean conditions does not necessarily imply a proportional decrease in extreme values.

The snowpack conditions in the Haute-Maurienne valley are rather specific: a relatively high valley floor compared to typical alpine valleys and substantial precipitation associated with eastern flows. Therefore, the spatial generalization of these results should be addressed in future work. In particular, future analyses should explicitly account for the influence of elevation on avalanche activity, given its major role in modulating the impacts of climate change on snowpack conditions (Dumont et al., 2025) and extreme snowfall events (Le Roux et al., 2023). Based on documented changes at low elevations (Giacona et al., 2021) and high elevations (Ballesteros-Cánovas et al., 2018), a conceptual framework describing elevation-dependent avalanche activity has been proposed by Eckert et al. (2024). However, further investigations are required to provide quantitative assessments of these elevation-dependent effects.

5 Conclusions

This study aims to better understand the impact of climate change on avalanche activity in the European Alps using an exemplary alpine valley. This study brings methodological improvements and provides quantified trends for several indicators of avalanche activity.

This study focuses on the upper part of the Haute-Maurienne valley in the French Alps, a region offering several advantages for the analysis of avalanche activity. In particular, it benefits from a long-term avalanche observation dataset documenting numerous events on well-identified avalanche paths over several decades, collected using a methodology that is as homogeneous as possible over time. To model avalanche activity, we combined meteorological and snowpack modeling, using the SAFRAN reanalysis for past conditions and an ensemble of climate simulations for future projections, together with the Crocus snowpack model. This framework provides consistent meteorological and snowpack variables at the massif scale, which are used as predictors of avalanche activity.

This study introduces a machine learning regressor to predict the daily number of avalanches, going beyond previous works based on binary classification. Unlike classifiers, it provides insight into the magnitude of avalanche activity and enables the use of quantitative indicators for trend analysis. The model was trained only on post-2006 data to avoid biases from inconsistent historical observation records, and the training step accounts for uncertainties in the reported avalanche release date, which significantly reduces prediction biases and errors. The proposed chain model is generic and can be transferred to other regions where simulated meteorological and snowpack data, together with systematic avalanche observations, are available.

This study shows a clear long-term decline in avalanche activity linked to climate warming. In the past, both reanalysis and historical climate simulations indicate a clear downward trend in the annual number of avalanches. The reanalysis exhibited a trend of −4.8 % per decade over the 1958–2023 period (relative to the 1975–2005 baseline). In the historical climate simulations (1950–2005), the median trend is estimated at −5.9 % per decade. This decline is seasonally heterogeneous, mainly driven by a sharp reduction in spring avalanches (March–May), while winter activity (December–February) remains relatively stable. The 30-year return level of the weekly maximum number of avalanches shows a decline of −1.4 % per decade for the reanalysis from winter seasons 1959 to 2023. This suggests a moderate decrease in extreme avalanche events, although less pronounced than the decline in total avalanche numbers.

For the future, we used an ensemble of 18 GCM/RCM couples for the RCP4.5 and 8.5 emission scenarios. The projected evolution of avalanche activity over the historical period (1950–2005) is consistent with the decline observed between 1958 and 2005. Under RCP4.5 and RCP8.5 scenarios, the annual number of avalanches is expected to decrease by 4.6 % per decade and 8.8 % per decade, respectively (relative to the 1975–2005 baseline). By the end of the century, this corresponds to a 37 % reduction under RCP4.5 and 65 % under RCP8.5 compared to 2025. Again, while winter activity remains relatively stable, spring avalanches (March–May) are projected to nearly disappear, leading to a substantial shortening of the avalanche season. Extreme events are also expected to decrease, but at more moderate rates, with return levels expected to decrease by 2.2 % per decade and 5.5 % per decade under RCP4.5 and RCP8.5 scenarios relatively to the 1975–2005 winter seasons.

Our results show a marked decrease in avalanche activity in the Haute-Maurienne valley, both in terms of frequency and intensity, particularly driven by a sharp reduction in spring avalanches. However, unlike lower-elevation ranges such as the Vosges (Giacona et al., 2017) where avalanche activity has nearly disappeared in the 20th century due to natural and anthropogenic warming, avalanches remain a frequent and significant phenomenon in the Haute-Maurienne. This highlights the persistence of avalanche risk in high alpine environments, even under ongoing climate change. Future work should explore how evolving snow and weather regimes may continue to reshape the spatial and temporal patterns of avalanche activity in mountainous regions.

Appendix A: Hyper-parameters

The gradient boosting model's optimal hyper-parameters depend on its predictors. The choice of hyper-parameters was informed by expert knowledge. For instance, the maximum depth of the trees was set to 6, the default value that prevents overfitting while maintaining sufficient model complexity. The subsample and colsample by tree parameters have been set to 0.7 and 0.8 respectively, to prevent overfitting. Regularization tends to reduce the predicted probability of rare events. Setting lambda to zero allows the model to predict rare events. They are summarized in Table A1. Furthermore, the results show that the scores (defined in Sect. 2.4.4) are not greatly affected by the chosen hyper-parameters, demonstrating that the model is robust.

Table A1 Hyper-parameters values chosen in the machine learning model.

Figure A1 Distribution of the climatic trends based on ADAMONT simulations. Panel (a) represents the distribution of the slope associated with the number of avalanches per year. Panels (b) and (c) represent the distribution of the slope of the number of avalanches in December, January, and February (DJF) and March, April, and May (MAM). Panel (d) shows the distribution of the slope relative to the 30-year return level corresponding to the maximum number of avalanches during 1 week.

Appendix B: Generalized Extreme Value distribution

According to extreme value theory (Coles, 2001), the Generalized Extreme Value (GEV) distribution is the suitable model for block maxima (e.g., annual maximum rainfall or here the annual maximum number of avalanches considered as a continuous variable). For example, if Y is a maximum over a block of one year, the probability that Y is below y is given by the following distribution:

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}& P(Y<y)\\ & =\left\{\begin{array}{ll}\exp \left\{-{\left[\mathrm{1}+\mathit{\xi }\left(\frac{y-\mathit{\mu }}{\mathit{\sigma }}\right)\right]}^{-\mathrm{1}/\mathit{\xi }}\right\},& \text{for}\mathit{\xi }\ne \mathrm{0},\\ \exp \left\{-\exp \left(-\frac{y-\mathit{\mu }}{\mathit{\sigma }}\right)\right\},& \text{for}\mathit{\xi }=\mathrm{0}.\end{array}\right.\end{array}\end{array}$$

where:

• μ is the location parameter, determining the center of the distribution.

• σ is the scale parameter. This positive number controls the spread of the distribution.

• ξ is the shape parameter, governing the tail behavior.

This distribution cannot be applied directly to non-stationary processes, such as snow-related variables in a changing climate. To model such processes, a first approach consists in considering smaller time-windows over which the stationarity assumptions may hold, and compare the results over the different sub-periods. However this approach reduces the data quantity used for each GEV fit, resulting in significant uncertainties. We therefore chose another approach, which consists in defining the GEV parameters as functions of time. As done in Le Roux et al. (2022), we assumed that μ and σ vary linearly with time, while ξ is kept constant.