Data

In order to establish a link between regional avalanche danger and modelled snow-cover variables, a test region was chosen where both modelled snow-cover data and an estimate of the regional avalanche danger were available. We selected the region of Davos in the eastern Swiss Alps (225 km²), which covers typical avalanche-release zones with an elevation range of 1600–2800 m a.s.l. The study plot in this test region indicates that the automated weather station Weissfluhjoch, at 2540 m a.s.l., shows an average maximum snow depth of ∼2.2 m for the nine winters (1999–2007; 1229 days) when the required data were available.

The forecast of the regional avalanche danger level (Fig. 1), which is issued every day at 0800 h, was used as a proxy target parameter since an accurate measure of the avalanche danger is not available (Reference MeisterMeister, 1995). The frequency of the danger levels in the chosen time period (black columns) can be seen in Figure 1 (1: ‘Low’, 2: ‘Moderate’, 3: ‘Considerable’, 4: ‘High’ and 5: ‘Very high’; ‘Very high’ did not occur in this time period). The danger levels ‘Moderate’ and ‘Considerable’ were most frequent. The frequency in the chosen time period was substantially different from that in the periods used in the studies of Reference Brabec and MeisterBrabec and Meister (2001) (dark grey) and Reference Schweizer and FöhnSchweizer and Föhn (1996) (light grey). For the danger levels ‘Low’ to ‘Considerable’ the avalanche danger was characterized by very high persistence. The probability that the avalanche danger tomorrow will be the same as today was ∼80%.

Fig. 1. Relative frequency of the avalanche danger levels in the region of Davos. Black columns show the forecasted levels during the time period considered in this study. Dark grey columns show the forecasted levels for the study of Reference Brabec and MeisterBrabec and Meister (2001) and light grey show the verified levels used in the study of Reference Schweizer and FöhnSchweizer and Föhn (1996).

As input variables, modelled snow-cover data were generated for the location of Weissfluhjoch. The relevant processes influencing the regional avalanche danger (e.g. new snow, wind or weak layer formation) are assumed to be represented in the chosen study plot. The snow-cover model SNOWPACK was used to model settling and layering of the snow cover, as well as its energy and mass balance (Reference Bartelt and LehningBartelt and Lehning, 2002; Reference Lehning, Bartelt, Brown, Fierz and SatyawaliLehning and others, 2002a, Reference Lehning, Bartelt, Brown and Fierzb). This model requires meteorological data as input. We focused on dry-snow situations since dry-snow avalanches are the main threat over most of the winter. For forecasting wet-snow avalanches models need to be trained separately. For each winter a date was determined for the beginning of wet-snow conditions (i.e. 2 weeks before the snow cover modelled for the Weissfluhjoch study plot became isothermal), which is usually around the beginning of April. We furthermore restricted the dataset to days with snow depth >75 cm, because we think that the stability part of SNOWPACK produces more reliable results with this restriction. This restriction excluded only a limited number of days because the avalanche danger forecast for our region started at about the time when that snow depth was reached – for most winters by the end of November.

A stability index (Reference Schweizer, Bellaire, Fierz, Lehning and PielmeierSchweizer and others, 2006) defined the potential weak layer interface in the modelled snow cover. Motivated by a study which evaluated stability from observed snow stratigraphy (Reference Schweizer and JamiesonSchweizer and Jamieson, 2003), characteristics of the weak and adjacent layer and of the slab were considered. These model variables were completed withmeasured and calculated meteorological and snow-surface variables (e.g. wind velocity, outgoing longwave radiation or surface albedo).

All statistical methods described here were tested on a dataset in order to assess their classification capability. This test dataset (Fig. 1) was used for the study of Reference Schweizer and FöhnSchweizer and Föhn (1996) and contained a verification instead of a forecast of the regional avalanche danger level and observed, rather than modelled, snow-cover data as input variables. This dataset was more reliable than that covering the years 1999–2007, and therefore adequate to test the statistical methods for this specific problem. The test dataset covered ten winters between 1985 and 1994. Unfortunately no modelled SNOWPACK data can be produced for that time period, because no automated weather station existed before the 1990s. This meant that modelled and measured snow-cover data and their explanatory power could not be directly compared.

