Avalanche runout distance and return period estimates are essential to snow avalanche risk assessment for hazard mapping and mitigation design. We present a validation of the Space-Time (ST) model, a statistical model that expresses runout distance as a function of return period. The validation is based on field observations of tree damage and tree ring data of infrequent (1:5 to 1:100 year return period) avalanches from 34 paths in the Canadian Rocky and Columbia Mountains. While the ST model has been applied successfully for a longer return period (>100 year) runouts in one path, our validation showed that it does not independently estimate infrequent runout distances with sufficient accuracy for hazard mapping and mitigation design. As with estimating extreme runouts, it should be used in combination with other methods. The model was found to perform differently across mountain ranges, and tended to estimate runout distance downslope of those observed in the field.

## Introduction

Snow avalanches threaten people and infrastructure in mountain regions worldwide. Between 100 and 200 fatalities occur per year in Europe and North America (Hadmoko and Mauro, 2012). Fatalities also occur in South America, China, Japan, New Zealand, Afghanistan, and Iran. Data on annual worldwide costs of avalanche mitigation and damage are not systematically assembled. In Canada, mitigation, damage, and downstream costs exceed $100 million per year (Stethem et al., 2003). In Europe, the costs are substantially higher because of dense populations in mountain regions.

Avalanche risk assessment and mitigation relies on estimates of avalanche magnitude and frequency. In this paper, magnitude is described only by the distance the avalanche runs (runout). Other definitions of magnitude can include destructive potential, size of deposit, and impact pressure. Frequency can be described by the number of avalanches that reach a certain distance over a given period of time. Both frequency and magnitude typically decrease as one moves down the centerline of an avalanche path. Standard practice for determining avalanche magnitude and frequency along a path includes surveys of vegetation damage, analysis of historical occurrence records and aerial photographs, terrain analysis using digital elevation models, discussion with local experts, and the application of dynamic and statistical models. Each method contributes to the final risk assessment. This study focused on the Space-Time model (McClung, 2000), a statistical model that estimates runout distance (Space) as a function of avalanche return period (Time) and vice versa.

The Space-Time (ST) model has many potential benefits. First, it is simple and based on a widely used statistical model, the Runout Ratio model. The Runout Ratio model is a statistical model that estimates runout distances based on empirical evidence of a set of extreme runouts for a particular mountain range. Both the Runout Ratio, and another closely related model, the α-β model, estimate the runout point (α point) using its relationship with the point where the slope first reduces to 10° (β point). Second, it is based on probability, meaning that consultants can calculate runout distance based on likely non-exceedance limits. Third, the model can calculate runout distance for different frequency of avalanches; both frequent (<5 year), infrequent (5- to 100-year), and more extreme (>100-year return periods). Avalanches with these different frequencies are undetermined by previous models such as the α-β and Runout Ratio, which estimate extreme (100- to 300-year) runout distance only (Bakkehøi et al., 1983; Lied and Bakkehøi, 1980; McClung and Lied, 1986). Estimates of 10- and 30-year runout distances are useful for mitigating risk to infrastructure and development within the avalanche path, such as ski lift towers, industrial roads, residential development, and even recreational route planning (Canadian Avalanche Association, 2002; Statham et al., 2006). Estimates of 1000-year runouts are sometimes part of wider geohazard analyses.

McClung (2000) successfully applied the ST model to a path in Bleie, Norway, to calculate return periods from 50 years to 2000 years, but its application to a wider set of paths has yet to be established. The goal of this study was to validate the ST model using a wide range of paths, thereby improving its value for engineering risk assessment and mitigation. The validation was focused on infrequent events (∼1:5 to 1:100 year return), since the 10- and 30-year runout distances appear often in engineering risk mitigation applications. Note that very frequent events (<1 year return) were excluded from this study. The objectives were to (1) compare modeled runout distances with field-observed runout distances for infrequent avalanches, and (2) assess model performance across two mountain ranges.

### The Space-Time Model

Estimating the magnitude and frequency of a physical phenomenon is a known problem in the geohazards field. A widely used statistical approach is to apply a Peak Over Threshold model, which combines one distribution that models event magnitude with another that models event frequency. Events smaller than a given magnitude (threshold) are excluded (Benjamin and Cornell, 1970). In the case of the Space-Time model, runout distance is modeled by a Gumbel distribution and avalanche arrival rate is modeled by a Poisson distribution. The following paragraphs describe these distributions in relation to avalanche behavior.

#### Runout Distance (Space)

The runout component of the ST model is based on the Runout Ratio model (Fig. 1) (McClung and Lied, 1986; McClung et al., 1989; McClung and Mears, 1991).

The Runout Ratio model relates the extreme (100- to 300-year) runout position (α point) with the point where the avalanche path slope first reduces to 10° during descent (β point). This relation is given by the dimensionless ratio, *r*, where

*X*

_{β}= horizontal distance from the top of the start zone to the β point, and Δ

*x*= horizontal distance between the β point and the α point.

Since *r* is assumed to be Gumbel distributed, the probability of not exceeding a given value of *r* is given by the Gumbel cumulative distribution function:

*r*= given runout ratio;

*u*= Gumbel location parameter; and

*b*= Gumbel scale parameter.

Since the Gumbel parameters have been shown to vary with location, the Runout Ratio model has been calibrated to different mountain ranges in North America, Norway, and Japan (e.g., Fujisawa et al., 1993; Johnston et al., 2012). The variation is thought to reflect terrain shape rather than climate regime (McClung et al., 1989; McClung and Mears, 1991), because extreme avalanches are expected to share similar characteristics (large and dry) regardless of whether they are in a maritime, continental, or transitional climate.

Calibration involves collecting α and β point locations for a set of paths within a particular mountain range through field work or inspection of maps and aerial photos. The non-exceedance probabilities of each path are ordered, and then plotted using Hazen, Weibull, or a custom plotting position (e.g., McClung and Mears, 1991). The Gumbel parameters *u* and *b* are then determined by fitting to the ordered positions.

