2.1 Avalanche data

Throughout the Swiss Alps, about 80 trained observers report avalanches in their area to the national avalanche warning service at the WSL Institute for Snow and Avalanche Research SLF, which is responsible for issuing the avalanche warnings for all regions in the Swiss Alps and Jura mountains. Observations are made on a daily basis during the forecasting season (typically from December to April). Avalanche observations include, among other parameters, an approximate location and the date of the avalanche release - and sometimes the time, the elevation and slope aspect of the release area, the LWC of the snow in the release area of the avalanche (dry or wet) (SLF, 2020), and the size of the avalanche, which is estimated according to the European avalanche size classification ranging from 1 (small) to 5 (extremely large) (Table 1; EAWS, 2022b). In most cases, the exact time of release is unknown and estimated, as the actual avalanche release was not observed. Avalanches can be reported as an individual avalanche event, or by reporting several avalanches in an avalanche summary report. If no recent avalanches were observed, this is also reported.

Table 1. Avalanche size s according to EAWS (2022b) and weight w used for the calculation of the avalanche activity index (Eqn. (2))

The avalanche data base includes avalanche records since 2001. During this time, the reporting system changed once. Initially, avalanche records used to be linked to the place of observation rather than the actual location of the avalanche release (IFKIS platform until summer 2019; Bründl and others, Reference Bründl2004). Today the reporting system is map-based (SLFPro, since Oct 2019). A further important difference between previous and present recording standards concerns the way avalanche summaries are reported. Before October 2019, the summary form permitted to report both dry- and wet-snow avalanches in the same report. As a consequence, in case of a summary report, it is not possible to unambiguously extract the number of avalanches of a certain size and liquid water content, nor the aspects and elevation ranges where these avalanches released from the data base (see also Mitterer and Schweizer, Reference Mitterer and Schweizer2013).

For this study, we considered avalanches that released between 1 December and 30 April.

2.2 Weather and snowpack data

We built the training and test dataset based on data from AWS used by the operational avalanche warning service in Switzerland (Section 2.2.1). We also investigated the influence of using NWP output as input data for wet-snow avalanche predictions. NWP data are briefly described in Section 2.2.2. Further details about AWS and NWP setups are provided in Pérez-Guillén and others (Reference Pérez-Guillén2022).

2.2.1 Measurement-based SNOWPACK simulations (‘nowcast’)

Snowpack simulations for the purpose of avalanche forecasting for the Swiss Alps use a station-based setup (Lehning and others, Reference Lehning1999; Morin and others, Reference Morin2020). Currently, the network of automatic weather stations (AWS) of the Intercantonal Measurement and Information System (IMIS) network (Fig. 1) consists of 186 AWS, located at the elevations of potential avalanche starting zones. We only used, for our work, the 124 stations that measure snow properties in addition to weather properties. The measurements are recorded every 30 minutes (Fig. 1; see Appendix A for a list of parameters measured). These measurements are used to drive the simulation of snow stratigraphy with the 1-D snow cover model SNOWPACK (Bartelt and Lehning, Reference Bartelt and Lehning2002). We used the operational SNOWPACK setting (Snowpack DEFAULT version 20210925.e77eeeb; Wever and others, Reference Wever, Fierz, Mitterer, Hirashima and Lehning2014; Morin and others, Reference Morin2020) where liquid water percolation is modeled with a bucket-type approach (Lehning and others, Reference Lehning1999; Vionnet and others, Reference Vionnet2012; Lafaysse and others, Reference Lafaysse2017). Alternatively, water transport can be modeled with the Richards’ equation. However, this transport scheme is more demanding and not used in the operational setting for avalanche forecasting in Switzerland.

Figure 1. Location and elevation of the 124 automated weather stations (IMIS network) that measure snow properties.

Stratigraphic profiles are simulated every 3 hours at the location of each station, including simulations of the snowpack for four virtual slopes (aspects asp: N, E, S, W) with an incline of 38$^\circ$ (Morin and others, Reference Morin2020), resulting in 8 profiles per day and aspect.

If not stated otherwise, we always refer to this nowcast setup.

2.2.3 Variable definition