Performance measures

The quality of a method was assessed by a cross-validated hit rate (HR) for all danger levels and by the true skill score (TSS) for each of the four danger levels separately (Reference Doswell, Davies and KellerDoswell and others, 1990; Reference WilksWilks, 1995).

Since both input and target variables were autocorrelated (e.g. a weak layer might have an influence over a long time period), random cross-validation turned out not to be a useful method, giving unrealistically high hit rates. Therefore, each winter was forecasted using a model and a variable selection based on the remaining eight winters, i.e. both model parameters and input variables could change in this cross-validation scheme.

It seems useful to draw special attention to days on which the danger level changed: first, these days might be the most important to predict correctly; second, because of the high persistence of the target parameter, an obviously useless forecast which always predicts the level of the previous day, would show good hit rates for all days. Introducing HR and TSS for only those days on which the danger level changed helped to identify methods which tended to a persistent forecast. Similarly, the HR and TSS were also considered for the days on which the modelled danger level increased or decreased.

Variable selection

Since the snow-cover model delivers a large number of variables at high temporal resolution, variable selection is useful, firstly for data reduction, which increases the speed of the final statistical forecasting model, secondly to achieve an overall improved performance and thirdly to understand which variables are important (Reference Guyon, Elisseef, Guyon, Nikravesh, Gunn and ZadehGuyon and Elisseef, 2006).

For many variables, it may make sense to also consider their sum, mean, extreme values or rate, for different time intervals or time lags. For example, the 24 hour change of the air temperature may be more correlated with the regional avalanche danger than the midday value, or the previous day’s snowfall more correlated than that of the current day. This leads to a rapid increase in the number of possible variables. A simple univariate variable selection, a Fisher’s discriminant analysis (e.g. Reference BishopBishop, 2006), was performed to determine for each variable the two most important derived variables. After this step, a large number (300) of variables still remained. Subsequently a rating of this variable set was carried out with Fisher’s discriminant analysis. An alternative ranking was obtained with univariate classification trees (Reference Breiman, Friedman, Stone and OlshenBreiman and others, 1998). To avoid overfitting, the best tree was determined by a pruning algorithm based on cross-validation, and the number of terminal nodes was limited to ten. The cross-validated HR of these univariate pruned trees delivered a ranking between the variables. Only variables that were not pairwise linearly correlated (r ² < 0.6) were considered. Since correlated variables are not redundant per se, the redundancy of excluded variables was visually examined with scattergraphs (Reference Guyon, Elisseef, Guyon, Nikravesh, Gunn and ZadehGuyon and Elisseef, 2006). For support vector machines (SVM) a variable selection algorithm was used, based on the Fisher’s discriminant analysis in combination with SVM (Reference Chen, Lin, Guyon, Nikravesh, Gunn and ZadehChen and Lin, 2006).

Additionally, for the nearest-neighbour method (KNN) a combination with a genetic algorithm (GA) was used to rank variable relevancy (Reference Li, Weinberg, Darden and PedersenLi and others, 2001). Several variable subsets (each of ten variables) were tested with the KNN method. ‘Good variable subsets’ (see below for the selection criteria) were stored in a final pool. Once the final pool was filled to a certain threshold, the cardinality of each single variable was interpreted as a relevancy ranking: the more often a single variable was selected as a member of a good variable subset, the more relevant the single variable was assumed to be. Since testing all possible ten-member variable subsets of 300 different variables would cost too much computing time, a GA was used as a search tool to create ‘good parameter subsets’: the GA maximized a fitness function which depended on the selected variable subsets. The fitness function was a combination of the three different hit rates described above (i.e. for all days, for the days on which the avalanche danger level changed and for the days on which the modelled danger level changed). In more detail, the GA was initiated with 100 randomly selected ten-member variable subsets. The GA then determined a maximum by continuous mutation and stored it in the final pool. This procedure was reiterated until the final pool contained 100 variable subsets, which was enough to achieve reproducible results.

