Introduction
In the Columbia Mountains, western Canada, most fatal avalanches are skier-triggered dry-slab avalanches. Many natural avalanches are a direct result of snowfall events, whereas many skier-triggered avalanches occur during fair weather. Strength and stability indices of buried persistent weak snowpack layers are highly correlated with skier-triggered avalanches, but forecasting models have made little use of stability indices, partly because incorporation requires daily strength values and partly because spatial variability of snowpack properties may compromise extrapolation from the study plot to the regional scale (Reference Hägeli and McClungHägeli and McClung, 2001; Reference Landry, Birkeland, Hansen, Borkowski, Brown and AspinallLandry and others, 2002). In current computer-assisted forecasting models, mainly datasets for natural avalanches or a combination of natural and skier-triggered avalanches are used (e.g. Reference McClung and TweedyMcClung and Tweedy, 1994; Reference Purves, Morrison, Moss and WrightPurves and others, 2003). A recent empirical model based on shear frame measurements for the Columbia Mountains forecasts the strength of surface hoar layers up to 8 days ahead using snow profile data (Reference Chalmers and JamiesonChalmers and Jamieson, 2003). This paper uses a similar empirical model to forecast the strength of layers of faceted crystals based on load and loading rate (Reference Zeidler and JamiesonZeidler and Jamieson, 2002). A skier stability index, Sk38 (Reference Jamieson and JohnstonJamieson andJohnston, 1998), is calculated daily using these two weak-layer models. Sk38, slab thickness (H) and slab load (Load) on the weak layer are incorporated into a nearest-neighbour model to verify that stability indices and slab properties improve the number of successfully forecasted days on which skier-triggered avalanches occurred or were likely. The aim of this paper is (1) to show the predictive potential of the snow-pack properties Sk38, H and Load using rank correlation and tree analysis and (2) to incorporate Sk38, H and Load into a nearest-neighbour model and verify that these stability measures improve the performance of a nearest-neighbour model for skier-triggered avalanches. Our focus is to better forecast skier-triggered dry-slab avalanches on persistent weak layers and does not take into account loose-snow avalanches, wet avalanches or dry-slab avalanches on weak layers of precipitation or decomposed and fragmented precipitation particles (non-persistent weak layers).
Literature Review
Nearest-neighbour models
For several years nearest-neighbour models have been used operationally and they have proven to be a useful forecasting tool (e.g. Reference BuserBuser, 1989; Reference Purves, Morrison, Moss and WrightPurves and others, 2003). These models are applied and verified mostly for natural avalanche activity and not as much for skier-triggered avalanches. Nearest-neighbour models calculate a numerical distance for each one of the input variables between the value from the current day and values from the historical database, and the days with the closest distances (nearest neighbours) are determined. A list of avalanches on the generally 10 or 30 nearest neighbours is displayed and interpreted by the forecaster. Input variables range from weather and snowpack data to elaborated variables such as settlement of storm snow (Reference McClung and SchaererMcClung and Schaerer, 1993, p. 157). Weighting of input parameters allows consideration of the physical understanding of avalanche release.
Reference Brabec and MeisterBrabec andMeister (2001) consider local avalanche forecasting as the task of predicting the probability of avalanche release for a small area (e.g. in a ski area of about 100 km2). Models like NXD2000 (Reference Brabec and MeisterBrabec and Meister, 2001) and Cornice (Reference Purves, Morrison, Moss and WrightPurves and others, 2003) are normally used as a tool in local avalanche forecasting. In this study we apply Cornice on a regional scale, which requires the extrapolation of studyplot measurements to the regional scale. This extrapolation is questioned by Reference Hägeli and McClungHägeli and McClung (2001) based on scaling theory, and requires less spatial variability than reported by Reference Landry, Birkeland, Hansen, Borkowski, Brown and AspinallLandry and others (2002) for the Rocky Mountains of Montana, U. S.A.
Forecasting skier-triggered avalanches
About 90% of the fatal skier-triggered avalanches in Canada occur on persistent weak layers (Reference Jamieson and JohnstonJamieson and Johnston, 1992). This makes it important to look at the properties of these weak layers and the overlying slabs when forecasting skier-triggered avalanches. Reference JamiesonJamieson (1995, p.215–221) showed that a snowpack stability index was the most significant indicator of skier-triggered avalanches on persistent weak layers. As pointed out by Reference Obled and GoodObled and Good (1980), one avalanche on a clear day in a tourist area is potentially more dangerous than several during stormy weather. Avalanches on a clear day are harder to forecast and the number of people exposed is usually greater as well. Jones and Reference Jamieson and JohnstonJamieson (2001) assessed the importance of meteorological variables and previous avalanche activity with respect to skier-triggered dry-slab avalanches in the Columbia Mountains and explored their interaction. They found previous avalanche activity, new snowfall, storm snow, precipitation, height of the snowpack, days since 1 December, relative humidity and barometric pressure as the variables highest ranked in terms of their correlations with the response variable, the largest size of triggered avalanches. In this study, we build on their analysis using a similar dataset and try to improve the forecast of skier-triggered avalanches by including a stability index based on weak-layer strength and slab properties. Instead of forecasting the destructive potential (Canadian avalanche size), this paper focuses on stability as indicated by occurrence or absence of skier-triggered avalanches and the end-of-the-day stability ratings from ski guides. This is partly because the response variable in commercially available nearest-neighbour models is avalanche day (true or false) and partly because in this study all avalanches on persistent weak layers are considered as potentially harmful to a skier. Also, Reference Obled and GoodObled and Good (1980) assume that in a small region, if one avalanche occurs it may be assumed that all similar slopes, or even the whole region, are suspect. If so, it can be argued that an observed small avalanche might indicate that larger avalanches are likely on a bigger slope. Reference Obled and GoodObled and Good (1980) also point out that many small avalanches are as dangerous as one bigger avalanche.
It is more important to forecast unstable days correctly than stable days since a day on which an avalanche occurred but was not correctly forecasted has greater consequences than a day without an avalanche where an avalanche was forecasted (Reference Jamieson and JohnstonJamieson and Johnston, 1993; Reference Blattenberger and FowlesBlattenberger and Fowles, 1995; Reference JamiesonJamieson, 1995, p. 140–141).
Data and Nearest-Neighbour Model
The study area lies in the Columbia Mountains near Blue River, British Columbia. Meteorological variables are available from a remote weather station at Mount St Anne at an altitude of 1900 m. Snowpack data and strength values were recorded once or twice per week in a level study plot within 300 m of the meteorological station. Stability indices calculations were possible for seven winters from 1995/96 to 2001/02. Snowpack parameters such as slab load (Load) and slab thickness (H) as well as the stability index Sk38 were extrapolated for days without manual snowpack observations to obtain daily values for a nearest-neighbour model. Avalanche activity data and stability evaluations were provided by a heli-skiing operation near Blue River. From our entire dataset we extracted time periods in which a persistent weak layer was present in study-plot snow profiles and tested with the shear frame. Typically testing starts when a persistent weak layer has been buried for 2– 5 days, and continues until the weak layer has been buried for approximately 6 weeks. The testing of the persistent weak layers did not always begin when avalanche activity on the layer was first reported. On days between the burial date and the first measurement we calculated Load by subtracting the average loading rate per day from the first measured value and estimated H from the daily measurements of the height of the snowpack and the height of new snow. Days were excluded when the calculated ski penetration used for Sk38 exceeded H (Reference Jamieson and JohnstonJamieson and Johnston, 1998).
In the case of more than one persistent layer in the study plot, we chose the weakest one since the shear strength of a weak layer correlates strongly with its stability (Reference Föhn and CamponovoFöhn and Camponovo, 1997). The dataset contains 515 days during which the shear strength of a total of 17 different weak layers was measured. The dataset is composed of four time series of faceted crystals as a weak layer (139 days) and 13 time series of surface-hoar weak layers (376 days).
Our study area is in total approximately 5000 km2, whereby most skiing over the winter is within 30 km of the study plot and the weather station (Reference Jones and JamiesonJones and Jamieson, 2001). In our case we apply Cornice, a nearest-neighbour model developed in Scotland (Reference Purves, Morrison, Moss and WrightPurves and others, 2003), as a tool for regional avalanche forecasting. As mentioned, this requires the extrapolation of study-plot measurements to the regional scale. In the Columbia Mountains this is a reasonable approach since many significant persistent weak layers form throughout most of the area (Reference Hägeli and McClungHägeli and McClung, 2002), and empirical studies showed good correlations between study-plot measurements and avalanche activity in the whole study area (Reference JamiesonJamieson, 1995, p. 211–221).
The main features of Cornice that facilitate use and testing of the results are the automatic weighting function of input variables using a genetic algorithm (Reference Purves, Morrison, Moss and WrightPurves and other, 2003) and a batch forecast. The batch forecast calculates the nearest neighbours for each day in the database and gives the number of correctly forecasted dayswith the response variable equal to zero or one, as well as the combined percentage.
Methods
Response variable
A daily skier instability index (DSI) was developed. DSI was assigned a value of one for each day in which one or more skier-triggered dry-slab avalanches on a persistent weak layer were reported, or the ski guides operating in the area rated the stability at the end of the day as fair–poor, poor or very poor (Canadian Avalanche Association, 2002, p. 74), and DSI = 0 for the other days. This index allowed us to include days on which skier-triggered avalanches were likely, but due to limited conditions for flying (and hence helicopter skiing) or the terrain selection of the guides, no dry-slab avalanches on persistent weak layers were skier-triggered. We defined 128 unstable days (DSI = 1). Skier-triggered avalanches in non-persistent weak layers or low stability ratings due to storm snow instabilities were defined as stable days (DSI = 0).
Predictor variables
The variables considered in this study are listed in Table 1.
Daily predictor variables
Three snowpack variables were incorporated into Cornice. These are Sk38, Load and H.
Sk38: This stability index is the ratio of shear strength of a weak layer measured with a shear frame to the shear stress due to the overlying slab and a hypothetical skier, taking into account the ski penetration estimated from slab thickness and density, and the slope angle (Reference JamiesonJamieson, 1995, p. 141–156; Reference Jamieson and JohnstonJamieson and Johnston, 1998). This index is not calculated when the calculated ski penetration is greater than the overlying slab. As shown by Reference Chalmers and JamiesonChalmers and Jamieson (2003), Sk38 has the potential to predict the regional stability of buried surface-hoar layers in the Columbia Mountains.
Load: Load is the slab density integrated over slab thickness. The strength of weak layers increases with slab load (Reference Jamieson and JohnstonJamieson and others, 2001). Lags of a few days are likely for persistent weak layers (Reference ChalmersChalmers, 2001, p.79–84).
H: Slab thickness is an important variable to be considered when looking at skier triggering. The shear stress induced by skiers on a weak layer decreases with increasing depth of the weak layer. Consequently, deeper buried weak layers are harder to trigger than shallower weak layers.
Slab thickness (Equation (1)) and Load (Equation (2)) were predicted on days with no snowpack observations:
where Hi is the slab thickness on the day to be forecasted, Hi - 1 is the slab thickness on the previous day, and 0.95 an average daily settlement factor
where Loadi is the slab load on the day to be forecasted, Loadi - 1 is the load on the previous day, and 0.062 kPa d-1 an average loading rate for tree-line elevation in the area around Blue River (Reference Zeidler and JamiesonZeidler andJamieson, 2002).
We use average values of settlement and loading rate, which are available in the morning when route selection decisions are being made, rather than measured values, which would not be available until hours later when guides and technicians might visit the study plot.
Calculating Sk38 required shear strength values of weak snowpack layers.On dayswith snowpack observations these were measured using a shear frame (Reference Jamieson and JohnstonJamieson andJohnston, 2001). On days without manual snowpack observations, daily shear strength values were calculated by applying two empirical models: Reference Chalmers and JamiesonChalmers and Jamieson’s (2003) model to forecast the strength of surface hoar layers, and Reference Zeidler and JamiesonZeidler andJamieson’s (2002) model to forecast the strength of faceted layers. Reference ChalmersChalmers (2001) and Reference Chalmers and JamiesonChalmers andJamieson (2003) use load, slab thickness, height of snowpack, thickness of weak layer, temperature of weak layer, temperature gradient across weak layer and the minimum grain-size of the weak layer for their predictions. The forecast for the strength of faceted layers is based on slab load and loading rate. With values of the shear strength of the weakest persistent snowpack layer, overlying load and slab thickness, daily values of Sk38 can be calculated. This allowed us to compare the performance of Cornice for predicting unstable days on a daily basis using, and not using, Sk38, H and Load.
Statistical methods
Statistical methods are applied in this analysis to assess the association of snowpack properties, especially Sk38, Load and H, with DSI.
We apply two methods that do not assume the distribution of the variables since the response variable DSI has only two values (0 or 1) (Reference Jarrett and KraftJarrett and Kraft, 1989, p. 600):
Spearman rank correlations allow us to determine the degree of association of each predictor variable with the response variable (Reference Johnson and BhattacharyyaJohnson and Bhattacharyya, 1996, p. 632).
Classification trees are used to determine the importance of predictor variables associated with DSI when meteorological and snowpack variables are used in combination and to understand the interactions. Classification trees were constructed using all variables in Table 1 and using all these variables except Sk38, Load and H. Due to missing values in the dataset, the number of days used in the tree analysis is 411.
Classification trees consist of nodes that are each split into two subsets using a critical value of one variable (e.g. as seen in Figure 1 where the first node, including all data points, is split into two subsets using a critical value of Sk38 ≤1.33). The resulting child nodes are split until they become terminal nodes, either when only cases belonging to the predicted class are present in the node or when a node contains a minimum number of cases as defined by a stopping rule. In our analysis, stopping occurs when there are five or fewer cases in a node. Applying stopping rules in classification-tree analysis prevents the tree from being grown to a perfect fit for the dataset whereby low-order splits would have little predictive value for other datasets based on similar conditions (Reference Jones and JamiesonJones and Jamieson, 2001).
Classification tree using all available predictor variables (N = 411). Gray boxes show terminal nodes. The value in the boxes is the predicted DSI class. The number above the box is the number of days to be split.
We used discriminant-based univariate splits to determine the best splitting variable. For each predictor variable a significance level (p value) is calculated and the variable with the lowest p value is used to split the cases in a node into two subsets.
The right-sized tree should be as simple as possible, but should be sufficiently complex to account for patterns in the data likely to occur in similar datasets. We used minimal cost–complexity cross-validation pruning, since our limited dataset did not allow us to exclude an independent test sample of adequate size. Here the originally grown tree is pruned upwards by deleting branches of the tree until the child node is reached. We used a three-fold cross-validation where the dataset is split into three samples of which two in turn are used to calculate auxiliary classification trees and the third sample is used as a test sample to predict the accuracy of the computed classification trees. For each of the subtrees the cross-validation costs and errors are calculated following Reference Breiman, Friedman, Olshen and StoneBreiman and others (1984, p. 66). The subtree with the lowest classification cost is thought to be the right-sized tree.
The misclassification cost is generally the proportion of misclassified cases. The advantage of calculating the misclassification cost instead of looking at the misclassified cases is that it is possible to consider that a more accurate prediction for one class is sometimes desired. Setting priors on each class influences the cost calculation and with this the development of the tree. Our focus is on correctly classifying instability (DSI = 1) more often than stability (DSI = 0). Considering this and the unbalanceddataset between classes (stable days = 304, unstable days = 107), we changed the default from estimated priors, whichareproportional to class size, to equal priors to account for the desire to forecast unstable days more correctly. Setting an even higher prior on unstable days, which tends to decrease its misclassification rate (Reference Breiman, Friedman, Olshen and StoneBreiman and others, 1984, p. 112), resulted in a questionable tree structure.
As mentioned, cross-validations use the number of misclassified cases to estimate the probability of misclassifying a case (Reference Breiman, Friedman, Olshen and StoneBreiman and others, 1984, p. 73). We applied the global cross-validation function to test the tree performance. The classification-tree analysis is replicated a specific number of times, in our case three times, each time holding out a fraction (one-third) of the dataset, which is used as a test sample to cross-validate the classification tree. In the end, each case is tested once in a classification tree. The misclassification cost and error are calculated, giving an indication of the tree performance.
Further we looked at the importance ranking of the predictor variables. Here all predictor variables are ranked in accordance with their potential effect on the classification (Reference Breiman, Friedman, Olshen and StoneBreiman and others, 1984, p. 41, 147). This is necessary since the final tree structure alone is often not indicative of the potential importance of the predictor variables.
Analysis
Spearman rank correlations
The results of the Spearman rank correlations are listed in Table 2. Significant correlations (p < 0.05) are highlighted in bold. Table 2 also shows the Spearman correlation coefficient R.
Rank correlations of predictor variables with daily instability index (DSI). Significant correlations are in bold
The physical explanations for the significant correlations with DSI are as follows:
Positive correlations
Humidity: Higher values of RH5, RHmn Y and RHmx Y are associated with stormy weather (precipitation) and slab formation (Reference McClung and SchaererMcClung and Schaerer, 1993, p. 161) and consequently with avalanche activity and therefore instability.
Precipitation: Higher values of precipitation (Pcp Y), new-snow (HNY) and storm-snow (Strm) accumulations indicate recent loading on the weak layers.
Wind: When wind speed (WSa) is greater, more snow is transported, loading the weak layer faster, and the formation of wind slabs is promoted (Reference McClung and SchaererMcClung and Schaerer, 1993, p. 157–158).
Skier-triggered dry-slab avalanches on previous day: Avalanches often occur on consecutive days because loading due to a snowstorm often continues for several days and because persistent weak layers require days to adjust to increased overlying load (Reference ChalmersChalmers, 2001, p.79–84).
Negative correlations
H: Having a negative correlation implies that a thinner slab thickness is associated with skier-triggered dry-slab avalanches. Skier stress generally decreases with increasing depth (Reference FöhnFöhn, 1987; Reference Schweizer and CamponovoSchweizer and Camponovo, 2001). Consequently, deeper weak layers are less often triggered by skiers. Reference FöhnFöhn (1987) found that since the calculated static skier stress at 1m is only 10% of the stress due to gravity acting on the slab, skiers are not efficient triggers where the slab is >1m thick. Aswell Johnson (Reference Johnson2000, p. 57) found a positive correlation of slab thickness with the strength of faceted layers, implying that deeply buried weak layers of faceted crystals are usually stronger.
Load: The negative correlation can be interpreted in the same way as slab thickness. Load is the primary snow-pack variable that affects the shear strength of layers of faceted crystals. Thicker slabs typically overlie stronger facet layers (Reference JohnsonJohnson, 2000, p. 57; Reference Zeidler and JamiesonZeidler and Jamieson, 2002). Additional slab load causes increased densification (Reference KojimaKojima, 1967; Reference Conway and WilbourConway and Wilbour, 1999) and increased bonding.
Sk38: This result was expected since lower values of Sk38 indicate lower stability and increased probability of skier triggering (Reference JamiesonJamieson, 1995, p. 148–158, 215–221; Reference Jamieson and JohnstonJamieson and Johnston, 1998).
HS: Winters with a deeper snowpack typically have less clear weather, and consequently the surface hoar layers that form typically consist of smaller crystals, which gain strength faster (Reference Jamieson and JohnstonJamieson and Johnston, 1998). Also, the slabs that bury those persistent weak layers are likely thicker sooner. Consequently, the stress induced by skiers might not penetrate deeply enough to cause the weak layer to fail. When the snowpack is relatively thin, densification takes longer and weaknesses are preserved for a longer time.
Even though temperature and wind are thought to be important factors in avalanche forecasting, we did not find significant correlations in our dataset. For temperature this may be because we focus on dry-slab avalanches from December to March, or because only daily temperature variables were considered, when hourly changes might be more relevant. The non-significant correlations for ridge-top winds (WS5, Wrun Y) may be due to the relatively sheltered location of the Columbia Mountains.
Classification tree
A classification tree was developed once using all variables listed inTable 1 and once using all variables except Sk38, H and Load. Figure 1 shows the tree developed using all variables.
The tree consists of three splits and four terminal nodes. Two splits are determined by Sk38 and one by load, suggesting that Sk38 is the most predictive factor, followed by Load, to classify DSI. When excluding Sk38, Load and H; the tree was developed mainly using the height of the snowpack and storm-snow accumulation.
Table 3 summarizes the global cross-validation results of the tree analysis. Sk38, Load and H improved the classification of both unstable days (DSI = 1) and stable days (DSI = 0). The proportion of misclassified cases and the misclassification cost were reduced by including the three snowpack properties.
Global cross-validation results from classification tree (N = 411)
The importance ranking shows that Sk38, Load, HS, the number of previous triggered avalanches and H are the most important factors for forecasting skier instability (DSI) in the area near Blue River.
Including Sk38, Load and H into a nearest-neighbour forecasting model for the Columbia Mountains near Blue River is justified based on the results of the rank correlations and the classification-tree analysis.
Incorporation of Snowpack Variables into a Nearest-Neighbour Model
We used several different combinations of meteorological predictor variables to find the best fit for the nearest-neighbour model, including variables found by statistical methods, all available variables and variables found in prior studies based on conventional avalanche forecasting. Two approaches yielded the best performance:
Case 1: (a) All variables including Sk38, Load and H; (b) All variables without Sk38, Load and H;
Case 2: (a) Variables from Spearman rank correlation including Sk38, Load and H ; (b) Variables from Spearman rank correlation excluding Sk38, Load and H.
The forecast of Cornice is said to be correct on a particular day if three or more of the most similar days have DSI = 1 and DSI = 1 on the forecast day, or if less than three of the most similar days have DSI = 1 and DSI = 0 on the forecast day. Table 4 summarizes the results. Due to missing input data, 37 of 515 days were not forecasted.
Performance of nearest-neighbour model
Sk38, H and Load were all included since the results improved by an additional 2% when using all three variables instead of using Sk38 alone.
In both cases, including Sk38, H and Load improved the performance of the nearest-neighbour model.
Further we looked at misclassified days to understand how the model could further be improved. As an example we chose a time series of a surface hoar layer buried on 16 February 2002, and compared case (2a) with case (2b). Sk38 values were calculated for 38 days.
Incorporating Sk38, H and Load improved the number of correctly classified unstable days on the 16 February 2002 surface hoar layer (Table 5). The number of correctly fore-casted stable days stayed the same. In this small sample, false stable misclassifications were reduced from four to one.
The model without Sk38, H and Load misclassified more unstable days than it classified correctly. This suggests that these three excluded variables are important predictors of skier-triggered avalanches in a nearest-neighbour model.
Summary
In this study, three predictor variables, Sk38, H and Load, were explored for their potential to predict skier-triggered dry-slab avalanches on persistent weak layers and included into a nearest-neighbour model. We defined a daily skier instability index as a response variable to be able to include days when triggering of dry-slab avalanches on a persistent weak layer was likely, but not reported.
Spearman rank correlations were applied to assess the relative importance of each predictor variable individually. Humidity, precipitation, wind, previous triggered avalanches, slab thickness, slab load, Sk38 and height of snow-pack were the most significant. The physical relationship of the significant predictors to instability was discussed.
Classification-tree analysis allowed us to assess the importance of the stability and slab properties of persistent weak layers observed in a centrally located study plot, but did not identify an optimal set of predictor variables that could be used in a nearest-neighbour model. The importance ranking placed Sk38 as the most important variable when forecasting daily instability for skier-triggered dry slabs on persistent weak layers in the Columbia Mountains. Slab load and slab thickness were ranked second and fifth in importance, respectively, of all variables considered for potentially predicting unstable days.
Performance of Cornice for 16 February 2002 surface hoar
Cornice, a nearest-neighbour model developed in Scotland, was configured using two sets of predictor variables. Both sets were run once with and once without Sk38, Load and H. In both cases the forecast was improved by including the snowpack properties.
It is more important to correctly forecast unstable days since false stable predictions have greater consequences than false unstable predictions. The improvements for correctly forecasting unstable days were 7% for the case of using all available predictor variables and 3% when using variables that were significant in the Spearman rank correlation analysis. Looking at a test series, the first week after burial with low strength measurements was forecasted more accurately using Sk38, Load and H than when excluding these parameters.
Acknowledgements
For discussions on computer-assisted forecasting, we thank R. Purves, J. Schweizer and M. Gassner. For their careful snowpack measurements, we are grateful to many people including J. Hughes, L. Allison, A. Cooperman, O. David, K. Black, R. Gallagher, N. Irving, P. Langevin, J. Shepherd andJ. Olson. Our thanks to Mike Wiegele Helicopter Skiing for providing weather data, avalanche occurrence reports, logistical support and a stimulating environment for field studies. A. Zeidler is sponsored by the Gottlieb Daimler- und Karl Benz-Stiftung/Germany. This paper is a result of an ongoing University–Industry project funded by the British Columbia Helicopter and Snowcat Skiing Operators Association, the Natural Sciences and Engineering Research Council of Canada, Canada West Ski Areas Association and the Canadian Avalanche Association.