From the snowpack simulations, we extracted weather parameters and properties describing the state of the snowpack. As the temporal resolution of the avalanche observations was often restricted to the release day rather than an exact time, we derived 24 h averaged weather and snowpack features (00:00–24:00 local time, LT), as well as weather and snowpack properties corresponding to the time of the wettest profile of the day (Fig. 2). We defined the wettest profile as the profile with the highest LWC index (LWC_index; described below); in most cases, the profile at 15:00 LT was the wettest. If the LWC_index was zero for all time-steps, we used the data from the simulated profile at 15:00 LT. We chose this setup to account for daily weather and snowpack properties and to have a snapshot of the situation, when the snowpack was wettest (similar to Wever and others, Reference Wever, Vera Valero and Techel2018). In addition, we computed the differences to the values describing snowpack properties 1, 2 and 3 days before. The complete list of input features used for model development is reported in Appendix A. LWC-related features used were mostly based on Wever and others (Reference Wever, Vera Valero and Techel2018). We extracted in total 146 of features (Tables 5 and 6). Those included 36 features related to LWC, 9 of them were extracted from the wettest profile of the day and the remaining 27 describe changes in LWC over the last 1, 2 and 3 days. The other 110 features were meteorological and snowpack variables describing either daily averages (flagged with the subscript daily) or actual values at the time of the wettest profile (see Fig. 2).

Figure 2. Temporal characteristics of the input features: 24 h averaged weather features (gray area), weather and snowpack features at the time of the wettest profile (for example at 15:00 LT). Red stars indicate observed wet-snow avalanches.

From the simulated profiles, the LWC of each layer can be derived. It allows to compute the so-called LWC index, LWC_index, a feature that has been suggested to predict wet-snow avalanches (Mitterer and others, Reference Mitterer, Techel, Fierz and Schweizer2013). The index is the average of the LWC, in percent, of the entire snowpack weighted by layer height and divided by 3 %, i.e.

(1)$$\hbox{LWC}_{\hbox{index}} = {1\over 3}{\sum_{\hbox{layer}}h\left(\hbox{layer}\right)\hbox{LWC}\left(\hbox{layer}\right)\over \hbox{HS}}$$

where the sum carried out over all snowpack layers, with ${\rm HS} = \sum _{{\rm layer}}h( {\rm layer})$ the snowpack height, h the layer thickness and LWC the LWC in percent volume. Above a LWC of 3 %, it is assumed that excess water drains to the layers below since it is no longer held between grains by capillary forces, i.e. the threshold of 3 % corresponds to the transition from the pendular to the funicular regime. Therefore, ${\rm LWC}_{\text{index}} \ge 1$ may indicate a particular loss of strength (Mitterer and others, Reference Mitterer, Techel, Fierz and Schweizer2013), or a snowpack wet throughout.

2.3 Linking avalanche observations to SNOWPACK simulations

The occurrence of wet-snow avalanches depends on the snowpack conditions, which may differ strongly between slope aspect and elevation (Mitterer and Schweizer, Reference Mitterer and Schweizer2013; Wever and others, Reference Wever, Vera Valero and Techel2018). Therefore, we linked an observed avalanche in the vicinity of a station st and SNOWPACK simulation for aspect asp if:

1. 1. The avalanche was classified as a wet-snow avalanche, a size was estimated and the aspect and elevation were reported.

2. 2. The avalanche released within an elevation band of ±250 m of the elevation of the station, H _st. In case of avalanche summary reports indicating an elevation range (av _elev.min-max), avalanches were assigned to this station relative to the overlapping proportion between this elevation band (H _st ± 250 m) and the elevation band indicated in the avalanche summary report. For example, if H _st was at 2000 m (± 250 m: 1750 to 2250 m), and av _elev.min-max ranged from 2000 to 2500 m, the overlapping elevation band is 2000 to 2250 m, which is 50 % of the indicated av _elev.min-max. In this case, each of the avalanches in the summary report was weighted with 0.5.

3. 3. The avalanche was within a distance of at most 8.9 km, 17.8 km and 39.9 km corresponding to circle-shaped areas around the location of the station of 250 km², 1000 km² and 5000 km².

4. 4. The avalanche released in the same slope aspect as the simulation was performed for. For example, a wet-snow avalanche observed on a north-facing slope was only linked to the simulated stratigraphy on the virtual slope of northern aspect. If an avalanche was observed on a northeast-facing slope, this avalanche was assigned with half its weight w to both northern and eastern aspect slope simulations.

Figure 3 illustrates the matching requirements 2 to 4.

Figure 3. Schematic representation of how observed wet-snow avalanches on south facing slopes (red squares) were linked to input features derived from SNOWPACK simulations on a virtual south-facing slope at the location an AWS (station; red star). The pyramids represent different elevation bands and aspects for each observed avalanche and the virtual slopes at the station. The elevation bands show the ±250 m tolerance band. The circles represent the borders of areas surrounding the station. Green dots represent avalanches not assigned to the station because they did not satisfy the elevation and/or aspect criteria.