For all methods variables were scaled linearly to [−1, C. − C1], except for categorical variables which were translated to Booleans for each category. As input for the statistical methods, subsets of the 5, 10 and 30 most important variables of each algorithm were tested.

Statistical methods

In order to find an optimal link between input variables and predicted regional avalanche danger, we used the following statistical methods:

• classification trees (TREE),

• artificial neural networks (ANN),

• nearest-neighbour methods (KNN),

• support vector machines (SVM),

• hidden Markov models (HMM).

In the following each method is described in more detail.

A simple classification tree with only one measured variable, the 3 day sum of the new snow (HN3d_meas), performed very well and was therefore used as a benchmark for more complex models with more (especially modelled) input variables. For this tree and more complex trees using more input variables, generalization was achieved with a pruning algorithm using ten-fold cross-validation on the training sets.

Artificial neural networks, both feed-forward and recurrent (e.g. Reference ElmanElman, 1990; Reference BishopBishop, 2006), with different set-ups (hidden layer size, training algorithms, early stopping or Bayesian regularization) were trained for each winter. Since results depend on the initial weights, the networks were initialized five times and the mean of the results was considered. The results of ANNs discussed in the next section were gained with a recurrent network which used adaptive learning, 100 hidden neurons and 100 passes through the sequence.

The nearest-neighbour approach as used by Reference Brabec and MeisterBrabec and Meister (2001) was applied for a direct comparison to previous work. In addition, this method was modified in two ways. (1) The avalanche danger level of the previous day was used as an additional input to predict the current day. With this information the training dataset was reduced: only those days that showed the same avalanche danger level on the previous day were considered as possible nearest neighbours for the current day. This implicated an even more unbalanced classification problem due to the persistence of the target parameter. (2) The classification was modified. A decision boundary was used to determine whether the danger level should change to a new value on a particular day. Using ten nearest neighbours, three of these days must show a higher (or lower) danger level. This decision boundary was obtained by optimizing the TSS (Reference Heierli, Purves, Felber and KowalskiHeierli and others, 2004). Subsequently the danger level was determined with a majority vote between the neighbours showing higher (or lower) danger levels. For breaking a tie the nearest neighbour among them was used. Weights in KNN methods allow certain variables to have more influence while calculating the distance between nearest neighbours. Optimal weights were determined with a GA (Reference Purves, Morrison, Moss and WrightPurves and others, 2003).

SVM (Reference Schölkopf and SmolaSchölkopf and Smola, 2001) were applied to the problem using the 2001 LIBSVM software package of C.-C. Chang and C.-J. Lin (software available at www.csie.ntu.edu.tw/cjlin/libsvm). Gaussian radial basis functions with radius γ were used as kernel functions. Cross-validation on the training sets was performed to obtain γ and the penalty variable, C (Reference Chen, Lin, Guyon, Nikravesh, Gunn and ZadehChen and Lin, 2006).

In addition, hidden Markov models were used, since they are able to predict time series (Reference RabinerRabiner, 1989; Reference BishopBishop, 2006). To obtain a discrete input for the HMM, the continuous input vectors were mapped into a discrete codebook index with K-means vector quantization.

Variable selection