#### Arrival Rate (Time)

Avalanche arrival rate (also referred to as frequency) at a specified position is described by the Poisson probability density function:

where µ_{0}= expected annual arrival rate = 1/T

_{0}at the specified point;

*T*

_{0}= return period at the specified point (years); and

*n*= number of random events.

For example, the probability of one (*n* = 1) avalanche reaching the specified position with an expected return period of 10 years (*T*_{0} = 10) occurring in the next winter is given by *P*(0.1, 1) = 0.09.

The goodness of fit of avalanche frequency to the Poisson distribution has been tested in several studies. Föhn (1978) found avalanche arrival rate over 420 years in a single, unforested, low-frequency path fit a Poisson distribution (χ^{2} test, sig. level = 0.01). Smith and McClung (1997) applied Normal and Poisson distributions to 24 years of occurrence data from 43 high-frequency paths in Rogers Pass, British Columbia. The χ^{2} tests showed a better fit to the Normal distribution than the Poisson, particularly for higher frequency paths. Despite this, the Poisson was identified as the preferred distribution since it matched the conditions for a Poisson experiment (discrete, independent events in a continuous time period, each with an equal probability of occurrence).

#### Compound Distribution

When the two distributions are combined, the ST model can calculate either the runout distance for a given return period, or the return period for a given runout distance. In practice, inputs include:

known return period at the β point (

*T*_{0});horizontal distance from the start zone to the β point (

*X*_{β});Gumbel parameters appropriate for the mountain range (

*u*and*b*);return period(s) where the runout distance is desired or known (

*T*); and/or_{X}runout distance(s) where the return period is desired or known (

*X*)._{T}

The probability of an avalanche running a certain distance with a certain return period (Space-Time) is given by the compound cumulative distribution function (McClung, 2000):

where*P(*µ

_{o},

*n)*represents the Time component (see Equation 3), and [

*F(r)n*] represents the Space component (see Equation 2).

Although McClung (2000) chose a Gumbel distribution to represent runout (McClung and Mears, 1991) and a Poisson distribution to represent frequency (Smith and McClung, 1997), other distributions such as the Generalized Pareto Distribution could be used as a probability density function for extreme runout and would be theoretically preferable (Keylock, 2005).

The model works most simply when the return period at the β point is known; however, if the return period at β is unknown, the return period at a reference point (RP) further upslope or downslope can be used. This requires an adjustment to the location parameter, *u*, as follows:

*r*

_{RP}= (

*X*

_{RP}-

*X*

_{β})

*X*

_{β}

^{-1}is the Runout Ratio at RP, and

*X*

_{RP}= horizontal distance from the start zone to RP.

Note that other probabilistic models have been developed to describe avalanche return period and runout distance within a path (e.g., Harbitz et al., 2001; Barbolini and Keylock, 2002). Recently Bayesian approaches have been used to quantify uncertainty caused by the lack of site-specific data and limitations of dynamic and empirical models (e.g., Eckert et al., 2008; Cappabianca et al., 2008). Although these are valuable developments, this paper focuses on validation of the ST model only.

The Appendix contains two worked examples of the ST model. The first example calculates the return period at a given distance (*T _{x}*) when the return period at the β point is known. The second repeats the calculation for when the return period is only known at a reference point, RP.

## Methods

To validate the ST model, we compared ST model estimated runout distances against those observed in the field for 34 paths in the Canadian Rocky and Columbia Mountain Ranges.

The objective of the field validation was to establish a data set of avalanche paths (*n* = 34) where both the return period (*T*) and runout distance (*X*) were known for several key points along each path based on terrain and vegetation observations. The observed variables included *X*_{β}, *X*_{RP} (if applicable), *T*_{0}, and a series of runout distances (*X*) with return periods (*T*) inferred from vegetation damage.

This validation was undertaken in three stages. First, a review was conducted to select paths, identify trim lines, and estimate return periods based on aerial photographs and satellite imagery. This was followed by a field assessment, and finally a comparison of field assessments with model output.

### Path Selection and Pre-assessment

Over 100 avalanche paths within British Columbia and Alberta were inspected using digital elevation models (DEMs), topographic maps, aerial photographs, and satellite imagery. The initial criterion for identifying paths was to have good vehicle or walking access from main highways or roads. A number of the 100 paths (or path groups) were suggested by local experts, who perform avalanche control work in the regions and who have extensive backcountry knowledge of the mountain ranges. Other paths were suggested by avalanche researchers, who had visited the areas for earlier research projects. The remaining paths were found by inspecting drainages that intersected or ran parallel to main highways or roads, through Google Earth. Of these, 34 paths were selected for field-survey based on criteria that excluded paths with runup (avalanches running up the other side of the valley) and short vertical drop (<350 m). These two characteristics substantially affect avalanche behavior; runup results in shorter avalanche runout by reducing momentum, and paths with short vertical drop exhibit runout ratios that are different from taller paths due to length scale effects (Jones and Jamieson, 2004).

It was important to select paths that were not exposed to explosive avalanche control or did not contain defense structures because these human interferences affect avalanche frequency and/or runout distance. Since avalanche control is typically performed to protect roads and railways, paths were chosen from valleys behind the first line of mountains directly bordering these transportation corridors. All paths were accessed by foot.

Since limited historical records were available to determine return periods at given runout distances, we relied on vegetation boundaries (trim lines) to provide evidence of past avalanches. Trim lines were identified as boundaries between vegetation of constant height or constant age, which regrew after being destroyed by an avalanche (Fig. 2). Trim lines were initially identified using historical aerial photographs and satellite imagery. We did not attempt to field-truth the trim lines identified in the office, but rather selected the paths with three or more distinctive trim lines in preference to those without. This helped to select paths where trim lines were most likely to be easily observed in the field. In order to focus on larger events with sufficient snow to damage trees (Germain et al., 2010), only evidence resulting from avalanches greater than Canadian Destructive Size 1 were considered (Canadian Avalanche Association, 2014). Figure 3 presents a map of western Canada showing the study sites. Table 1 shows the approximate GPS coordinates for each study site (path).