The avalanche summary reports prior to October 2019 often did not allow to extract the number and size of avalanches by aspect, as dry-snow and wet-snow avalanches could be reported in the same form. Therefore, we could rarely use these data in our analysis.

To quantify avalanche activity in a region, the avalanche activity index (AAI) was used (Schweizer and others, Reference Schweizer, Jamieson and Schneebeli2003). For each station st, we determined the AAI for several spatial scales a (area) for the four slope aspects asp as:

(2)$$AAI_{st}^{asp}\left(a\right) = \sum_{s = 1}^{5} w( s) N_s^{asp}\left(a\right),\; $$

with avalanche size s ∈ {1,  2,  3,  4,  5}, weight w(s) as defined in Table 1, the number of avalanches of size s, $N_s^{asp}( a)$ reported in the area a ∈ {250,  1000,  5000 km²} and slope aspect asp ∈ {north,  east,  south,  west}. The computation of AAI was constrained by restricting it to a specific elevation band around a station: ±250 m.

3.1 Definition of the target variable: avalanche day and non-avalanche day

We defined a binary target variable to discriminate days with wet-snow avalanche activity (avalanche days, AvD) from days without any activity (non-avalanche days, nAvD). Even though avalanche observations are reported from throughout the Swiss Alps, there are regions where rarely any observations are reported. And even in the areas generally covered by observers, observations are not always possible. Therefore, the absence of an avalanche in the data base does not necessarily indicate that no avalanche occurred. Beside the challenge of detecting days when no wet-snow avalanches occurred, we are aware of two main sources of error or uncertainty relating to avalanche observations:

1. 1. The classification of the avalanche as either a wet-snow or a dry-snow avalanche. According to the observational guidelines (SLF, 2020), the presence of liquid water should be estimated based on the conditions in the release area. However, since the avalanche release area is often not accessible, the LWC of the snow in the release area is often an estimate from the valley floor.

2. 2. The release date. The uncertainty in the reported release date depends on the frequency of observations in an area, on weather conditions (e.g. visibility), but also on the overall quality of the observations provided by a specific observer.

To address these points, which are highly relevant for this study, we opted for an approach allowing us to check the plausibility of the release date and avalanche type classification by considering observations independently made on the same day by observers in neighboring areas. More specifically, we define an avalanche day or non-avalanche day, our target variable $Y_{st}^{asp}$, associated with data from a station st for a given slope aspect asp as:

(3)$$Y_{st}^{asp} \!= \!\! \left\{\matrix{ \! 0 \hfill & \! \hbox{if}\; AAI( Swiss Alps) = 0 \hbox{ for all elevations and aspects}\hfill \cr \! 1 \hfill & \! \hbox{if}\; AAI_{st}^{asp}( 250\, {\rm km}^2) \! \geq \! 0.1 \ \hbox{\amp}\ \hfill AAI_{st}^{asp}( 1000\, {\rm km}^2) \geq 0.4 \;\hbox{\amp} \ \hfill & \hfill \cr \hfill & \! AAI_{st}^{asp}( 5000\, {\rm km}^2) \geq 2 \ \hbox{\amp} \hbox{ gap check}\hfill \cr \! \hbox{NaN} \hfill & \! \hbox{otherwise,} \hfill }\right.$$

with NaN ‘not a number’. The gap check requirement was $AAI_{st}^{asp}( 5000 \, {\rm km}^2) > AAI_{st}^{asp}( 1000\, {\rm km}^2) \;\hbox{\amp}\; AAI_{st}^{asp} ( 1000 \, {\rm km}^2) \gt$ $AAI_{st}^{asp}( 250\, {\rm km}^2)$, ensuring that avalanche activity increased with increasing area and was not a local phenomenon.

The choice for selecting an AAI threshold of 0.1 was guided by the fact that an $AAI_{st}^{asp} ( 250\, {\rm km}^2) \geq 0.1$ either means that there was at least one potentially life-threatening size 2 event (Eqn. (2); Table 1) or 10 avalanches of size 1. Moreover, an area a of 250 km² corresponds approximately to the mean size of the warning regions, the smallest spatial units used for communication of avalanche danger in the Swiss avalanche bulletin (e.g. Techel and others, Reference Techel2018), and 5000 km² is approximately the size of a snow-climate region (e.g. Techel, Reference Techel2020, p. 50). Within a warning region there is often at most one observer providing regular avalanche observations; within a snow-climatological region there are typically between 5 and 20 observers. With the definition given above, we ensure that local avalanche activity with $AAI_{st}^{asp} ( 250\, {\rm km}^2) \geq 0.1$ is representative of avalanche activity in a larger area. We chose this approach to ensure that the target variable was not a local property, but included several independent observations, thus reducing noise due to erroneous dating or wetness classification. We considered this step as important, as $AAI_{st}^{asp}$ only considers observations within a specific slope aspect and elevation band. We used a similarly strict criterion to define a non-avalanche day (nAvD: $Y_{st}^{asp} = 0$) as the uncertainty related to observations of no avalanches is much higher. Thus, a day is considered only as nAvD if none of the observers in the Swiss Alps reported a wet-snow avalanche. Furthermore, we considered only nAvD for station-aspect combinations that had at least one AvD in the same season.

With this definition for the target variable, we wanted to separate days with widespread avalanche activity (AvD) of a certain magnitude from days with absolutely no avalanches (nAvD), and excluded days with rather local avalanche activity (NaN values). With regard to the development of the model, our strict definition led us to train and test the classifier on rather extreme cases, which are however comparatively reliable in terms of the quality of the class label.

3.2 Labeling and splitting the dataset

Input data were labeled as AvD and nAvD according to Eqn (3). We split the data into three uncorrelated datasets according to the winter the avalanches occurred (Fig. 4): Dataset1 contained AvD (Y = 1) and nAvD (Y = 0) labeled data for the winter 2019–2020 and historical AvD events between 2001 and 2014. Dataset2 contained data for winter 2020–2021. Operational tests were carried out during the winter 2021–2022 (dataset3). The dataset characteristics are summarized in Table 2. The stringent manner, we defined the target variable (Eqn. (3)), considerably reduced the number of AvD events. For instance, dataset2 and dataset3 originally contained 2595 and 1696 records with AAI ≥ 0.1, respectively. However, only 223 (dataset2) and 97 (dataset3) were labeled as AvD according to our target variable definition (Table 2).

Figure 4. Number of days labeled as AvD for dataset1 (blue), dataset2 (orange) and the dataset3 (green).

Table 2. Number of cases in the datasets used for training and testing

In parentheses, the number of individual calendar days are given; for example, during winter 2021–2022 there were 97 AvD events on 12 different calendar days.

In the following, we name ’test set’ the complementary of the training set, for example, if we used dataset1 as the training set, the test was dataset2.

3.3 Random forest model development

We trained a RF model, as it is one of the most common and robust machine learning models to classify tabular data (Breiman, Reference Breiman2001). We used the scikit-learn library in Python (Pedregosa and others, Reference Pedregosa2011).

The datasets are heavily imbalanced and contain only relatively few AvD events (Table 2). As the datasets are also rather small, we did not use a standard approach of splitting the data in a training, a validation and a test set. Instead, we only used a training and a test set.

We relied on the out-of-bag (OOB) score to assess model performance (Probst and others, Reference Probst, Wright and Boulesteix2019). Each tree of a RF model was trained separately using a training subset randomly sampled with replacement. These subsets are named bootstrap samples. The left out samples, i.e. samples not contained in the boostrap samples, are named the OOB samples and were used to compute OOB scores, such as the f1-score (see Appendix C for metrics definitions). The OOB samples size was around 1/3 of the size of the training set (Breiman, Reference Breiman2001).

The RF hyperparameters were tuned to maximize the out-of-bag f1-score of class AvD of the training set using a grid search (Table 7).

In order to reduce the complexity of the model and remove non-relevant input variables, we applied recursive feature elimination (RFE; Guyon and others, Reference Guyon, Weston, Barnhill and Vapnik2002). First, we fitted a model keeping all the original 146 features and computed impurity-based features importance ranking (for implementation details see Breiman, Reference Breiman2017). Then we discarded the less important features, re-fitted the model and repeated the process until only 4 features remained. At each step, we computed the OOB f1-score of AvD, which allowed us to select the best model as a function of the number of features. The best performance was obtained using 40 features (Fig. 11). We used this set of features for the remaining of the work.

The final prediction of the RF model is obtained by taking the majority vote of the predictions made by each individual decision tree, and the probability of a particular class, here AvD or nAvD, is determined by the percentage of trees that have voted for that class. In the following, we used also the AvD probability to assess model performance.