Since the variables were selected for each forecasted winter separately, we present in this section a summary of the results obtained for each winter. In Table 1 the most important variables are listed, selected by the three methods: Fisher’s discriminant analysis (Fisher); univariate classification trees (TREE); and the combination of genetic algorithm and nearest neighbour (GA/KNN) described above. In the last column the sign of correlation with the avalanche danger level is given to allow a plausibility check. Except for the variables strain rate of the weak layer, 3 hour rate of crust thickness and net longwave radiation, a physical interpretation can easily be given. Important variables were the new-snow depth (HN, HN3d) and the other new-snow variables with larger values at higher avalanche danger levels. Surface albedo, 3 hour rate of the slab thickness and relative humidity can also be related to new-snow situations, although they are not pairwise linearly correlated. The snow transport index (Reference Lehning and FierzLehning and Fierz, 2008) was also selected by all three algorithms with larger index values (more drifting and blowing snow) at higher levels. Higher maximum wind speeds in the last 24 hours were also correlated with higher avalanche danger levels. The deformation index (Reference Lehning, Fierz, Brown and JamiesonLehning and others, 2004) was also found to be important, despite the fact that a bug has recently been discovered in the description of the temperature dependence, although the index should still describe a relation between critical and actual stress in the bonds. As the definition suggests, small index values are correlated with high avalanche danger. The variable crust is positively correlated to the number of elapsed days without snowfall and therefore negatively correlated to avalanche danger. The modelled profile type ‘Four’ (Reference Schweizer and LütschgSchweizer and Lütschg, 2001) was related to the danger level ‘Considerable’. The profile type ‘Four’ is the most frequent in the dataset which describes a weak base. Although a weak base is not by itself conclusive (Reference Schweizer and WiesingerSchweizer and Wiesinger, 2001), it makes sense that a weak base is related to the danger level ‘Considerable’ as a sign of structural instability. Small bond size and low density of the weak layer might also be considered as signs of structural instability.

Table 1. Overview of the variable selection results

A result of a variable selection with the combination of GA/KNN methods can be seen in Figure 2, which shows every single variable on the x axis and their frequency in the final pool on the y axis. Variables with a high pick-frequency were interpreted as important. This method selected, amongst others, the 3 hour rate of the outgoing longwave radiation, which was highly correlated with the rate of the snow-surface temperature. Strong warming of the snow surface before midday was related to lower avalanche danger. It is worth noting that this method chose especially meteorological or snow-surface variables. This may be because the avalanche danger level of the previous day was introduced as an additional input. This already provides a certain level of stability to the system so that snowstratigraphy variables become less important. Accordingly, only weather or surface properties were important because they describe the change in danger level. To confirm this hypothesis, the GA/KNN method was performed without this additional information. However, the classification power of the KNN method without this additional information of the danger level on the previous day was too poor to reach a conclusion. Better methods were too computationally demanding for combination with the GA.

Fig. 2. Variable selection with the GA/KNN method for winter 1999/2000. Frequency (y axis) for each single variable (x axis) in the final pool. The 12 variables selected are marked; the two most important were (1) HN and (2) the 3 hour rate of outgoing longwave radiation.

Statistical methods

The a priori capability of the statistical methods was compared with the cross-validated HR of the test dataset (verified danger levels). Classification trees had the lowest performance, while SVM, ANN, HMM and the KNN method using the avalanche danger level of the previous day as an additional input achieved HRs of ∼60%. These results were comparable to the DAVOS4 model (63%) used by Reference Schweizer and FöhnSchweizer and Föhn (1996) and the Kohonen neural network (61%) used by Reference Schweizer, Föhn, Schweizer and HamzaSchweizer and others (1994).

Since the distribution of the forecasted avalanche danger level of the dataset covering the years 1999–2007 (Fig. 1) was substantially different to previous studies (Reference Schweizer and FöhnSchweizer and Föhn, 1996; Reference Brabec and MeisterBrabec and Meister, 2001), results were not directly comparable. Therefore the method used by Reference Brabec and MeisterBrabec and Meister (2001) with measured meteorological data as input (BRABEC) was applied to the time period of this study.

In the following, results will be given for a selection of models, which represent the range of model quality. No results for SVM are presented since they delivered no additional improvements. The selected models and their input variables are:

Table 2. Input variables used in the BRABEC model (Reference Brabec and MeisterBrabec and Meister, 2001)