This pre-assessment stage was key to reducing the number of samples of vegetation required to establish return periods in the field.

### Field Assessment

As avalanches destroy or damage sections of forest, they leave behind changes in vegetation patterns (trim lines). These trim lines reveal changes in return periods down the length of the path. During field work undertaken between May and September 2012, key trim lines were identified and recorded in the field using both vegetation and terrain observations according to the following method.

#### Field Pre-assessment

A path-wide assessment was undertaken from a vantage point. The assessment involved matching path characteristics (roads, creeks, trim lines) to the maps, aerial photographs, satellite imagery prepared during pre-assessment (see Path Selection and Pre-assessment), and notes from conversations with local guides and avalanche experts. The objective was to develop initial (rough) estimates of the

beta point location and return period range,

reference point location and return period range (if applicable), and

three to six key trim lines and associated return periods.

Return period ranges (e.g., 10 to 30 years or 100 to 300 years) were estimated based on path-scale observations of vegetation patterns. These included observations of stands that shared a similar tree height, age, species, and stem density, having regrown since being cleared by previous avalanche events. These vegetation characteristics were first identified by Schaerer (1972) and then by McClung and Schaerer (2006) as a method for establishing a return period range (Table 2). In many cases, these patterns were more discernible from afar, than from within the path. After initial assessment, a plan was made to target particular zones along the centerline of the path.

## TABLE 1

Coordinates for each study site, recorded as UTM Zone 11 eastings and northings.

#### Detailed Survey—Terrain Observations within Path

A segmented survey of the runout was conducted for each path. The survey included walking the likely centerline of the path that a large, dry avalanche would follow through the runout zone, and recording the incline, slope distance, elevation, and coordinates of each segment. Segments ranged from 10 m to 90 m depending on visibility through vegetation and terrain shape. Vegetation observations were also made along the centerline and are explained further in the next section.

In determining the point where the path slope first reduces to 10° (β point), characteristics of each individual path were considered. For example, large avalanches tend to run farther into the runout when the path has a confined track (Bakkehøi et al., 1983; McClung, 2003). In these cases, benches high in the runout zone may not have as much effect on the runout distance as paths with non-confined or planar tracks. In all cases, benches in the track and benches shorter than 3% of *X*_{β} were ignored because they were assumed to have negligible effect on avalanche flow when compared with overall avalanche size.

Each β point was confirmed with the highest accuracy and resolution DEM available, which ranged from 1 m to 15 m. For paths where discrepancy existed between field and DEM methods, judgment was used to determine an appropriate β point (Sinickas and Jamieson, 2014).

Equipment included an inclinometer (accuracy of ±0.5°), hip chain (accuracy of ±2%) (Johnston et al., 2012), and Global Positioning System (GPS) receiver. Two different GPS units were used; a Garmin GPSMap 62s (average horizontal accuracy ±7 m) and—for a limited number of paths—a NovAtel DL-V3 (average accuracy ±<1 m). Multiple waypoints were recorded when stopped between segments using the Garmin GPSMap 62s to allow for averaging during analysis. Elevation values were recorded using an aneroid altimeter (accuracy ±5 m) calibrated on the day of the survey, which provided more reliable data than elevations recorded by the GPS units.

## TABLE 2

Examples of vegetation as an indicator of avalanche frequency. After McClung and Schaerer (2006).

#### Detailed Survey—Vegetation Observations within Path

Sampling was designed to estimate changes in return period at each trim line, identified during the field pre-assessment, down the centerline of the avalanche path. The first objective was to isolate avalanche-related damage from damage caused by other traumatic events such as snow creep, debris flow, animal damage, infestation, logging, or fire. The second objective was to balance a sufficient number of samples from within each path to estimate return period ranges, with a sufficient number of paths across the mountain range to validate the ST model.

Two types of observations were gathered:

*External*evidence of avalanches (Fig. 4), which included observations of scars, split trunks, trees bent in the direction of avalanche flow, decapitation of trees, elimination of branches, or trees bent in the direction of avalanche flow and uninjured survivor trees (Stoffel et al., 2010). These observations were recorded in addition to initial observations regarding common tree height, age, and density.*Internal*evidence of avalanches (Fig. 5), which refers to annual tree rings and to reaction wood. Tree rings indicate tree age, where the combination of one dark and one light ring represents one year of growth. Reaction wood develops when the upper growth of a damaged tree veers back toward vertical growth after a traumatic event such as an avalanche, snow creep, animal damage, or fire. Since reaction wood occurs within a year or so following an avalanche event, the approximate date of each event can be determined in the field by counting tree rings. An average return period can then be calculated by dividing the age of the tree by the number of avalanche occurrences in its lifetime. For example, a 50-year-old tree that shows evidence of five avalanche events has an approximate average return period of 10 years. Because avalanches decrease in frequency with distance down the path, reaction wood becomes less common, and average return periods increase. Paths can have abrupt increases in return period depending on the terrain and existing vegetation. These changes in return period were typically identified during the field pre-assessment step.

To expose internal evidence, we took cross-section slices (for trees with diameters of less than approximately 15 cm) and increment cores (for trees with diameters of more than 15 cm, or trees within National or Provincial Parks). In order to isolate avalanche-related events, from other events, we did the following:

Targeted avalanche paths that had well-defined tracks and runout zones during the pre-assessment (Path Selection and Pre-assessment).

Estimated return period ranges based on aerial photographs and satellite imagery during the pre-assessment (Path Selection and Pre-assessment).

Targeted specific lines down the centerline of the path, as identified during the pre-assessment (Path Selection and Pre-assessment), and the assessment conducted from a vantage point at the start of each field survey (Field Pre-assessment).

