A method for long-term simulation of snow avalanches is developed, based on coupling statistical interpretation of triggering snowfall process and regional avalanche data. The case study area is the Alta Valtellina region, in the northern Italian Alps. Therein, a 21-year-long series of daily snowfall data from 21 snow stations is used to calibrate a daily point snowfall statistical model. Then, a data set including 68 avalanche events from six historical avalanche sites are used to evaluate regionally valid features of avalanche release probability, geometry, and runout. These findings are then used to set up a model for the occurrence of avalanches. One particular case study site is considered, the Vallecetta mountain, of interest because of the considerable number of avalanche events occurring there. Long-term simulation of daily snowfall is performed, which is then fed into a model of snow avalanche occurrence. Snow avalanche simulations are then carried out, resulting in synthetic statistics of avalanche geometry, volume, and runout for a return period of 300 years. These are compared with regionally observed statistics in the considered area, resulting in acceptable agreement. The proposed model allows long-term simulations of avalanche occurrences for evaluation of snow avalanche volume and runout, usable for ecological and geomorphologic purposes. Integration with an avalanche dynamics model would provide long-term avalanche hazard assessment for land use planning purposes.
Introduction
Influence of snow avalanches on morphology and ecology of mountain areas is paramount. Snow avalanches redistribute snow mass deposited during winter, considerably affecting snowpack distribution on a local scale (e.g. Elder et al., 1991; Bozhinskiy and Losev, 1998), reducing accumulations on steep slopes (e.g. Elder et al., 1998), and impacting mass balance of snow-covered and glacierized areas (e.g. Kuhn, 2003; Wagnon et al., 2007). While airborne avalanches provide little geomorphic effect (e.g. Bozhinskiy and Losev, 1998; Jomelli and Bertran, 2001), flowing avalanches result in soil removal in the release and flowing zone (e.g. King and Brewster, 1978; Heckmann et al., 2005), removal of shrubs and trees (e.g. Muntán et al., 2004; Casteller et al., 2007; Reardon et al., 2008), sediment transport at the avalanche-to-soil interface (e.g. Jomelli and Bertran, 2001; Schrott et al., 2003), formation of impact pool-mound complexes (e.g. Smith et al., 1994), and modified landforms in alpine areas (e.g. Jomelli and Francou, 2000). Avalanches may influence the erosional features of headwater streams of active river channels, cascading into sediment yield from upper catchments (e.g. Ackroyd, 1987). The frequency and magnitude of avalanches influence soil composition and sediment properties (e.g. De Scally and Owens, 2005), sorting in the flowing and deposition area (e.g. Luckman, 1977; Gardner, 1983; Schaerer, 1988; Bell et al., 1990; Blikra and Sæmundson, 1998; Heckmann et al., 2002; Decaulne and Saemundsson, 2006), and also soil chemical composition (e.g. Freppaz et al., 2006). Due to considerable soil removal on their track, avalanches may result in enhanced soil instability (e.g. Heckmann et al., 2002; Freppaz et al., 2006).
Also, the occurrence and size of snow avalanches affect biophysical diversity (e.g. Geertsema and Pojar, 2007), influence forest structure (e.g. Patten and Knight, 1994; Bebi et al., 2001; Kulakowski et al., 2006) and yield of dead vegetation to streams (e.g. Fetherston et al., 1995; Bartelt and Stöckly, 2001; Benda and Sias, 2003), and also interact with the dynamics of forest fires (e.g. Veblen et al., 1994). Furthermore, the structure and composition of avalanche tracks provide unique habitat for animal and plant species (e.g. Erschbamer, 1989; Mace et al., 1996; Krajick, 1998; McLellan and Hovey, 2001). Accordingly, modeling of flowing avalanche occurrence, volume, and runout and the geomorphologic effect therein is of interest.
The snow depth in the avalanche release zone is often assumed to coincide with the snow depth precipitation in the three days before the event, H72 (e.g. Burkard and Salm, 1992). Statistical analysis of H72 is therefore used to evaluate snow depth (and volume) at an avalanche start for given return periods (Hopf, 1998; Barbolini et al., 2002, 2003; Ancey et al., 2004; Bocchiola et al., 2006; 2008a). For instance, the Swiss procedure (Sp; e.g. Salm et al., 1990), adopted with few changes also in Italy (Barbolini et al., 2004), provides mapping criteria for dense snow avalanches with return periods T = 30 and T = 300 years which uses as an input the T-years value of H72 for T = 30 and T = 300, based on the assumption that the avalanche return period matches that of the triggering snowfall. Accordingly, unusually extreme (i.e. greatest annual) avalanches, i.e. those with high return periods, are of interest. However, estimation of snow mass budget and avalanche driven sediment transport might require the description of avalanche dynamics at a smaller time scale.
Long-term simulation based on climatic input, either observed or simulated, is widely adopted, for instance to evaluate flood frequency (e.g. Cowpertwait, 1994, 1995; Rahman et al., 2002; Rulli and Rosso, 2002), soil erosion (e.g. Rulli and Rosso, 2005), and even avalanche occurrence aimed at hazard mapping (e.g. Ancey et al., 2004). Here, a regionally valid method for long-term simulation of snow avalanches is proposed and tested.
The index value approach is adopted, which implies that values of a variable that are scaled, i.e. divided by an index value, have identical frequency distributions across all sites within a given homogeneous area, or region. This creates increased sample dimensionality for distribution fitting. This approach was developed in the field of flood hydrology (e.g. De Michele and Rosso, 2001, 2002; Katz et al., 2002; Bocchiola et al., 2003) and was successfully applied, among others, to evaluate H72 (e.g. Barbolini et al., 2002; Bocchiola et al., 2006, 2008a) and daily Snow Water Equivalent (herein SWE; e.g. Bocchiola and Rosso, 2007a). Here, the index value approach is tentatively adopted to interpret scarce avalanche data from different sites, so increasing sample dimensionality for model estimation.
The case study area is the Alta Valtellina region in northern Italy, where meteorological data as well as avalanche data from six historical avalanche sites are available (see Bocchiola and Medagliani, 2007; Arena Lo Riggio et al., 2009). A 21-year-long series of daily snowfall data is used to calibrate a statistical model able to simulate daily snowfall. This provides the meteorological factor for triggering of avalanche H72 evaluated at the daily scale, H72d. A probabilistic framework, conditioned on the magnitude of H72d, is sketched to trigger avalanche events (e.g. Bozhinskiy and Chernous, 1986; Ancey et al., 2004), based on the findings from the avalanche sites. The geometry of the avalanches, including release depth, width, and length, is interpreted using regional statistical distributions evaluated from the observed data. Then, the avalanche runout needs to be evaluated.
The current approaches to avalanche modeling are mainly aimed at the assessment of runout and flow velocity for hazard mapping (e.g. Salm et al., 1990; Barbolini et al., 2004) and risk assessment (e.g. Keylock and Barbolini, 2001; Fuchs et al., 2004; Fuchs and McAlpin, 2005).
Avalanche runout can be estimated according to two main approaches, i.e. empirical (e.g. McClung et al., 1989; McClung and Mears, 1991) and dynamical (e.g. Barbolini et al., 2000; Salm, 2004). The former include, among others, statistics of extremes based on topographic approaches (e.g. α-β model: Lied and Bakkehøi, 1980; Martinelli, 1986; and the runout ratio: Smith and McClung, 1997; McClung, 2001; Keylock, 2005), and empirical regression against terrain features and avalanche size (e.g. Bovis and Mears, 1976; Gruber and Sardemann, 2003; Maggioni and Gruber, 2003; McClung, 2003). The latter include numerical models, either based on the center of mass hypothesis (e.g. Voellmy, 1955; Perla et al., 1980) or continuous (e.g. Bartelt et al., 1999; Christen et al., 2002; Naaim et al., 2004), also including mass uptake (e.g. Sovilla and Bartelt, 2002; Naaim et al., 2003; Eglit and Demidov, 2005; Sovilla et al., 2006, 2007). Empirical methods often rely on regional calibration (e.g. Lied and Bakkehøi, 1980; McClung and Mears, 1991) based on similar topography, used to increase otherwise scarce sample dimensionality (e.g. Keylock, 2005). Similarly, dynamic models are often calibrated against observed runout (while along-track depth and velocity are seldom available) for a number of avalanche sites inside a climatically and topographically homogeneous region (e.g. Barbolini et al., 2003). Dynamic models usually apply to extreme (i.e. of considerable size, and hence return period) avalanches (e.g. Salm, et al., 1990; Bartelt et al., 1999; Christen et al., 2002), while dynamics of small sized (i.e. more frequent) avalanches, which can contribute significantly to geomorphic action (e.g. Luckman, 1977; Ackroyd, 1987), are generally less investigated. Also, proper validation of dynamic models requires evaluation of along-track depth and velocity, which is not available here (and is generally seldom available). Here, a simple regionally valid rule is found, based on the available avalanche data, to evaluate the runout distance. One particular avalanche site is then considered for model application, the Vallecetta mountain, close to the city of Bormio (Sondrio Province, SO), where a considerable number of historical avalanche data is available and that is of interest because of the considerable avalanche events observed therein (e.g. Riboni et al., 2005; Bocchiola and Rosso, 2007b). Simulation of H72d and occurrence of snow avalanches is performed, resulting in estimates of avalanches with a return period up to T = 300. The statistics of the geometry, volume, and track length of the resulting avalanches are then compared with the statistics of the avalanche runout and volume observed in the considered site.
Study Area and Available Database
HISTORICAL RECORDS OF SNOW DEPTH
The study area (Fig. 1a) covers the mountainous area of the Lombardia region, in the central Alpine and pre-Alpine area. In this area, at least 40 snow measurement stations are present (see Bocchiola et al., 2006; Bocchiola and Rosso, 2007b), mainly adopted for avalanche warning purposes. These are property of the Interregional Association for Snow and Avalanches (AINEVA, 33 stations), located at the Snow and Meteorological Center of Lombardia Region in Bormio city (see Fig. 1a) and of the Regional Agency for the Protection of the Environment (ARPA) of the Lombardia region (7 stations, concentrated in Valmalenco, shaded area in Fig. 1a), located in Sondrio city. In a recent study, the authors used daily snowpack Hs data from this network to demonstrate the regional homogeneity of the frequency distribution (i.e. growth curve) of the greatest annual dimensionless values of H72, namely H72* = H72/μH72, with μH72 local (at a site) average or index value in the considered area (Bocchiola et al., 2006). Notice that the Sp suggests evaluation of H72 as the “increase in snowpack depth during a period of 3 consecutive days of snowfall” (e.g. Sovilla, 2002; Barbolini et al., 2004), thus in practice including changes in snowpack density due to settling. For consistency with the Sp, here the daily Hs data from one station (Bormio, BOR) are used to evaluate the daily snow depth increase Hd as the daily positive variation in Hs. This is different from what was done by Bocchiola and Rosso (2007a), where daily SWE was evaluated by freshly fallen snow depth Hn and its density ρn, the difference being due to snowpack settling. We choose here to model the daily Hd instead of directly modeling H72. This allows us, on one side, to describe the snowfall process at the most intuitive daily scale (i.e. the same as the available data), while, on the other side, this might allow for consideration of snowfall dynamics for avalanche simulation at a different time scale than the three day one (e.g. one day or more), in those cases where the hypothesis of avalanche triggering based on H72 would not be accurate (e.g. in the Maritime Alps of France, where H48 is used; M. Naaim, personal communication, February 2008). BOR station is the most centered with respect to the considered avalanche sites, and also the closest one to the Vallecetta site, which will be investigated in detail. Also, this snow station shows one of the longest available data sets, because there are daily snow depth data available for 21 years (winter 1986–2006). This site is located at a relatively high altitude (1960 m a.s.l.), as compared with other stations in this area (see e.g. Bocchiola et al., 2006). However, reference control is carried out by considering the data from three other stations, namely Cancano (CAN), Santa Caterina Plaghera (SCA), and Livigno San Rocco (SRO), featuring also 21 years of data each.
Bocchiola and Rosso (2007b) showed that in the investigated region, scaling of the index value μH72 with altitude A shows a less homogeneous behavior than the dimensionless quantiles H72* (Bocchiola et al., 2006). This occurs because the value of μH72 interprets the local variability of the processes influencing the average magnitude of the observed values of three days snow depth. This issue is often observed in the field of flood hydrology, when regional homogeneity of the growth curve does not necessarily imply homogeneity of the index flood values (i.e., their relationship with geomorphologic attributes, particularly with drainage area), thus requiring a further step in the regionalization procedure (“hierarchic approach” to regionalization; see Gabriele and Arnell, 1991). Bocchiola and Rosso (2007b), using cluster analysis (CA; see Baeriswyl and Rebetez, 1997), divided the considered region into two sub-regions where homogeneous scaling of μH72 against A is observed, indicated in Figure 1a and referred to as S-W and N-E.
To simulate Hd in the release area of avalanches, usually located at high altitude where no snow gages are available, one can use the dependence of its average value μHd upon A (Bocchiola et al., 2008a). In such cases, evaluation of average snow depth from altitude A is reasonable, because A is possibly the factor that mostly influences the distribution of snowfall in space (see e.g. Barbolini et al., 2002, Bocchiola and Rosso, 2007b; Bocchiola et al., 2008a). Because here different scaling of μH72 with A is observed in the two sub-regions S-W and N-E, it is expected that also μHd behaves accordingly. Therefore, data of Hd from the 21 available stations in region N-E are also included in the analysis, to evaluate scaling of μHd against A.
HISTORICAL RECORDS OF SNOW AVALANCHES
The available avalanche data (Fig. 1b) cover the eastern part of Valtellina Valley. Therein, six historical avalanche sites were investigated for 68 observed avalanche events. The six sites are Val Cerena (herein CER), Val Gandolina (GAN), Val Novalena (NOV), Rastello (RAS), Vallaccia (VCA), and Vallecetta (VCE). Avalanche data were retrieved mainly from the database of the Rangers of Sondrio City (Corpo Forestale dello Stato, CFS) and also from the archives of the Italian Association for Snow and Avalanche (AINEVA). Inside these databases, sparse avalanche events are found dating back to 1886, but more accurately measured in the last 30 years or so (see Bocchiola and Medagliani, 2007). Snow avalanche data included avalanche main cause (i.e. heavy snowfall, rain, and temperature, while anthropic action is not considered as a main cause, but rather as a triggering effect), type of avalanche (e.g. slab, loose, superficial, and bottom), geometry at release (i.e. altitude, fracture depth, release area) (sometimes mapped on a chart at a scale of 1∶25,000). In most cases, the area of release is evaluated by remote observation (e.g. use of binoculars). In situ measurement is carried out only in case of relevant avalanche size. This introduces some degree of uncertainty in the evaluation of the avalanche geometry. Also, deposition data are reported, including morphology (open slope, channelized, gully, etc.), deposition altitude, deposition geometry (also on a chart at a scale of 1∶25,000), and volume. In situ measurements are carried out in the deposition area if the avalanche event involves people (or property), while if not, remote assessment usually is carried out. Unfortunately, no measurements are available of the deposited snow density, of interest in evaluating the avalanche type and properties. The presence of avalanche defense structures, including rakes, dikes, wedges, diversion structures, etc., are recorded where available. A summary of the database is reported in Table 1.
Table 1
Main features of the selected avalanches. CER—Val Cerena, GAN—Val Gandolina, NOV—Val Novalena, RAS—Rastello, VCA—Vallaccia, VCE—Vallecetta. A0 is release altitude, h0 release snow depth, W0 release width, L0 release length, S0 release area, V0 release volume, Ar altitude at deposition, L runout length, Vr volume at deposition. Bold indicates length and volumes at release that required estimation using GIS and avalanche track maps data to be used in the analysis. Italic indicates volumes that are found to be less accurate and were not used in the evaluation of growth index. Cause indicates snowfall (SF), rainfall (R), or temperature (T) due avalanches. Types, first letter indicates either Slab (S) or Loose (L) avalanches, second letter indicates either Surface (S) or Bottom (B) avalanches. Missing values are cases when either a feature was not measured, or it was not reported in the avalanche forms, or was unintelligible.
Methods: Snow Depth Data
DISTRIBUTION OF DRY AND SNOWY PERIODS
To provide synthetic climatic input for avalanche long-term simulation, the available 21 years (1986–2006) of observed series of daily snowfall at BOR station was used to tune a statistical model (see e.g. Ancey et al., 2004). The winter period was investigated, from 1 November to 30 April, because no noticeable snow precipitation was observed before or after this period in the study area. In this preliminary approach, seasonal (e.g. monthly) difference was not considered, and the snowfall series was treated as homogeneous during winter. In the future, more accurate description could be achieved by seasonal partitioning. First, the duration of dry and snowy periods, τd and τs, was investigated (e.g. Kottegoda et al., 2000). Analysis of the database resulted in a total of 556 wet and snowy spells, i.e. of 2366 dry (1949) or snowy days (417). First, (linear) correlation analysis was carried out to evaluate the dependence between dry and snowy period duration, which might be included in some cases for precipitation simulation purposes (e.g. Salvadori and De Michele, 2006). For the BOR station, this resulted in negligible correlation (ρτd,τs = 0.05, not significant for a reference significance level α = 5%). Because no sensible correlation exists, duration of dry and snowy periods can be treated independently for present purposes. Following Kottegoda et al. (2000), a geometric (GE) distribution was tentatively used here for distribution fitting of snowy and dry periods. Use of GE distribution stems from modeling the number of dry and snowy days according to a series of Bernoulli trials, each one with given probability of snowing. In some cases, use of log series (LS) might be recommended for τd (e.g. Kottegoda and Rosso, 1997; Kottegoda et al., 2000). Other distributions are suitable for the length of snowy and dry periods, including Exponential (EXP) and Weibull (WE; see e.g. Veneziano et al., 2002). Here the four mentioned distributions were tested for accommodation of both τd and τs. Concerning τd, the GE distribution gave considerably better fit than the WE and EXP, and slightly better than LS (not shown; see Pagani and Sala, 2007). Concerning τs, all the proposed distributions gave in practice equivalent fitting. For simplicity and consistency with τd, here GE distribution is adopted. The parameters of the mentioned distribution are reported in Table 2. In view of these results, the snowfall process is here modeled as a sequence of snowy and dry periods, with duration independent of each other, and evaluated according to a proper GE distribution. As reported in the section Historical Records of Snow Depth, reference control is carried out by considering the data from three other stations, namely CAN, SCA, and SRO. The results reported here concerning τd and τs for BOR station was found to be in practice equivalent for all four stations. Based on these findings, it is tentatively assumed that the distribution of τd and τs is valid for the considered region and can be therefore used to simulate snowfall process in the considered avalanche sites.
Table 2
Proposed statistical distributions for snowfall and avalanches. Dist. = distribution; see text for definitions.
DISTRIBUTION OF DAILY SNOW DEPTH Hd
First, the dependence between τs and cumulated snow depth Hd was tested. In fact, a direct dependence can exist between the two (e.g. Kottegoda et al., 2000; Salvadori and De Michele, 2006); however, no correlation was found therein (ρHdτs = 0.02, not significant for α = 5%). Also, the independence between duration of the dry periods τd and subsequent snow depth Hd in the first snowy day (Kottegoda et al., 2000) was verified (ρHdτd = −0.04, not significant for α = 5%). Then, the dependence between Hd in two consecutive snowy days was investigated (ρHd,Hd+1 = 0.08, not significant for α = 5%). Eventually, Hd distribution was evaluated from the available database of 417 events and treated as independent according to these findings. Particularly, the distribution of the Hd values, made dimensionless with respect to their local index value, Hd* = Hd/μHd was investigated. Distribution fitting of Hd* is tentatively carried out using GA, LN, WE distribution (e.g. Bocchiola and Rosso, 2007a), and EXP distribution, suitable for precipitation depth (e.g. Bacchi et al., 1994). The best visual fit (see Pagani and Sala, 2007) is given by the GA and WE distribution, while LN and EXP fail to meet the confidence limits according to the Kolmogorov Smirnov goodness of fit test (herein KS; see e.g. Kottegoda and Rosso, 1997) for a level α = 5%. The GA distribution is chosen here (reported in Table 2) because of its frequent use to model rainfall (e.g. Katz, 1999; Wilby and Wigley, 2002) and also daily snow depth Hd (e.g. Skaugen, 1999). Also for Hd*, reference control is carried out by considering the data from three other stations, namely CAN, SCA, and SRO. The results reported here for BOR station was found to be in practice equivalent for all the four mentioned stations. It is therefore assumed that the distribution of Hd* is valid for the considered region and can be used to simulate snowfall process in the considered avalanche sites. However, to fully simulate snowfall depth Hd in an unmeasured site, e.g. in the release area of an avalanche, one needs to know the index value μHd. Using the daily snow depth database from the 21 snow stations in the N-E region reported in Figure 1a, scaling of μHd against altitude A was investigated (Fig. 2). Using weighted (i.e. on their estimation variance, depending upon the number of available data) linear regression, a significant dependence of μHd upon A was found, as reported in Table 3, in spite of the somewhat low determination coefficient, R2 = 0.65. This method is used to evaluate the index value μHd in ungauged avalanche release areas.
Table 3
Proposed regression equations.
DISTRIBUTION OF H72
The distribution of H72d is then evaluated by the sum of Hd on a three day moving window. Notice that because Hd is distributed according to a GA distribution, its three day sum is also distributed as a GA distribution (see e.g. Skaugen, 1999). Bocchiola et al. (2006) investigated the regional distribution of the greatest annual H72 for a wider region, including the one here considered, finding in practice a Gumbel (GU) distribution. Also Bocchiola et al. (2008a) provided regional frequency curves for H72 for Switzerland, finding a satisfying fitting of the Gumbel distribution for six homogeneous regions out of seven.
This is consistent with the findings here, because the greatest (i.e. extreme) values of a variable displaying GA distribution display in turn a GU distribution (e.g. Kottegoda and Rosso, 1997). Several tests (see Pagani and Sala, 2007) indicated consistent statistics (i.e. statistically equivalent variance and average) of the synthetically simulated series of both Hd and H72d with the observed ones.
Figure 3 shows the comparison, of interest for the purpose of extreme avalanche evaluation, of the greatest annual values of H72d, H72dmax in BOR station, against those estimated according to the regional growth curve (Bocchiola et al., 2006), used here because it is more reliable than the single site one for the considered return periods. Also, the locally observed values are reported for comparison. Synthetic simulation of H72d for a period of 195 years is carried out using the Monte Carlo method (see Kottegoda and Rosso, 1997), resulting in an estimated return period of 300 years (according to the plotting position APL, F = (ord − 0.65)/ntot, where ord is the ranking position of the particular value of H72d in the sample, ordered in decreasing order of magnitude and ntot = 195). To avoid sample effects, particularly for high return periods, 100 repetitions were carried out, for a total of 1.95E4 simulated values of H72dmax. Then, the average (i.e. average of 100 values for each return period) plotting position so obtained was put in Figure 3. Apart from some slight underestimation, the synthetic series seems to fit well the regional one, and is well contained inside the confidence limits (α = 5%; see Bocchiola et al., 2006, 2008a, about their evaluation), indicating that the proposed model represents well the snowfall process, also for considerable return periods.
Methods: Avalanche Data
ASSESSMENT OF AVALANCHE DATA
For those avalanche events that were reasonably well mapped, the main avalanche features were evaluated. These are avalanche release altitude A0 and geometry (greatest width W0, surface S0, snow depth h0 at release, and in some cases maps at 1∶25,000), and runout altitude Ar, length L, and volume Vr, and in some cases maps at 1∶25,000. Frequently, some of the data were missing and therefore several avalanches could not be fully analyzed. A screening was carried out for the presence of manmade structures (walls, rakes, dikes, wedges, diversion structures, etc.) that might have influenced avalanche properties (e.g. release shape or volume, or runout distance, etc.). However, for the studied cases, no considerable influence was observed.
Also, a preliminary screening was carried out to highlight features that are clearly inaccurately estimated. These include erroneous release areas or depth, evaluated e.g. by consistency with released and deposited volume, or erroneous estimates of snow volumes at deposition, evaluated by consistency with snow depth and deposition geometry. Also, GIS-based topographic analysis was carried out to check accuracy of the release and deposition areas and of the indicated topography in the deposition area. This was based on a 20 m cell size DEM of the considered area. The length of the release area L0 was also calculated, based on the available maps, also verified through topographic analysis. In Table 1 the main avalanche properties are reported, together with a judgment of the quality of the different features. The assessment of the database as here carried out surely implies a degree of subjectivity, but it is deemed to avoid erroneous interpretation of the avalanche process due instead to evident mistakes in data retrieval.
AVALANCHE TRACK GEOMETRY
Figures 4a–4c report the track of the avalanches, including main avalanche profile, track width, and the release and deposition altitude for the avalanche events here mapped. The given profile is the most indicative one, as defined after a preliminary analysis based on dynamic simulation of the greatest (i.e. with the lowest runout altitude) observed avalanche (historical end mark) using a simple dynamic model (Voellmy Salm model, VS; see e.g. Bocchiola and Medagliani, 2007). The six sites show a similar morphology, featuring a wider release area, followed by a channelized flowing zone, and a depositional fan at the valley bottom. The average slope of the release, channelized (i.e. flowing), and deposition, or runout areas, are reported in Table 4. Also, the percentage length of each zone as compared to the whole track length (from the highest altitude to the lowest deposition altitude) is given. Release zone is defined coupling slope analysis (i.e. 30°–60°; e.g. Maggioni and Gruber, 2003) and observed release area. We labeled as flowing area the channelized zone, until start of the alluvial fan. The runout zone is defined between the end of the flowing zone and the historical end mark (i.e. longest runout). The six avalanche sites display reasonably similar values of the slope in the three areas so defined, and also a similar percentage length, as indicated by the corresponding low coefficient of variation CV[.]. This seems to suggest that a similar behavior might be observed in terms of avalanche geometry and runout, provided the proper scale factors are considered (e.g. for different altitude and track length), and therefore index value approach might be tentatively adopted.
Table 4
Avalanche track geometry. E[s0] is average track slope including the release area, E[sf] is average track slope including the flow area, E[sf] is average track slope including the runout area (average track slope). L0′ is percentage of track length including release area, Lf′ is percentage of track length including flow area. For other abbreviations, see Table 1.
PROBABILITY OF RELEASE
First, probability of avalanche release was investigated, as compared to snow depth at release. Preliminary analysis (reported in Bocchiola and Medagliani, 2007) showed that for 68% of the observed avalanches the triggering factor was linked to heavy precipitation in the 72 hours preceding the event (similar to e.g. Decaulne and Saemundsson, 2006). Therefore, it seems reasonable to investigate the probability of release as connected to the three day snowfall depth, H72d. It is therefore assumed for simplicity that the release depth coincides with H72d at the release site, h0 ≈ H72d. Then, the regionally valid distribution of H72 evaluated by Bocchiola et al. (2006) was used to estimate the frequency of non-exceedance of h0, Fh0,H72d. This was done by taking the dimensionless values of h0, as h0* = h0/μH72, where μH72 is the local average, or index value of H72. This is possible because the variable h0* so obtained can be sketched as distributed according to a regionally valid distribution, as shown in Bocchiola et al. (2006). Then, the frequency of exceedance of the observed release depth h0, Fh0 has been evaluated using empirical plotting position. So doing, the ratio F0 = Fh0,H72d/Fh0 gives an estimate of the probability that an avalanche occurs due to a given snowfall depth. The values of F0 are reported in Figure 5. In spite of some sparseness, also due to possible uncertainty in the remote evaluation of snow depth at release, the observed distribution is well fit using a uniform distribution, UN. This is shown together with confidence limits for α = 5%, according to KS test. Notice from Table 1 that here no avalanche event was observed for h0 < 0.2 m. Therefore, if H72 < 0.2, the probability of an avalanche should be set to zero. The proposed distribution F0 is reported in Table 2 and will be used here to trigger avalanche release conditioned on snowfall H72d for simulation purposes.
RELEASE ALTITUDE
The observed avalanche altitude (i.e. altitude of snow fracture, or avalanche crown line) at release is variable, as shown in Figure 4. In each site, greatest (±) deviations are observed in the order of 10 to 15% or so with respect to the average value. The assumption is made here that while the average altitude at release is linked to local topography, a random component exists that results in changing of the release altitude for each single event. If one tentatively assumes that such a random component is similar in each avalanche site, a regional distribution of the dimensionless release altitude A0* = A0/μA0, with μA0 average, or index local value, can be estimated by grouping the whole observed sample. The calculated sample frequency is shown in Figure 6. It can be seen that this has a relatively homogeneous behavior at the considered sites. The sample frequency is reasonably well accommodated using a Normal (henceforth NR) distribution, reported in Table 2.
AVALANCHE SHAPE AT RELEASE
Here, the shape of the avalanche release zone is investigated. The observed values of release width W0 and length L0 were adopted in the hypothesis of simplified avalanche geometry (i.e. release area featuring a rectangular shape with length L0 and width W0). Then, the statistical distribution of W0 and L0 was investigated. First, the existence of a dependence between h0 and either W0 or L0 was investigated, because it might be that a link exists between the snow fracture depth and release area shape. However, no significant dependence was found between these variables, neither by visual assessment nor after correlation analysis (correlation coefficient on the whole sample ρW0,h0 = 0.13 and ρL0,h0 = −0.09, not significant at a reference α = 5% confidence level). Therefore, avalanche width and length at release seem not linked to snow fracture depth, at least for these sites. Mutual dependence of L0 and W0 is also preliminarily investigated, showing in practice no significant correlation (ρL0,W0 = 0.16, not significant for α = 5%). It is also made here the assumption that the average conditions at release are linked to local topography, while a random component exists, which is regionally valid. The sample frequency is then calculated of the dimensionless values, W0* = W0/μW0 and L0* = L0/μL0. These are reported in Figures 7 and 8. The frequencies of both W0* and L0* seem quite homogeneous in the considered sites and reasonably well accommodated by Gamma distributions (GA; the Lognormal distribution, LN, seems less fit), reported together with uncertainty bounds for confidence level α = 5%, according to the KS test. The parameters of the GA for W0* and L0* are reported in Table 2.
AVALANCHE RUNOUT VS. RELEASE VOLUME
Avalanche runout is linked to snow mass and volume (e.g. Hutter, 1996; Maggioni and Gruber, 2003) and avalanche dynamics models are sensitive to release mass (e.g. Voellmy, 1955; Bartelt et al., 1999), and greater avalanches run farther. However, dynamic models require proper parameter tuning, which is usually carried out for extreme (i.e. with high return periods, 30 to 300 or so in Switzerland; e.g. Salm et al., 1990) events, which are more likely to be mapped, while accurate information about more frequent (i.e. with low return periods) events is seldom available, thus making use of dynamic models less accurate. Further, dynamic model estimation and use for long-term simulation requires some probabilistic assessment of model parameters (e.g. Ancey et al., 2003; Naaim et al., 2004), which goes beyond the scope of the present paper. Here, the relation between snow mass at release and runout is directly investigated. In view of the different geometry and size of the considered avalanche sites, again dimensionless values are used with respect to the average, site specific, values, namely V0* = V0/μV0 and L* = L/μL. The hypothesis is made that a regionally homogeneous mechanism is present that relates avalanche volume and track length. This is similar to what is done in empirical methods based on topography (e.g. Lied and Bakkehøi, 1980; McClung et al., 1989; McClung and Mears, 1991), where regional analysis is used to increase sample dimensionality for model estimation (see e.g. Keylock, 2005). In Figure 9, it is shown how a linear regression fits well to the logarithms of the observed values of L* against those of V0* for the considered sites. In Table 3 the calculated (power) law coefficients and their significance is reported. As expected, avalanches with greater volume tend to travel farther, while smaller ones have shorter track length. Albeit the explained variance is not very high (R2 = 0.74), the proposed equation has the merit of providing a rapid evaluation of the expected runout distance of the avalanche. These proposed equations will be used for the assessment of avalanche runout here. Notice, however, that if a proper model parameterization is available, evaluation of the avalanche runout, depth, and velocity (i.e. pressure) can be obtained feeding the topographic and geometric properties of the avalanches to a dynamic model (e.g. Bocchiola et al., 2008b).
AVALANCHE VOLUME AT RUNOUT
It is of interest to estimate the avalanche volume at runout, Vr. Usually, the latter is greater than the volume at release, V0, due to snow entrainment (e.g. Sovilla and Bartelt, 2002; Sovilla et al., 2006, 2007). A growth index Ig can be estimated (e.g. Sovilla, 2004) giving the ratio between the snow mass at runout to that at release. This requires evaluation of both snow density and volume. Here, because no snow density measurements were taken in the post-event field survey, the simple hypothesis is adopted that snow density remains unchanged and the growth index is defined as Ig = Vr/V0. Ig is here estimated for a number of 27 events showing acceptable evaluation of both Vr and V0, according to the preliminary analysis reported above. The average observed value on the whole sample is of E[Ig] = 3.5, with a considerable variability, CV[Ig] = 1.32, ranging from the lowest value of 0.15 (CER15; see Table 1) to the greatest 15.5 (RAS7; see Table 1) (compare e.g. with Sovilla, 2004, giving an average value of Ig = 4.6 for a number of avalanches in Switzerland, with a considerable variability from about 1 to 12). Sovilla (2004) found a substantial independence between Ig and avalanche size and terrain characteristics along avalanche track. Ig is a complex function of snow depth and area at release, along-track erodible snow cover, and snow properties (e.g. Sovilla et al., 2006). Erodible snow cover is in turn linked to deposited snow cover before avalanche release, i.e. H72 along avalanche track and to snowline altitude (e.g. Bianchi Janetti and Gorni, 2007; Bianchi Janetti et al., 2008), and therefore depends on altitude (Bocchiola et al., 2006; Bocchiola and Rosso, 2007b). However, evaluation of eroded area and, therefore, snow mass is not straightforward and requires use of properly developed dynamic models (e.g. Sovilla et al., 2006; Bianchi Janetti et al., 2008).
Here, a number of cases (i.e. 8 cases; see Table 1) were found with Ig ≪ 1, indicating along-track deposition. Apart from mistakes in volume evaluation that might still be present in spite of the preliminary assessment here carried out and explained in section 4.1, in some cases (i.e. CER7, CER12, CER15, NOV5), low values of Ig might be due to the occurrence of an avalanche in late spring (i.e. after mid April), when no considerable snow cover is available for along-track entrainment at low altitudes, and instead deposition might occur due to ground roughness. In other cases (i.e. CER13, CER14, NOV6, NOV9), avalanches occur during late autumn to winter (i.e. November to February), but still along-track deposition might occur.
A preliminary analysis showed that Ig is in practice uncorrelated to snow depth, avalanche width, and length at release, as well as to avalanche runout length. Plotting of Ig against release volume V0 is reported in Figure 10. This can be roughly split into two parts. For low values of V0 (i.e. V0 < 105 or so), Ig ranges from about 1 to more than 15, with no apparent link to V0. For V0 > 105 or so, in practice Ig = 1. The five avalanche events when V0 > 105 occurred in early to late spring (see Table 1; NOV4, NOV8, VCE10, VCE11, and RAS8 [the date of RAS8 is not known exactly, but it happened in late spring]), when little snow cover was available at low altitudes, thus resulting in little snow uptake, if any. Here, it is assumed for simplicity that for V0 < 105, Ig is in practice a random variable, in view of its clear independence from V0, while for V0 ≥ 105, Ig = 1. Albeit this approach is very simple in respect to the complexity entailed in understanding and modeling avalanche mass uptake, it seems that it might be used to capture the statistical properties of the phenomenon and to evaluate in a first approximation the statistics of avalanche volumes at deposition. It was observed here that for V0 < 105, the dimensionless growth index Ig* = Ig/μIg, with μIg average local value, can be fairly well accommodated according to an exponential (EXP) statistical distribution, as shown in Figure 11. Therein, it is seen that a pretty homogeneous distribution is obtained by grouping the single site values of Ig*. This indicates that a regionally similar random component seems to exist that results in variability of the growth index for each event. Therefore, we introduce here random value of Ig* as from the proposed EXP distribution, reported in Table 3.
Long-term Simulation of Avalanche Series
SIMULATION STRATEGY AND APPLICATION TO VALLECETTA AVALANCHE SITE
Avalanche dynamics simulation is carried out daily. First, to simulate snowfall depth at the release altitude in the Vallecetta site, the average value of Hd needs to be scaled for the release altitude A0, as from the proposed equation in Table 3.
For the Vallecetta site, use the average value of the release altitude, E[A0] = 2697 results in μHd = 12.1 cm. Using this value, every day the increase in snowpack depth Hd is simulated according to the model proposed in the section Methods: Snow Depth Data and a simulated series of H72d is obtained. Based on H72d, the probability pav that an avalanche occurs is evaluated using the UN distribution in Figure 5. Extraction of a random value from a UN distribution is used to evaluate occurrence of an avalanche according to a binomial distribution (yes/no with parameter pav). If an avalanche occurs at day d, its release depth is set to h0 = H72d. No snow drift load is applied in the avalanche release area, as suggested by the Sp. Analysis of the avalanche database, including qualitative meteorological information, indicated presence of strong wind in the 3 days before the event only on one date (1 November 1986) out of 21. So, one may assume here that wind upload is not preponderant for avalanche release. Further, neither estimation of wind snow load is available here, nor do we know of any accurate formula for its estimation for this area. However, because we estimate release probability using observed snow depth at release, which may include wind load, we in some way account for such phenomena. Notice that in those cases where snow load can be estimated, it can be used in the procedure. In case of an avalanche at day d, evaluation of H72 in the next day d + 1, H72d+1 accounts for avalanche occurrence (i.e. in practice, H72d+1 = Hd+1). The second step is the random simulation of the release (crown) altitude A0*, according to the NR distribution proposed in the section Release Altitude, and A0 = μA0ṡA0*, with the further constraint that A0 is included into the release zone as deduced in the section Avalanche Track Geometry. Then, dimensionless release width and depth W0*, L0* are randomly extracted from the Gamma distributions reported in Table 3, and W0 = μW0ṡW0*, L0 = μL0ṡL0*, with the boundary conditions given by topography (i.e. greatest width and length). Avalanche release volume is then estimated as V0 = h0′ ṡL0 ṡ W0, where h0′ is the release depth as corrected for local slope following the Sp procedure (e.g. Salm et al., 1990). Then, dimensionless runout L* can be evaluated as a function of dimensionless release volume V0* as explained in the section Avalanche Runout vs. Release Volume, and L = L*·μL. Then, runout volume Vr is evaluated according to V0. If V0 < 105, Ig* is extracted from the EXP distribution reported in the section Avalanche Volume at Runout, and Vr = V0·μIg ṡIg. If V0 ≥ 105, Vr = V0. A sample year of simulation is reported in Figure 12. The sequence of H72d is reported, together with information of avalanche occurrence, volume at release V0, and runout. We report here the absolute runout position measured from the start of the avalanche track, R, rather than the track length, L, which is instead relative to the starting point, which changes from case to case. In this particular year, seven avalanche events were synthetically simulated. Here, simulation was carried out for a number of years ntot = 195, consistent with the snowfall simulation in the section Methods: Snow Depth Data. Again, to decrease sample effects, a number of 100 simulations were carried out, for a total of 1.95·104 years.
RESULTS
In Table 5, some explanatory statistics are reported concerning the avalanche release volume, V0, absolute runout, R, and deposition volume, Vr. The full sample average is reported (including the whole set of simulated avalanches, e.g. E[V0]), together with the sample average of the extreme avalanches (i.e. the greatest annual values E[V0max]) and the greatest observed value, of interest for land use planning (e.g. MAX[V0]). These are compared with the same statistics as obtained from the available sample of data. Notice that in view of the relatively small sample size, the confidence limits σ ± for the average values are considerably wide. Also, notice that for the full sample, average comparison against the observed values is particularly difficult, because the observed statistics are considerably skewed towards the greatest values (i.e. those more likely to be observed), thus resulting in likely overestimation. In view of the small available sample, higher order statistics (e.g. coefficient of variation) were not considered for the comparison.
Table 5
Vallecetta mountain. Some observed and simulated statistics of avalanche data. Symbol σ ± indicates sigma bounds (i.e. mean ± standard deviation). Notice the relative width of the sigma bounds, in view of the small sample size. Bold indicates values falling outside the sigma bounds.
Also considering the relevant degree of uncertainty of the proposed statistics, one notices that the proposed model for long-term simulation results in estimated avalanche features that are compatible with those observed at the investigated sites.
Average yearly number of avalanches was nav = 3.1. Comparison with the observed database is probably senseless, because several (minor) avalanches are likely to be neglected. As a rough comparison, Naaim et al. (2004), referring to a “well known avalanche path” reported of a number of observed avalanches of about 153 in 72 years. On average, here no avalanche occurrence was found in 4.2% of the years (i.e. the average number of years with no avalanche is n0 = 8.2 out of 195), giving the probability that in a given year no avalanche is observed at the Vallecetta site. In Figures 13 to Figure 1415 there are reported the synthetically obtained plotting positions of the greatest annual values of the runout, release, and deposition volume, Rmax*, V0max*, Vrmax*, made dimensionless with respect to their average values in Table 5. These are compared with the dimensionless sample values from the six investigated sites, used to increase sample dimensionality and to allow more likely estimation of the frequency of occurrence for the considered events. To clarify this issue, one could consider e.g. the VCE11 event (Rmax = 4930, Rmax* = 1.56), which occurred in 1886. According to the empirical plotting position from the single site analysis (i.e. using the available local data, for a number of 15 years) its estimated return period (using APL plotting position, see the section Distribution of H72) would be TVCE11,site = 23 years. One can attain unbiased estimation of the related return period according to the historical flood (here, avalanches) approach, giving Thist = (n + 1)/(n′ + 1), with n years of observation and n′ number of years when a given value is exceeded (e.g. Kottegoda and Rosso, 1997; Bocchiola et al., 2003). Here, RmaxVCE11* was attained once over the last 123 (from 1886 to 2007) years, therefore giving TVCE11,hist = 124/2 = 62 years, with sigma bounds (see Bocchiola et al., 2003, Equation 11) T ± VCE11,hist = 35–198 years. Using the regional plotting position of the dimensionless observed runout values, one obtains TVCE11,reg = 75 years (i.e. yT = 4.32; see Fig. 14), closer to the unbiased value as reported than the return period for single site plotting position. The synthetically obtained plotting position of V0max* (Fig. 13) seems to fit well the observed one. Some discordancy is observed for T < 2 years or so. However, this could be due to improper labeling of small avalanches as greatest annual ones, in turn possibly resulting from either lack of observation of small avalanche phenomena, or poor assessment of the release volumes for such small avalanches, as reported. Runout Rmax* in Figure 14 seems quite well represented for the whole range of considered return periods. Eventually, the deposition volumes (Fig. 15) seem well estimated, with a slight discordance for T < 2 or so, again possibly related to poor assessment of small avalanches properties.
DISCUSSION
The acceptable fitting of the simulated avalanche properties to the observed values as reported seems to indirectly confirm the capacity of the proposed model to reasonably represent an avalanche regime in a given site. In spite of the small sample size, suggesting care in interpretation of the results, the proposed approach seems valuable for its application in the field of geomorphology and ecology of mountain ranges, where snow avalanche magnitude, runout, and frequency are acknowledged to play a predominant role.
The present model provides a description of long-term avalanche occurrence, to be used as an input for geomorphologic theories. If a simple, data-driven parameterization would be introduced to link, say avalanche properties with erosional features (as in e.g. Bell et al., 1990), the present approach would provide an estimation of these properties also for low frequency of occurrence, thus allowing long-term evaluation of morphological effects of snow mass movement.
The present approach does not explicitly include use of a snow dynamics model, the parameterization of which goes somewhat beyond the scope of the present paper, while runout and snow uptake are more simply assessed.
However, in case a dynamics model would be available, including i.e. sediment transport from snow avalanches, or vegetation removal (as e.g. in Bartelt and Stöckly, 2001), still the input given by the model would be valuable for long-term simulation. Because H72 can be used as a condition for erodible snow cover for mass uptake simulation (Bianchi Janetti and Gorni, 2007; Bianchi Janetti et al., 2008), the proposed model may give valuable input for long-term avalanche dynamics including mass uptake (e.g. Naaim et al., 2003; Eglit and Demidov, 2005; Barbolini et al., 2005; Sovilla et al., 2006, 2007). Furthermore, the present approach may be modified to include an assessment of medium to long-term scenarios of changes in snowfall dynamics on the avalanche frequencies. Because changes in avalanche regime affect Alpine ecosystems in many ways (e.g. Erschbamer, 1989; Mace et al., 1996; Krajick, 1998; Kulakowski et al., 2006), long-term avalanche scenarios are of interest.
The proposed model accounts for snowfall-triggered avalanche events, usually resulting in dry, or moderately wet flowing avalanches. Albeit these seem to cover the vast majority of the observed avalanche events here (see also e.g. Decaulne and Saemundsson, 2006), and are in widespread use in practice to model avalanche occurrence for long-term simulation, model tuning, and hazard evaluation (Salm et al., 1990; Burkard and Salm, 1992; Bartelt et al., 1999; Ancey et al., 2004; Barbolini et al., 2004), different triggering mechanisms can be observed. These include more frequent wet avalanches, or less frequent (T ≈ 40 years or so) slush avalanches, including mud flows, both usually occurring at spring thaw (May to June, see e.g. Jomelli and Bertran, 2001). Note that the approach here shown accounts in practice for both dry and wet avalanches, without an explicit distinction. This may be introduced e.g. by considering snowpack state modeling in time (e.g. for snow depth and density), to be included in the model. As reported, for the considered avalanches no measurements of snow density were available, which is also of interest to evaluate avalanche type (e.g. Freppaz et al., 2006). For the considered area, the density of freshly fallen snow, Hn, is in practice a LN distributed random variable, with mean slightly decreasing with altitude (e.g. Bocchiola and Rosso, 2007a; Medagliani et al., 2007). Its observed average is ρn,av = 123 kg m−3 (and ρn,BOR = 103 kg m−3; Bocchiola and Rosso, 2007a). The density ρs of the snowpack Hs notably changes with the season (i.e. with day in the winter season d; see e.g. Elder et al., 1991). For instance, Martinelli et al. (2004) used data from a number of 196 snow pits in the Lombardia region, mainly concentrated in the BOR area and ranging from November to May for the period 1997–2002, and found ρs = 161–618 kg m−3, depending linearly on d and averaging to ρs,av = 324 kg m−3. Notice that the Sp suggests use of ρs = 300 kg m−3 for avalanche pressure calculation (Salm et al., 1990; Barbolini et al., 2004), thus considering in practice moderately wet avalanches. Because here the daily variation of Hs, i.e. Hd, is used according to the Sp, use of ρs might be suggested. Also, one could explicitly evaluate snow settling during the three days before the avalanche events to calculate the density of H72 (as e.g. in Martinec and Rango, 1991). Here, the reasonable matching of the statistics of the simulated values of H72d with those observed in the BOR snow series seems to indicate that no considerable settling effect needs to be accounted for when using snowpack series for H72 simulation.
Notice further that the proposed snowfall simulation model with slight changes is suitable for more complex approaches. For instance, it could be used to provide an input for simulation of snowpack dynamics using a more refined model (e.g. SNOWPACK; Bartelt and Lehning, 2002), say by including long-term simulation of fresh snow depth and density based on suitable distributions (e.g. Bocchiola and Rosso, 2007a).
Conclusions
Here, a relatively simple framework based on regional analysis is drawn that allows the building of statistically consistent series of daily avalanche occurrence and magnitude in a given site. Because the proposed approach combines statistical interpretation of climatic forcing with regional analysis of observed avalanche tracks, more accurate description is achieved of the avalanche occurrence process with respect to the presently adopted methods solely based either on runout analysis of extreme (i.e. with high return periods) events, or on derivation of extreme avalanches from extreme snowfalls. In turn, more comprehensive data are required, including avalanche volumes and geometry, as well as snowfall records. However, provided a reasonable regional homogeneity and topographic similarity holds, data from avalanches occurring at different sites inside a given region can be used to increase sample dimensionality. Further improvement could include explicit evaluation of different types of avalanches, and the introduction of different snow density, according to the available literature and the season. Also, further developments will include strategies for the evaluation of the index values of the observed avalanche features starting from some a priori measurable indexes, including e.g. local topography. This would allow use of the proposed approach also in those places where no avalanche data are available to evaluate the required index values. Coupling with an avalanche dynamics model could provide a more physically based description of long-term avalanche dynamics for avalanche hazard assessment. In the future, we will try to extend the study area, so hopefully increasing the avalanche database, to reduce uncertainty in our findings.
Acknowledgments
The authors kindly acknowledge Eng. Alice Pagani and Eng. Giovanni Sala for their contribution to the present research, developed during their MS thesis. The authors kindly acknowledge the personnel of AINEVA, in Bormio, and the Rangers of Sondrio City (Corpo Forestale dello Stato) for providing the snow and avalanche data from their archives and for helpful suggestions. ARPA Lombardia is acknowledged for providing snow data. Funding for the research presented in the present paper was granted by the European Community, through the EU projects AWARE (EC contract 012257) and IRASMOS (EC Contract 018412) and by the CARIPLO foundation of Italy through the project CARIPANDA.