Figure 3 shows the performance of the selected methods considering all days. All new methods reached a higher HR than the BRABEC method (55%). Even the simple TREE method achieved a remarkable HR of 65%. However, the TREE method produced only the most frequently used danger levels ‘Moderate’ and ‘Considerable’. Thus the TSS values of the other danger levels were 0%. For most of the nine winters the split value for 3 day sum of new snow (HN3d_meas) was 12 cm, which discriminated between the danger levels ‘Moderate’ and ‘Considerable’. The variable HN3d_meas was also amongst the input variables for the BRABEC method (Reference Brabec and MeisterBrabec and Meister, 2001). Using prior probabilities, classification trees are, in general, also able to handle unbalanced target parameters as exist in our case (Fig. 1). Whereas this improved the TSS value for the danger levels ‘Low’ and ‘High’, it produced insufficient overall hit rates. Also more variables did not improve the results.

Fig. 3. Comparison of the various model performances for all days. Abbreviations are explained in the text. The black column is the cross-validated HR of all danger levels, the grey columns the cross-validated TSS value for each of the four danger levels (lightest to darkest indicate ‘Low’ to ‘High’).

The other methods achieved better hit rates and better true skill scores for each single danger level. With an HR of 73%, the KNN method performed best.

Considering all days, the KNN_vg method (as KNN, but with measured, not modelled, input variables) achieved similar results to the KNN method. However, Figure 4 indicates the weakness of that method; it shows the same methods and performance measures as in Figure 3, but only for those days on which the target parameter changed. The low values in HR and TSS indicate that the KNN_vg method did not predict these important days satisfactorily. All other methods showed similar hit rates. For the KNN method it is interesting that the TSS values were greater for the higher avalanche danger levels. It seems to be easier for this model to predict an increase in avalanche danger than a decrease. The HMM and ANN methods best predicted a decrease in the avalanche danger, which can be seen by the more balanced TSS values between each danger level. Also the simple TREE method was reasonably able to predict the days on which the target parameter changed. As the variable HN3d_meas is used as input, it may not be surprising that the TREE method was able to predict an increase in danger level, but a decrease was also predicted correctly in ∼50% of all cases.

Fig. 4. The black column is the cross-validated HR of all danger levels, the grey columns the cross-validated TSS value for each of the four danger levels (lightest to darkest indicate‘Low’ to ‘High’).

Since the detection of days on which the avalanche danger decreased or increased is important information, models were trained for the binary criteria ‘Increase or not’ and ‘Decrease or not’. However, no improvement was achieved for these, at first glance, simpler classification problems compared to the methods which forecast the four avalanche danger levels. In the case of decreasing danger, one possible explanation is that it might be very difficult for human forecasters to decide consistently whether the avalanche danger level should decrease on a certain day and not on the day(s) before or after. Therefore, it might also be difficult to find reasons in the parameter space of the input variables presented to the statistical methods. In the case of increasing danger, the reasons for the change in danger level might differ between danger levels, especially between an increase from ‘Low’ to ‘Moderate’ and from ‘Considerable’ to ‘High’. This information is lost if the models are trained only on the binary criteria.

As we saw an improvement when adding the avalanche danger level of the previous day as additional information to the KNN method, this procedure was also implemented for the other methods. Due to the persistence of the target value, this implicated a highly unbalanced classification problem. Although most methods provided a formal method to include this additional information, (e.g. prior probabilities (Reference Breiman, Friedman, Stone and OlshenBreiman and others, 1998) or different penalties in the training error function), none of the methods showed an improvement.

A final test was whether the spatial variability of the input variables had additional classification power. Simulations of virtual slopes of the four main aspects and the flat field at two automated weather stations (Weissfluhjoch (2540 m a.s.l.) and Hahnengretji (2490 m a.s.l.)) in the test region were used as input for the statistical methods. The variability was implemented either (1) as the maximum difference between the ten simulations of each variable or (2) with additional categorical data describing the station and the aspect of the simulations. Neither option improved the performance of the statistical methods.

Figure 5 shows a comparison between the KNN method and the forecasted avalanche danger for a representative winter.

Fig. 5. The avalanche danger forecasted and modelled with the KNN method for winter 1999/2000.