Targeted trees that best represented a group of trees bent in the direction of avalanche flow, or that shared scars on the uphill side of the tree. These tree types have been shown to contain the most reliable evidence of past avalanches (Germain et al., 2010; Johnson, 1987).

Recorded an average of 6.7 observations per trim line of external and/or internal evidence of past avalanches.

Noted, in particular, any evidence of logging or human activity from oral histories, or from field-based evidence. In all cases, where evidence of logging was encountered, it was in the lateral boundaries of the runout zone, not in the centerline of the avalanche path.

Noted in particular, any evidence of animal damage (e.g., bear scratching). In all cases, animal damage was observed outside of the centerline of the avalanche path, not accompanied by other signs of avalanche damage and isolated to one or two trees, rather than many trees at a trim line.

Excluded paths that showed evidence of debris flow along the centerline of avalanche flow (Luckman, 2010). These paths were removed from the data set due to insufficient detail and time to sample, map, and separate the damage caused by avalanches from debris flow.

Sampling methods were designed to target avalanche-related damage (see above), and to be as efficient as possible in estimating return period ranges down the centerline of each path. Both the minimum sample size to establish an avalanche event (tree-scale) and the minimum sample size to characterize avalanche activity within a track (path-scale) vary in the literature.

The suggested minimum percentage of trees exhibiting the same evidence required to establish an avalanche event ranges from 10% (Butler and Sawyer, 2008) to 40% (Reardon et al, 2008). The suggested minimum number of trees to establish avalanche activity across a path ranges from 10 trees (Butler and Sawyer, 2008) to 100 (Corona et al., 2012).

Such a large sample size was impractical for this project because of the tradeoff between the number of observations within a path and the number of paths within the study. On average, one path was surveyed each field day, and 34 paths were surveyed in total. In total, 859 external observations were recorded, and 245 internal (163 slices and 82 cores) observations were recorded across the 34 paths (Table 3). Between 15 and 46 observations were collected per path, which is in line with Germain et al. (2010), who suggested that diminishing returns to determining avalanche chronologies at a sample size is >40. To achieve this, and collect sufficient quality data, we did the following:

Targeted coniferous species, which tend to show annual rings and damage more clearly than deciduous species (Luckman, 2010). Reaction wood caused by avalanche impact tends to develop on the downhill side in conifers, and the uphill side in deciduous trees (Burrows and Burrows, 1976).

Excluded reaction wood from the first five years since these smaller trees tended to exhibit reaction wood from simply growing through the snow on a slope.

Estimated return periods as ranges, rather than precise avalanche event dates. This allowed us to reflect uncertainty as a wider range. Ranges were calculated by averaging events over time, in different combinations. For example, if a tree cross section showed reaction wood from 60, 45, 32, 27, 10, and 3 years ago, then the return period would be calculated using:

This method was applied in the field as a quick check to identify extraordinary events if one range was much larger or smaller than another. In these cases, evidence was gathered to give surveyors more confidence in their observations and calculations. Wider ranges were recorded in cases where the extra evidence did not provide more certainty.Collected samples along the centerline of the avalanche track from top to bottom at the trim lines initially identified during the field pre-assessment. This way, we were able to track common events, and the location in the runout zone where these events did not appear to cause any lasting damage (scars, bent trees of similar age, leaders, reaction wood) in the trees.

Sample output from a typical field survey is shown in Figure 6.

### Comparison of Modeled and Field-Assessed Runout

The ST model was used to estimate the runout distance to each trim line (*X _{T}*) identified in the field (

*n*

_{trims}= 129 for 34 paths). Inputs for each path included the location of the β point (

*X*

_{β}), the reference point (

*X*

_{RP}), and their return periods (

*T*

_{0}). For each path, the highest trim line recorded in the field was chosen as the reference point. For 16 out of the 34 paths, the first trim line coincided with the β point so that

*X*

_{RP}=

*X*

_{β}.

Because the return period at *X*_{RP} was recorded as a range (*T*_{0,1}, *T*_{0,2}), the ST model runout distances were also estimated as a range. For example, a return period at the reference point of between 1 and 5 years might translate to a runout distance of between 1000 m and 1200 m at a given trim line. Descriptive statistics of the input values for the 34-path data set and subsets are shown in Table 4. Note that *T*_{0} refers to the return period at either at the β point or RP. The Gumbel parameters (*u*, *b*) were 0.079, 0.070 for the Rocky/Purcell ranges (McClung and Mears, 1991) and 0.105, 0.083 for the Columbia range (Johnston et al., 2012).

## TABLE 4

Descriptive statistics for the 34 paths.

#### Comparing the Field Survey to the ST Model

The ST model estimates of runout distance were compared with the field-observed runout distance (*X*_{T}*), through calculation of a residual (*X _{T}* –

*X**). Positive values represented scenarios where modeled runout distance (

_{T}*X*) was downslope of observed runout distance (

_{T}*X**), which means the model was too conservative. Although the model estimates runout distance as a range, the residual was calculated in favor of the model. For example, in Figure 7, a 30-year trim line is observed in the field at 1500 m horizontal distance from the top of the start zone (starting position). Because the ST model estimates the 30-year trim line somewhere between 1600 and 1800 m, the residual calculation is given by

_{T}*X*-

_{T}*X** = 1600 – 1500 = +100 m.

_{T}For scenarios where the observed runout (point) was within the modeled runout (range), the residual was recorded as zero.

This method minimized the residual (i.e., used 100 m rather than 300 m in Fig. 7) to simplify the validation process; if the model was shown to be inadequate even in its best light, then no further analysis would be required.

Residuals were scaled by *X*_{β} to account for different path lengths using (*X _{T}
* -

*X**)

_{T}*X*

_{β}

^{-1}. Note that residuals were small enough compared with

*X*

_{β}that error in determining

*X*

_{β}had negligible effects on the results (Sinickas, 2013). Residuals were then grouped into <10 year, 10 to 20 year, 20 to 50 year, 50 to 100 year, and 100 to 300 year trim lines to reflect return periods typically used in avalanche risk analysis. Ranges increase at higher return periods with the scale effect of uncertainty (higher uncertainty at longer return periods). This provides a larger chance for the model to be correct.

## Results

Box and whisker plots of scaled residuals for all paths are shown in Figure 8, and then divided into the Rocky/Purcell range and the Columbia range in Figure 9.

Plots can be interpreted by visualizing an avalanche path with the start zone above the top of the plot, and the runout zone at the bottom. Note that the vertical axis has been reversed so that positive residual values appear in the lower portion of the plot. Perfect agreement between observed and modeled runout distances lies on the zero line. As an example, for a path where *X*_{β} = 1000 m, a residual of +0.2 represents a scenario where the modeled runout distance (*X*_{T}) is 200 m farther downslope of the observed runout distance (*X*T^{*}).

Three findings arise from the results. First, median residual values are positive in both Figure 8 and Figure 9, meaning that modeled runouts are longer than observed runouts. Second, model performance varies across mountain ranges. In Figure 9, the Rocky/Purcell scaled residuals are significantly smaller than the Columbia residuals for all return period ranges (see Table 5 for significance test results); Rocky/Purcell medians are limited to less than 0.1, while Columbia medians range from 0.17 to 0.32 for *T* < 100 years. Third, the interquartile range remains relatively constant with increasing return periods in Figure 9.

Figure 10 shows scaled residuals calculated for paths where the *T*_{0} was observed at the β point (simple model), and where *T*_{0} was observed at a reference point away from the β point (adjusted model).

Figure 11 shows the Q-Q plot for path steepness, represented by the β angle. In this data set, however, there was no significant difference in the residuals between flatter and steeper paths (*t*_{stat} = 1.6 < *t*_{crit} = 1.7, *p* = 0.05).

## Discussion

The distribution of residuals presented in Figure 8 and Figure 9 implies that the ST model is generating longer runout distances (*X _{T}*) than what was observed in the field

*(X**). The distance (or residual) between the modeled and observed runout (

_{T}*X*- X

_{T}_{T}*) is scaled by

*X*

_{β}to allow for comparison across paths with different lengths, where

*X*

_{β}represents the horizontal distance between the start zone and the β point. Consistently positive medians for the ratio (

*X*-

_{T}*X**)

_{T}*X*

_{β}

^{-1}imply systematic errors occurred during the validation process, in the ST model, or in a combination of both.

## TABLE 5

Significance test results (Mann-Whitney U, 2-tailed) between residuals for Rocky / Purcell paths and Columbia paths.

### Validation Error

Systematic error in the validation could have arisen from primarily using vegetative damage to determine avalanche frequency. While analysis of changes in vegetation patterns and tree rings has been accepted as part of standard practice in North America, starting with Potter (1969) and Schaerer (1972), it is subject to high uncertainty (Reardon et al., 2008; Stoffel and Bollschweiler, 2006) for several reasons: (1) Vegetation is highly varied in itself. It can grow faster in some areas, and if subject to unseasonal temperature fluctuation can create false (intra-annual) growth rings, or skip growth rings entirely (Johnson, 1987). (2) Large avalanche events can destroy trees, which contain valuable evidence of previous events (Germain et al., 2009). (3) Not all avalanches leave evidence in the vegetation, particularly very frequent avalanches, which have limited lateral extension and can be more influenced by vegetation than extreme avalanches. (4) Vegetation can be damaged by natural and human events. Natural events include fires, wildlife, infestations, and extreme wind and ice events. Human events include forestry, farming, and avalanche control.

The methods described were chosen to balance this uncertainty with collecting observations from a sufficient number of paths to validate the model:

A minimum avalanche size was established.

Avalanche paths with well-defined tracks and runout zones were targeted.

Vegetation observations were targeted (skilled) rather than performed in transect (unskilled). They were collected along the centerline of avalanche flow, at key trim lines identified during the field pre-assessment that matched aerial and satellite imagery, local knowledge, and terrain shape to the path. In all cases, internal evidence was collected based on prior observations of external evidence.

Paths known to be affected along the centerline by human activity (mining, avalanche control, logging) or other natural hazards (bushfire, debris flow, animal scratching) were excluded from the data set.

Return periods were recorded as ranges rather than a point value, to reflect uncertainty.

While this study applied additional controls to the standard field methods, avalanche frequency could have been underestimated.

There was a tradeoff between collecting sufficient observations within each path and surveying a sufficient number of paths to validate the model. For paths in the Columbia data set, we recorded 295 external observations, 104 slices, and 26 cores over 14 paths. For paths in the Rocky/Purcell data set, we recorded 564 external observations, 59 slices, and 56 cores over 20 paths. Fewer slices were sampled within the Rockies/Purcell data set because 14 of the paths were within provincial and national park boundaries, where only cores were permitted.

Suggested sample sizes for establishing avalanche activity are varied in the literature, ranging from ten trees (Butler and Sawyer, 2008) to 100 (Corona et al., 2012). Our internal (cores and slices) observations were at the lower end of this range (between 6 and 19) for paths outside of the national and provincial parks, and below this range (between 2 and 5) for paths within the parks, where we were not permitted to slice trees. We relied on external evidence, rather than internal evidence for these paths (397 observations over 14 paths in parks). We expected to see such a small sample size of internal evidence to affect the validation results in two ways. First, if we had missed identifying avalanche events due to insufficient sample size, we would expect to see underestimation of avalanche frequency in both data sets, and therefore positive residuals (true *X _{T}** downslope of observed

*X**). This was observed in Figure 9, where both data sets exhibited positive residuals. Second, we would expect to see larger residuals in the Rocky/Purcell paths because fewer slices were sampled within this data set, potentially leading to larger underestimation of frequency. This was not observed in Figure 9, where Rocky/Purcell residuals were significantly smaller than Columbia residuals. Without historical records for comparison, it is difficult to estimate the influence of sampling design and size on the model validation.

_{T}The ST model may not be applicable to independently estimating shorter return periods (*T* < 100 years). In demonstrating the ST model, McClung (2000) calculated return periods up to 2000 years using reference return periods of 50 and 100 years for the Bleie path in Norway. In this validation, we calculated return periods from 2 to 300 years, based on reference return periods of 1 to 10 years.

For frequent (1 to 10 year) avalanches, total mass, controlled by entrained mass, has been shown to define runout distance (Steinkogler et al., 2014). It is expected that the ST model will lose relevance with shorter return periods, because the influence of snow supply for nonextreme conditions is not included in the statistical model. Since the ST model could be useful for estimating 10-year and 30-year events for engineering applications, and since the boundary between *frequent* and *extreme* is not clearly defined, we chose to validate the ST model for the following return periods: 2 to 10 years, 10 to 30 years, 30 to 50 years, 50 to 100 years, and 100 to 300 years.

### Model Error

Systematic error in the model could have been caused by several factors: The underlying assumption of a Poisson distributed arrival pattern may be inappropriate. In several field surveys, it was observed that avalanches could have “memory,” meaning that avalanche arrival may not be independent between winters. For example, in one path, an avalanche destroyed a 50-year-old stand of forest in 2008. By clearing the vegetation, a 2009 avalanche was able to travel farther than its predecessor, even though it was smaller in volume and destructive potential. The 2009 avalanche was therefore dependent on the 2008 avalanche and violates the requirements of a Poisson distribution. In the literature, avalanche dependence within years has been investigated (Smith and McClung, 1997) but not between years.

Avalanche runout distance is unlikely to be the only characteristic of avalanche magnitude. For example, one path in the data set showed evidence of an extreme event that had destroyed 300-year-old trees on the path periphery, but only traveled to the 30-year runout distance along the center-of-flow. In this case, although the avalanche was destructive and extreme, it was not adequately described by runout distance. The model may consistently underestimate the frequency of such extreme events, which could cause modeled runouts to be farther downslope of the true runout.

The adjustment of ST model's

*u*parameter (see Fig. 12 where*u*represents the peak of the Gumbel distribution) when using a reference point other than the β point has likely skewed the results. The simple version of the ST model is based on the initial return period*T*_{0}, at the β point. The more complex version allows the user to nominate the initial return period,*T*_{0}, at a different point,*RP*, by adjusting the*u*parameter according to Equation 5. This is a convenient and practical characteristic of the model (since the return period is not always easily defined at the β point), but is likely responsible for systematic overestimation of modeled runouts, particularly in the Columbia Mountains.Figure 10 shows the Columbia path residuals to be substantially smaller when using the simple, unadjusted form of the ST model. The difference may be due to parameter calibration for terrain shape and snow supply; however, since sample sizes were small for each time period (i.e., average number of samples was 4 or 5), more data are required to confirm the cause. Although not shown here, the difference between the simple and adjusted residuals for paths in the Rocky/Purcell Mountains is less pronounced, although paths using the simple ST model produced narrower interquartile and total ranges.

Path steepness may dominate runout estimates. The Runout Ratio, upon which the ST model is based, has been shown to be conservative for flat runout zones and short slopes (McClung, 2001; Jones and Jamieson, 2004). Johnson (1987) also investigated path steepness and found that the return period increased more slowly down paths with steeper slopes than flatter ones. In this data set however, there was no significant difference (

*t*_{stat}= 1.6 <*t*_{crit}= 1.7,*p*= 0.05) in the residuals between flatter and steeper paths, where path steepness, represented by the β angle, was assumed to be normally distributed (see Q-Q plot in Fig. 11).The

*u*and*b*parameters show the peak and spread respectively of the Gumbel distribution representing the Runout Ratio, which underpins the ST Model (Fig. 12). These are calibrated for extreme (>100 year) events. These parameters are expected to lose relevance when the model is applied to much shorter return periods (e.g., 2–10 years, 10–30 years), as snow supply starts to have a greater influence on runout distance.Last, the estimation of

*u*and*b*values for the Columbia and Rocky/Purcell Mountain Runout Ratio models were determined through extreme value analysis of field data (McClung and Mears, 1991; Johnston et al., 2012). These field data are affected by the factors described above.

### Performance across Mountain Ranges

The difference in ST model performance between the Columbia and Rocky/Purcell Mountains (Fig. 9) can be related to the calibration of the model parameters, *u* and *b*.

As mentioned above, a different *u* and *b* were derived for each mountain range through analysis of field observations. Figure 12 shows the Gumbel distribution, where *u* describes the location of the peak of the distribution and *b* describes the spread.

In his interpretation of the Gumbel parameters, McClung (2000) showed that *b* is dependent on terrain steepness, where flatter paths produce wider spreads, and that *u* is dependent on snow supply.

In this data set, there was no significant difference in path steepness between the Rocky/Purcell and Columbia Mountains (*t*_{stat} = 0.3 < *t*_{crit} = 1.8, *p* = 0.05), despite the different *b* values. This could explain the difference in ST model performance for the two regions. Further investigation of the paths used for calibration in McClung and Mears (1991) and Johnston et al. (2012) is required to determine whether this data set is unbiased.

McClung et al. (1989) suggested that regardless of snow climate (maritime or continental), large, dry avalanches will dominate the extreme runout distance over an observation period of 100 years. For frequent avalanches (1 to 10 years), total mass, controlled by entrained mass, has been shown to define runout distance based on five avalanches from a single path at the Vallée de la Sionne (VDLS) test site (Steinkogler et al., 2014). Sovilla et al. (2006) also showed that, for a path at the Mount Pizzac test site, and for the VDLS path, mass dominated runout for a data set of 12 events. The return periods for these events were not published; however, the data set included one extreme event, defined as the largest known avalanche. While snow cover affects runout distance at these shorter return periods for single paths (e.g., VDLS, Mount Pizzac) and is expected to change across different snow climates, there is limited literature on the effect of snow supply on avalanches with return periods between these two bounds (e.g., 10 to 100 year events) across many paths in multiple mountain ranges.

## Conclusion

Hazard mapping often requires estimates of runout distance for various return periods. Common estimates include the 10-year, 30-year, 100-year, and 300-year runouts. The ST model is a statistical model that estimates runout distance as a function of return period down the centerline of an avalanche path and vice versa. Although intended for estimating extreme events, this was the first study to assess its applicability to return periods of less than 100 years.

The goal of this study was to validate the ST model based on field observations of runout distances and return periods between 2 and 300 years. The primary objective was to compare modeled runout distances with those observed in the field. The secondary objective was to compare model performance across mountain ranges.

The results of this validation show that the ST Model does not estimate infrequent (1:5 to 1:100 year) runout distances with sufficient accuracy for independent application to hazard mapping and mitigation design. However, model estimates may be useful when combined with path-specific records (e.g., human records, trim lines) of avalanches with comparable return periods. Validation error was derived from using observations of vegetation damage to estimate infrequent return period ranges.

Two additional findings were made. First, modeled runout distances were consistently downslope of observed runout distances. This implies that the ST model is conservatively biased. However, the validation method (vegetative analysis) may systematically underestimate avalanche frequency, which would produce the same result. Second, median residuals were smaller for the Rocky/Purcell data set (3% to 11% of *X*_{β}), than the Columbia data set (17% to 35% of *X*_{β}). This implies that although the model adjusts for different terrain and climate characteristics, it does not perform equally across mountain ranges.

## Recommendations

Users cannot place confidence in the ST model as an isolated tool for hazard mapping and mitigation design of infrequent (1:5 to 1:100 year) avalanches. The ST model should only be used in combination with other path-specific estimates of runout or return period such as dynamic models, and analysis of vegetation damage and historical records. This is in line with how statistical models such as the α-β and Runout Ratio are used as one estimate for extreme avalanche runout.

When using the ST model to provide an estimate of runout period for infrequent avalanches, users should expect

conservative runout estimates,

estimate accuracy to change across mountain ranges, and

higher accuracy for the simple version of the ST model, where the initial return period is nominated at the β point, rather than at a different reference point.

To implement the ST model more widely, research could focus on applying the model to paths with consistent historical records. In this study, vegetative records were used in lieu of selecting paths with consistent long-term human records (chronicles) of runout.

Additional research could focus on the effect of using the return period at RP rather than at β on estimated runout distances. This would help to identify the source of bias, or lack thereof, in the ST model.

The calibration of *u* and *b* values also warrants further investigation. First, the parameters could be developed to include path-dependent adjustments based on terrain steepness and snow supply, rather than mountain range-wide adjustment. Second, the parameters could be developed to incorporate longer and shorter return periods than 100 years.

## Acknowledgments

Thank you to National Sciences and Engineering Research Council, Parks Canada, HeliCat Canada, the Canadian Avalanche Association, Mike Wiegele Helicopter Skiing, Canada West Ski Areas Association, Backcountry Lodges of British Columbia Association, the Association of Canadian Mountain Guides, Teck Mining Company, the Canadian Ski Guide Association, Backcountry Access, the B.C. Ministry of Transportation and Infrastructure Avalanche and Weather Programs, the Canadian Avalanche Foundation, and Tecterra for their support of the Applied Snow and Avalanche Research group at the University of Calgary. Special thanks to Chris Argue, Andrew Mason, and Christine Sinickas for volunteering to assist with fieldwork. Thank you to David McClung, Walter Steinkogler, Katherine Johnston, and Alan Jones for in-depth discussions of the model and research approaches. Thank you to Ryan Buhler, Mike Conlan, Shane Haladuick, Simon Horton, Michael Schirmer, and Scott Thumlert for assisting with fieldwork. Thanks to the following avalanche forecasters for sharing their local knowledge: Doug Wilson, Mike Koppang, Marc Deschenes, Phil Hein, Mark Vesely, Brad White, Jim Phillips, Dave Healey, and Jay Chrysafidis. Thanks to Doug Feely and Rick Kunelius for suggesting study sites, and Lisa Larson (Teck Coal Ltd.) for providing LiDAR DEMs of the Elk Valley.

## References Cited

*Journal of Glaciology*, 4: 24–29. Google Scholar

*Natural Hazards and Earth System Sciences*, 2(3/4): 239–245. Google Scholar

*Probability, Statistics, and Decisions for Civil Engineers*. New York: McGraw-Hill. Google Scholar

*Procedures for the Study of Snow Avalanche Chronology Using Growth Layers of Woody Plants*. Boulder, Colorado: Institute of Arctic and Alpine Research, University of Colorado, Occasional Paper 23. Google Scholar

*Natural Hazard Earth System Science*, 8: 303–309. Google Scholar

*Land Managers Guide to Snow Avalanche Hazards in Canada*. Jamieson, J. B. , Stethem, C. J. , Schaerer, P. A. , and McClung, D. M. (eds.). Revelstoke, British Columbia: Canadian Avalanche Association, 25 pp. Google Scholar

*Observation Guidelines and Recording Standards for Weather, Snowpack, and Avalanches*. Revelstoke, British Columbia: Canadian Avalanche Association, 90 pp. Google Scholar

*Cold Regions Science and Technology*, 54(3): 193–205. Google Scholar

*Cold Regions Science and Technology*, 74–75: 31–42. Google Scholar

*Stochastic Environmental Research and Risk Assessment*, 22(2): 185–206. Google Scholar

*In*Proceedings, Seminar of the International Union of Forestry Research Organizations, Davos, Switzerland, 241–254. Google Scholar

*Annals of Glaciology*, 18: 239–244. Google Scholar

*Climatic Change*, 92(1–2): 141–167. Google Scholar

*In*Stoffel, M. , Bollschweiler, M. , Butler, D. R. , and Luckman, B. H. (eds.),

*Tree Rings and Natural Hazards: A State-of-the-Art*. Heidelberg and New York: Springer, 51–73. Google Scholar

*Handbook of Hazards and Disaster Risk Reduction*. Milton Park , Oxfordshire: Routledge. Google Scholar

*Annals of Glaciology*, 32(1): 290–298. Google Scholar

*Ecology*, 68(1): 43–53. Google Scholar

*Canadian Geotechnical Journal*, 49(11): 1309–1318. Google Scholar

*Annals of Glaciology*, 38(1): 363–372. Google Scholar

*Cold Regions Science and Technology*, 42(3): 185–193. Google Scholar

*Journal of Glaciology*, 26(94): 165–177. Google Scholar

*In*Stoffel, M. , Bollschweiler, M. , Butler, D. R. , and Luckman, B. H. (eds.),

*Tree Rings and Natural Hazards: A State-of-the-Art*. Heidelberg and New York: Springer, 27–34. Google Scholar

*Canadian Geotechnical Journal*, 37(1): 161–170. Google Scholar

*Canadian Geotechnical Journal*, 38: 1254–1265. Google Scholar

*Arctic, Antarctic, and Alpine Research*, 35(1): 82–90. Google Scholar

*Cold Regions Science and Technology*, 13(2): 107–119. Google Scholar

*Cold Regions Science and Technology*, 19(2): 163–175. Google Scholar

*The Avalanche Handbook*. Seattle , Washington: The Mountaineers. Google Scholar

*Annals of Glaciology*, 13: 180–184. Google Scholar

*In*Schumm, S. A. , and Bradley, W. C. (eds.).

*U.S. Contributions to Quaternary Research*. Boulder, Colorado: Geological Society of America, Special Paper 123, 141–165. Google Scholar

*Arctic, Antarctic, and Alpine Research*, 40(1): 148–160. Google Scholar

*Mountain Geomorphology: Geomorphological Process in the Canadian Cordillera*, 215–222. Google Scholar

*Field-Based Statistical Modelling of Snow Avalanche Runout*. M.Sc. thesis, Department of Civil Engineering, University of Calgary, Canada. 128 pp. Google Scholar

*Cold Regions Science and Technology*, 104–105: 23–32. Google Scholar

*Journal of Glaciology*, 43(143): 165–171. Google Scholar

*Journal of Geophysical Research*, 111: http://dx.doi.org/10.1029/2005JF000391. Google Scholar

*In*Proceedings, International Snow Science Workshop, Telluride, Colorado, U.S.A.. Google Scholar

*Cold Regions Science and Technology*, 97: 121–131. Google Scholar

*Natural Hazards, Special Issue on An Assessment of Natural Hazards and Disasters in Canada*, 28: 487–515. Google Scholar

*Earth Surface Processes and Landforms*, 31(11): 1424–1437. Google Scholar

*In*Stoffel, M. , Bollschweiler, M. , Butler, D. R. , and Luckman, B. H. (eds.),

*Tree Rings and Natural Hazards: A State-of-the-Art*. Heidelberg and New York: Springer, 3–23. Google Scholar

## Appendices

### Appendix

#### Sample Calculations

*Example 1. Calculate the return period at a specified runout position, when the return period at the* β *point is known*.

An avalanche path in the Canadian Rocky Mountains starts at 2500 m elevation and runs to 1600 m elevation, covering a total horizontal distance of 2200 m (Fig. A1). The β point is 1500 m from the start zone, and has a return period, *T*β, of 5 years. The following steps show a method for calculating the return period at a horizontal runout distance of 1850 m.

**Step 1:** Calculate probability of not exceeding the desired runout distance (i.e. *X _{T}
* = 1850 m).

Use the cumulative density function of the Gumbel distribution (Eq. 2)

**Step 2:** Calculate the probability of not exceeding the desired runout distance (i.e., *X _{T}
* = 1850 m) given a return period of 5 years at the β point (i.e.,

*T*0 =

*T*

_{β}= 5 years).

Use the cumulative density function of the compound distribution (Eq. 4)

**Step 3:** Calculate the return period at desired runout distance (i.e., *X _{T}
* = 1850 m)

**3a)** Calculate exceedance probability, *E* = 1 - *F*(*v*) = 1 - 0.979 = 0.021

**3b)** Calculate return period, *T*_{1850} = 1 / *E* = 1 / 0.021 = 48.4 ≈ 50 years

*Example 2. Calculate the return period at a desired distance, when the return period at the* β *point is unknown, but the return period at a reference point is known*.

Using the same geometry as Example 1, this example assumes that the return period at the β point is unknown, but that the return period is 1 year for a reference point further upslope at 1320 m (Fig. A2). Changes to calculation method are shown in **bold**.

**Step 1:** Calculate probability of not exceeding the desired runout distance (i.e., *X _{T}
* = 1850 m).

Use the cumulative density function of the Gumbel distribution (Eq. 2)

**Step 2:** Calculate the probability of not exceeding the desired runout distance (i.e., *X _{T}* = 1850 m) given a return period of

**1 year at the reference point (i.e.,**

*T*_{0}= 1 years).Use the cumulative density function of the compound distribution

where*v*= variable which describes spatial and temporal probability

**Step 3:** Calculate the return period at desired runout distance (i.e., *X _{T}* = 1850 m)

**3a)** Calculate exceedance probability, *E* = 1 - *F*(*v*) = 1 - 0.980 = 0.020

**3b)** Calculate return period, *T*_{1850} = 1 / *E* = 1 / 0.020 = 50.502 ≈ 50 years