Natural Hazards and Earth System Sciences Regional evaluation of three day snow depth for avalanche hazard mapping in Switzerland

The distribution of the maximum annual three day snow fall depthH72, used for avalanche hazard mapping according to the Swiss procedure ( Sp), is investigated for a network of 124 stations in the Alpine part of Switzerland, using a data set dating back to 1931. Stationarity in time is investigated, showing in practice no significant trend for the considered period. Building on previous studies about climatology of Switzerland and using an iterative approach based on statistical tests for regional homogeneity and scaling of H72 with altitude, seven homogenous regions are identified. A regional approach based on the index value is then developed to estimate the T -years return period quantiles of H72 at each single sitei, H72i(T ). The index value is the single site sample average μH72i . The dimensionless values of H ∗ 72i=H72i /μH72i are grouped in one sample for each region and their frequency of occurrence is accommodated by a General Extreme Value, GEV, probability distribution, including Gumbel. The proposed distributions, valid in each site of the homogeneous regions, can be used to assess the T -years return period quantiles of H ∗ 72i . It is shown that the value ofH72i(T ) estimated with the regional approach is more accurate than that calculated by single site distribution fitting, particularly for high return periods. A sampling strategy based on accuracy is also suggested to estimate the single site index value, i.e. the sample average μH72i , critical for the evaluation of the distribution of H72i . The proposed regional approach is valuable because it gives more accurate snow depth input to dynamics models than the present procedure based on single site analysis, so decreasing uncertainty in hazard mapping procedure. Correspondence to: D. Bocchiola (daniele.bocchiola@polimi.it)


Introduction
Land use planning in mountain ranges requires reliable avalanche hazard assessments (e.g.Salm et al., 1990;Barbolini et al., 2004) and management of the risk therein (e.g.Bründl et al., 2004;Fuchs et al., 2004;Fuchs and McAlpin, 2005;Keiler et al., 2006).The current approaches to avalanche hazard mapping couple avalanche dynamics modelling (Salm et al., 1990;Harbitz, 1998;Bartelt et al., 1999;Ancey et al., 2004) with statistical analysis of snow depth at the avalanche starting zone (Burkard and Salm, 1992;Hopf, 1998;Barbolini et al., 2002).In avalanche hazard mapping exercise, the T -years return period avalanches at a given site are computed and their run out zone and pressure are evaluated (e.g.Salm et al., 1990).The snow depth in the avalanche release zone is often taken as the snowfall during the three days before the event, H 72 , evaluated with respect to a flat area and then properly modified for local slope conditions and snow drift overloads (Salm et al., 1990;Burkard and Salm, 1992).The Swiss procedure (hereon Sp, e.g.Salm et al., 1990), also used as a reference in other countries (e.g. in Italy, Barbolini et al., 2004), provides mapping criteria for dense snow avalanches, based on the definition of a more dangerous red zone (impact pressure T ≤300 years and Pr≥30 kPa, or T ≤30 years and Pr≤30 kPa), and of a less dangerous blue zone (impact pressure 30≤T ≤300 years and Pr≤30 kPa).Therefore, avalanche hazard mapping based on these criteria requires as input for each avalanche site the H 72 T -years value, for T =30 and T =300, at least.
Estimation of H 72 T -years quantiles is often carried out by distribution fitting of the single site observed data, i.e. of the maximum annual observed values of H 72 for a range of observation years.According to the theory of extreme values (e.g.Darlymple, 1960;Hosking et al., 1985;Hosking and Wallis, 1993;Kottegoda and Rosso, 1997;De Michele and Rosso, 2001;Katz et al., 2002) to provide reliable estimates Published by Copernicus Publications on behalf of the European Geosciences Union.Gumbel variable y T =− ln(− ln((T −1)/T )) µ H 72i Single site (i) sample average (index value) of Single site (i) scaled mean square error of the single site sample average estimate Single site (i) standard deviation of the H 72 single site sample average σ H * 72i T Single site (i) standard deviation of the H * 72 T -years quantile σ *

T i
Single site (i) dimensionless standard deviation of the H 72 T -years quantile of the T -years quantiles using empirical distribution fitting, a least number of observations is required, in the order of n obs =T /2.For T =300 years, this amounts to 150 years of sampled data.In Switzerland, except for few cases, shorter series of observed snow depth are available, covering on average a period of 50 years or so (e.g.Laternser, 2002).The lack of observed data for distribution fitting of extreme values in hydrological sciences can be overcome by using regional approaches, including index value approach (Darlymple, 1960;Sveinsson et al., 2001).This implies that values of a hydrological variable that are scaled, i.e. divided by an index value (e.g.Bocchiola et al., 2003), have identical frequency distributions across all sites within a given homogenous area, or region (Hosking and Wallis, 1993;Burn, 1997).The regional distribution is usually named growth curve and is more accurate as compared to single site analysis (e.g.Madsen et al., 1997;De Michele and Rosso, 2001;Katz et al., 2002).The regional approach requires a preliminary analysis to verify the regional homogeneity of the considered process and its spatial extent (e.g.Hosking and Wallis, 1993;Burn, 1997;Castellarin et al., 2001;De Michele and Rosso, 2002;Bocchiola et al., 2006;Bocchiola and Rosso, 2007).
Here, the authors try to apply the regional approach to the mountainous area of Switzerland, covering about 2.4 E 4 km 2 , i.e. 60% of the country.Former studies have investigated the degree of homogeneity of Switzerland with respect to the distribution of rainfall and snow precipitation in gauged sites (e.g.Rohrer et al., 1994;Baeriswyl and Rebetez, 1996;Laternser, 2002).Based on these studies, the authors test here the degree of homogeneity of Switzerland concerning H 72 using snow depth measurements for a network of 124 snow measuring stations.First, a test for stationarity in time of the series of H 72 is applied for each snow station with at least 30 years of available data.Then, considering an increasing number of regions the degree of homogeneity of the latter is assessed using statistical tests based on observed values of H 72 , scaled to their index value, or sample average µ H 72 .Also, the relationship between µ H 72 and altitude A is considered to define homogeneity, because µ H 72 often depends on altitude inside homogeneous regions (e.g.Barbolini et al., 2003;Bocchiola et al., 2006).
After the evaluation of the regions, distribution fitting is carried out by grouping the single site values of H 72 in each region, scaled with respect to µ H 72 , namely H * 72 =H 72 /µ H 72 .A General Extreme Value (GEV) distribution, including EV1 (i.e.Gumbel) is adopted, suitable for extreme value analysis.Because µ H 72 can be estimated from the observed data at a single site, its accuracy depends on the number of observed years at the specific site.A framework is drawn here to quantify the accuracy of the T -years quantile estimates depending on the number of the available recorded years, so providing an accuracy driven sampling strategy.Also, regression equations against altitude are given were suitable, to indirectly estimate µ H 72 in ungauged sites.

Study area and available data base
The study area (Fig. 1) covers the mountainous area of Switzerland.In the Swiss Alps, over 200 snow measuring stations have been installed by SLF (Swiss Federal Institute for Snow and Avalanche Research) and by Me-teoSwiss (the Swiss Federal Office for Meteorology and Climatology) (e.g.Laternser and Schneebeli, 2003).For the present study only data collected by the manual stations network of SLF were used (Fig. 1), 124 stations with altitude ranging from 230 m a.s.l. in Bellinzona station (6BE) to 2910 m a.s.l. in Felskinn station (4FK) and an average altitude of 1522 m a.s.l.At the stations, manual measurements are carried out daily by seasonal personnel, around 7.30 a.m.The measured variables are, among others, the depth of the cumulated snow cover, or snowpack, HS and the depth of the daily fallen (i.e. in the last 24 h) new snow, HN, together with the new snow water equivalent SWE for HN>0.1 m.  winter 1984, 3UI, +140 m in 1983, 5OB, +330 m in 1973and −190 m in 1979, 7CO, +420 m in 1994.Because here elevation is investigated as one of the main driving factor of heavy snowfalls, a considerable change in altitude cannot be neglected.Therefore, the data from these stations were divided into two sub samples, before and after the relocation.This is reasonable, because stationarity in time of the distribution of H 72 is observed in practice (shown further on).The list of the stations and the related data availability is reported in Table A.1 in the Appendix A, where stations after relocation are indicated with an apostrophe after their name (e.g.2ME and 2ME').Therein, it is reported the number of available years, i.e. the number of observed values of H 72 .
3 Regional approach to frequency estimation of H 72

Estimation of H 72
For all the stations and the available years the greatest annual H 72 is evaluated, using automatic extraction procedures (all the mathematical elaborations required in the present paper have been carried out in Matlab ® environment, version 6.5).According to the Sp, the value of H 72 is calculated by evaluating the "increase in snowpack depth during a period of three consecutive days of snow fall".Thus H 72 is taken as the greatest observed value of the positive difference in snow depth calculated using a three day wide window, moving by one day steps.The measured values of H 72 in the whole region range from a least value of 0.01 m to a greatest value of 2.3 m.The average sample value for Switzerland is µ H 72av =0.56 m.In Table A1, the sample average for each site i, µ H 72i is reported.

Stationarity of H 72 in time
To carry out the regionalization procedure and distribution fitting of H 72 , it is necessary to hypothesize stationarity of the series in time.Stationarity of SWE and HS for Switzerland has been recently investigated by some researchers (e.g.Rohrer et al., 1994;Laternser and Schneebeli, 2003).Particularly, Laternser and Schneebeli (2003) considered data from the same set of snow stations as here, for the period 1931 to 1999.They investigated HS, continuous snowpack duration, and number of days of new snow HN greater than given thresholds (HN≥0, 10, 20 and 50 cm).They highlighted a visual decreasing trend of the seasonal average of HS, together with a decrease of the snowpack duration in the last twenty years or so, more pronounced at low altitudes.Also, they found a similar decrease of the number of days with HN above a threshold, more marked at low altitudes and for the smallest thresholds.Eventually, they considered the annual maximum sum of new snow in three successive days (named therein HN 3 ), coinciding in practice with H 72 , observing no relevant trend (Laternser and Schneebeli, 2003, p. 746), apart from few low lying stations (therein, p. 745).
Here, to statistically evaluate the degree of stationarity of H 72 , Mann-Kendall test is adopted, used as a valid tool for identification of trends in hydrological and climatological data (e.g.Burn, 1994;Valeo et al., 2003;Cislaghi et al., 2005;Jiang et al., 2007).The stations showing at least 30 consecutive years of data are considered for significance of the test.These are 75 stations, evenly distributed in space, with data set ending in 2006, apart from 11 stations, ending between 1996 and 2003.The test is performed as in the paper from Cislaghi et al. (2005), to which the reader is referred for full explanation.Here, for each of the 75 stations, the p-val of the MK test is evaluated and compared with a significance level for stationarity of the series p-val=5% and reported in Table A1.For 72 out of the 75 stations significant stationarity was found (i.e.p-val>5%).The three stations featuring low p-val are 6BE (A=230 m a.s.l., Y =61, p-val=1%), 7BR (A=800 m a.s.l., Y =60, p-val=2%) and 2TR (A=1770 m a.s.l., Y =66, p-val=2%), all showing a slightly decreasing trend.However, notice that the two former stations are placed at quite low altitudes, particularly 6BE station.This seems consistent with the findings from Laternser and Schneebeli (2003).However, eight stations out of ten below 1000 m a.s.l.(see Table A1) show p-val>5%, while eyeball analysis of p-val against altitude above 500 m a.s.l.do not show any clear influence of the latter.Station 2TR data set is analyzed for consistency, but no evident issues are found.It is concluded that the p-val therein reflects some local situation.
Eventually, the analyzed data base of H 72 seems reasonably stationary in time.The three stations showing non stationarity, 6BE, 7BR and 2TR are not excluded a priori from the regional analysis, because they carry out considerable information about heavy snowfalls, in view of the considerable number of sampled years (see Table A1).A judgment about their exclusion is taken also based on the results from the statistical tests for regionalization.

Index value approach
The index value regional method is widely adopted in the field of statistical hydrology (Burn, 1997;Katz et al., 2002).The index value can be calculated by the single site sample average (e.g.Bocchiola et al., 2003), or by the median (middle ranking value of the single site observed series, e.g.Robson and Reed, 1999) or by a site to site location parameter (e.g.Sveinsson et al., 2001).Here, the mean of the single site distribution is used, estimated by the sample average of the observed values at the specific site i, defined as where Y i is the number of years of observation and the suffix y indicates the year yth.The related standard error of estimation is with σ H 72i sample standard deviation of H 72 at the specific site i.The scaled value of H 72 at each specific site i is therefore defined as The symbol F i indicates the cdf of H * 72 at site i.The average of H * 72i is obviously 1 and the remaining moments need to be estimated from data.The key hypothesis for the regional approach is that F i is the same at each site i.

Identification of homogeneous regions
The first step in regional analysis is the definition of the homogenous regions, where homogeneity of the function F i holds.This requires first partitioning of the considered stations into sub-groups, or regions.This is done by grouping techniques, including Principal Components Analysis (see e.g.Rohrer et al., 1994;Baeriswyl and Rebetez, 1997;Bohr and Aguado, 2001), region of influence assessment (ROI, e.g.Castellarin et al., 2001), or seasonality measures (e.g.Burn, 1997;De Michele and Rosso, 2002).The results of such procedure should provide grouping of stations which are spatially consistent, in order to avoid senseless regions (e.g.leopard skin patterns, see Bohr and Aguado, 2001).Traditionally, Switzerland was divided into seven climatological regions, depicted considering morphological features and political boundaries (e.g.Baeriswyl and Rebetez, 1997;Laternser 2002;Laternser and Schneebeli, 2003).Laternser (2002) directly grouped Swiss snow stations applying a cluster analysis either to daily snowpack HS or to new snow depth HN.Here, the latter is analyzed because it is more correlated to H 72 .Laternser (2002) calculated correlation level Cor, defined as the correlation of time series from pairs of stations, and used it to define groups of maximum similarity.He proposed for HN a hierarchical division up to 14 regions, depending on the chosen value of Cor.Here, this division is used to iteratively evaluate the homogeneity of the H 72 distribution.A key step in regional analysis is the evaluation of the index value in sites that are not gauged, e.g. based on some proxy data from site morphology (see e.g.Bocchiola et al., 2003).As far as H 72i is concerned, regression analysis against altitude A i is used (e.g.Barbolini et al. 2002).Therefore, homogenous scaling of the average at site i µ H 72i with altitude A i is used as a further criterion for homogeneity.For robustness of the analysis, only stations featuring at least 10 years of measured data were considered (e.g.De Michele and Rosso, 2002), resulting in N=114 sites, evenly distributed in the region (see Fig. 1 and Table A1).
Starting from all Switzerland taken as homogeneous and considering an increasing number of regions according to the grouping proposed by Laternser (2002), homogeneity tests were iteratively carried out on the considered regions.We also carried out cluster analysis on H 72 for verification, finding in practice equivalent results as in Laternser (2002).The least number of regions was taken providing a reasonable degree of homogeneity according to the chosen criteria (i.e.homogeneity of function F i and scaling of µ H 72i with altitude A i ) and a spatially consistent pattern, as reported.
Eventually, 7 homogeneous regions were defined, as reported in Fig. 1 (Bianchi and Gorni, 2007).The region numbers from 1 to 5 are coincident with those from Laternser (2002, Fig. 4.7 and Fig. 4.8).This mirrors the clustering based on 5 groups as reported therein (Fig. 4.8, Correlation level Cor=0.33 or so).Some slight differences are introduced here for an increasing level of clustering.Laternser ( 2002) introduces one more region (i.e. six clusters, Cor=0.4) by splitting region 1 in two parts (S-W and centre plus N-E).Here, preliminary analysis based on statistical tests and scaling of µ H 72i with A i did not show considerable difference, suggesting to consider region 1 as a whole.For Cor=0.43 seven clusters are obtained, dividing region 5 in E and W parts. Again, preliminary tests resulted in leaving region 5 untouched.For Cor=0.48 a number of 8 clusters is taken, splitting region 4 into 4E and 4W, as here.In this case, considerable difference is seen in the preliminary analysis, suggesting division of the region in two parts (Bianchi andGorni, 2007). For Cor=0.51, Laternser (2002) suggests 9 clusters by dividing region 3 in E and W parts.Here, no such difference is spotted and region 3 is left unchanged.Notice that in region 3 only 6 stations are available, so the present results might be less robust than in the other regions.For a slightly higher value of Cor=0.53, 10 clusters are found by dividing region 2 in two parts, here indicated as 2W and 2E.This division is here accepted after preliminary tests.
The seven regions (1, 2E, 2W, 3, 4E, 4W, 5) were then tested for homogeneity, providing acceptable results.Therefore, no further division was introduced.The proposed regions are as follows.The northern slope of the Swiss Alps is divided into 4 independent areas: region 1, the north west belt crossing the country from west to east, region 2W, the Rhone Valley, region 2E, the Gotthard Range and region 5, the northern part of Grison.The southern slope is divided into 3 independent areas: region 3, the southern valleys of east Valais, region 4W covering Ticino and region 4E covering the southern part of Grison.D. Bocchiola et al.: Regional three day snow depth in Switzerland 3.5 Discordancy measures of the single sites L-moments To test the degree of homogeneity of the distribution F i , an approach based on L-moments can be used (Greenwood et al., 1979;Hosking, 1986Hosking, , 1990).L-coefficients, i.e. ratios between L-moments from order two to four (L-Coefficient of Variation, L-CV, L-Skewness, L-SK and L-Kurtosis, L-KU) are used in regional frequency analysis (e.g.Hosking and Wallis, 1993;Burn, 1997) where the suffix t indicates transposition, and The term ūi is a vector containing the average of the terms in vector u i in the N sites.Large values of D i (D i >3, Hosking and Wallis, 1993) indicate sites that are most discordant from the group and are worthy of investigation.
In Fig. 2 the L-coefficients charts are reported for the seven regions.The results of the test are reported in Table A2.Results are reported for the seven regions and also for Switzerland taken as a homogeneous region.We observed D i >3 for five stations, evenly distributed nationwide (1MN, 4FK, 5SE, 5VZ and 6BE).In region 1, it is D i >3 for one only station (Interlaken, 1IN, 580 m a.s.l., Y =50, D i =3.95) out of 36.Because this station is well inside the boundaries of region 1 (Fig. 1) this might be due to its low altitude.Station 5LQ (520 m a.s.l., Y =60), near to the border with region 5, shows D i =2.95.However, analysis of µ H 72 against altitude A suggests that 5LQ is suitable to be clustered in region 1.Also, in view of its low altitude, some uncertainty is expected.For regions 2E, 3 and 4 (E and W), it is always D i ≤3.Region 2W, including 11 stations, shows high D i in station 4VI (Visp, 660 m a.s.l., Y =61, D i =3.86), close to the boundary between region 2W and region 3, representing the north to south alps divide (Level 1, Cor=0 for clustering in Laternser, 2002, indicating a first strong division between north and south of the Alps).Because station 4VI is in a valley northern of the divide it seems that it should be clustered in region 2W.Region 5 shows D i >3 for two stations (Sedrun, 5SE, 1420 m a.s.l., Y =38, D i =3.84 and Valzeina, 5VZ, 1090 m a.s.l., Y =21, D i =3.06) out of 29.However, For the four stations closest to 5SE, one has D i ≤3 (5DI, D i =0.51, 5CU, D i =0.39, 5FU, D i =2.90, 5ZV, D i =1.96).The same applies for the four stations closest to 5VZ (5PU, D i =0.50, 5KU, D i =0.14, 5SA, D i =0.78, 5KR, D i =0.14).Also, analysis based on parent distribution fitting and Wiltshire test for region 5 did not show considerable heterogeneity of 5SE an 5VZ (in Sect.3.7, see Table A2).This suggests that the considered stations might be reasonably included in region 5.

Homogeneity test based on L-moments
To test homogeneity, the stationarity in space of the Lcoefficients can be evaluated (Hosking and Wallis, 1993).The hypothesis is that the value of the L-coefficients is the same for each station inside the region.This is verified testing whether the inter site sample variability of the Lcoefficients is compatible with the chance variation due to the statistical properties of H * 72i .First, one defines where L-CV av is the group mean value of L-CV.This variability index V 0 , extracted from the available sample of H * 72i , can be compared with its value when it is calculated from a sample of the random variable H * 72i at the same stations, in the hypothesis that they all feature the same value of L-CV.To provide a sample from such hypothetical population, numerical simulation is required.This is carried out in the hypothesis that L-CV is constant in space, to provide the expected average of V 0 , µ V 0 and its standard deviation σ V 0 .The test statistic is then According to Hosking and Wallis (1993), if H 0 <1, the value of L-CV in the considered region can be reliably viewed as constant, while H 0 ≥2 means secure heterogeneity in space.For 1≤H 0 <2, more uncertainty is expected.It is also possible to construct two other heterogeneity measures based on L-coefficients.
and a measure V 2 , based on L-SK and L-KU H 1 and H 2 can be calculated as in Eq. ( 8).Statistics based on V 1 and V 2 are weaker than V 0 in discriminating between homogeneous and heterogeneous regions (Hosking and Wallis, 1993), while V 0 is particularly relevant, because L-CV D. Bocchiola et al.: Regional three day snow depth in Switzerland Considering Switzerland as a whole, one has H 0 =17.5, H 1 =5.4 and H 2 =2.66, thus indicating strong heterogeneity.Region 1 shows a slightly negative (see e.g.Rao and Hamed, 1999) value of H 0 and values of H 1 and H 2 positive and smaller than 1.Region 2E shows values of H 0 , H 1 and H 2 positive and smaller than 2, while region 2W has H 0 negative and H 1 and H 2 positive and smaller than 1.Region 3 shows negative, but smaller than two values for H 0 , H 1 and H 2 .Region 4E shows a high value of H 0 =4.75, but with H 1 and H 2 positive and smaller than 2. Region 4W shows a negative value of H 0 (but smaller than one) and also negative values of H 1 and H 2 .Eventually, region 5 shows negative values (smaller than two) for H 0 , H 1 and H 2 (both smaller than one).In general, the results of the H 0 (H 1 , H 2 ) test seem to accept the hypothesis of homogeneity of the considered regions.Only for region 4E the value of H 0 , based on L-CV, suggests a certain degree of heterogeneity, while V 1 and V 2 suggest more homogeneity.

Distribution based Wiltshire test
Here, the distribution based Wiltshire test is performed (Wiltshire, 1986).This tests the accuracy in representing the single site values of H * 72i using a single parent, or regional distribution.It represents in practice a site to site cross-validation of the regional parent distribution based on the evaluation of the frequency (i.e. return period) of the at site observed samples.Considering a random sample of size m (x 1 , x 2 . . ...x m ) from a given distribution, the value of the non exceedance probability of each value k in the sample G k =G(x k ) can be calculated.The resulting values (G 1 , G 2 . . .G m ) should show a Uniform (0:1) distribution, apart from random deviations.If a new sample (y 1 , y 2 . . ...y m ) from a different, unknown distribution is at hand and the distribution G is used to infer the sample frequency G(y k ), the resulting set of frequencies (G 1 , G 2 . . .G m ) will differ from the Uniform (0:1) distribution by an amount greater than expected for random deviations only.
Here, if the single site distributions fit the parent GEV distribution, their sample frequencies fit a Uniform (0:1) distribution.A test statistic is introduced to evaluate this fitting, as follows.First, the sample frequencies of the observed values for the single site i and the year y, G i,y are evaluated by the parent distribution.Then, their deviation from the centre of gravity of the U(0:1) is evaluated as The test statistic is then where G i,av is the average of the term in Eq. ( 11) over the Y i years of station i and p i is its sampling variance, p i =1/(12Y i ).G av is the average of G i,av over all the stations.
The R statistic is distributed as a χ 2 , with N−1 degrees of freedom.In Table A2, the values of G i,av for the considered stations are given, while in Table A4 the average values G av are given for the regions.The p-val of the corresponding χ 2 distribution is also reported.Here and in the following, we choose a standard reference level for acceptance of the test's hypothesis of α=5%.Of course this is arbitrary, and different results may be obtained using other reference levels (e.g.α=10%).Some subjectivity is entailed into such approach, but this is a part of the regional approach and, more generally, of use of clustering techniques to interpret physically based spatial variables.For Switzerland, p-val<10 −4 , indicating unsuitability of taking the parent distribution at single sites.For all the regions but 4E one has p-val>5%, indicating a considerable confidence in taking the parent distribution to represent the at site variability of H 72 .For region 4E one has p-val=1%, so indicating a lower, albeit non negligible, degree of confidence.

Evaluation of the regional distribution
Here, the sample frequency of H * 72 is accommodated using a General Extreme Value (GEV, including EV1, or Gumbel) distribution.The GEV quantiles featuring T -years return period are evaluated as with y T Gumbel variable, y T =− ln(− ln((T −1)/T )).The parameters of the considered distribution, ε, location α, scale and k, shape, are estimated according to the L-moments approach (e.g.Kottegoda and Rosso, 1997) applied to the sample obtained by grouping the H * 72 samples in each region.Use of L-moments is consistent with the homogeneity tests in section 3.5, and further with the literature concerning estimation of extreme values distributions, particularly GEV, as they provide robust parameters estimation because they are less sensitive to rare events with very high values, which may hamper estimation using other methods (e.g.Hosking, 1990;Stedinger et al., 1992;De Michele and Rosso, 2001).
In Fig. 3, fitting of the tested distribution to the sample frequency curve is shown (APL plotting position F =(ord−0.65)/n tot , with ord ranking position in decreasing order and n tot size of the sample).Also, the confidence limits are reported, for α=5%.These are evaluated as H * 72α (T )=H * 72 (T )± α σ H * 72 (T ), where α is the 1-α/2 quantile of the standard Normal distribution, and σ H * 72 (T ) is the standard deviation of the estimated value of H * 72 (T ), calculated as (De Michele and Rosso, 2001) .( 14) The value of σ H * 72 (T ) depends on the whole sample size n tot .When single site distribution fitting is carried out, Eq. ( 14) still holds replacing n tot with the observed years, Y i .
Inter site dependence inside a region could result in poor distribution fitting, and consequently into decreased accuracy of the T -years growth curve quantiles, i.e. in a value of σ H * 72 (T ) greater than that in Eq. ( 14) (e.g.Hosking and Wallis, 1988;Bayazit and Önöz, 2004).This is turn would result in decreased accuracy of the at site T-years quantile H 72i (T ), i.e. increase of its standard deviation, σ H 72i T .However, for moderate inter site dependence, homogeneity tests are not substantially affected (e.g.Lu and Stedinger, 1992a) and GEV distribution parameters are well estimated (e.g.Madsen et al., 1997;Salvadori et al., 2007).Even in case of inter site dependence, σ H * 72 (T ) for regional analysis is much smaller than for single site analysis (e.g.Lu and Stedinger, 1992b;De Michele and Rosso, 2001;Katz et al., 2002).Bayazit and Önöz (2004) studied the increase of σ H * 72 (T ) due to inter site dependence for Gumbel distribution, mainly used here, showing that for moderate average inter site dependence (correlation coefficient ρ≤0.3 or so), σ H * 72 (T ) changes slightly (less than 20% or so), independently of T .
Here, preliminary correlation analysis for the highlighted regions indicated moderate levels of average inter site correlation (ρ≈0.2-0.3, or less), likely to result into possible underestimation of σ H * 72 (T ) in the order of 15% or less.Notice that in spite of the presence of inter site correlation, regional analysis based on extreme value distributions is used to evaluate quantiles of peak floods (e.g.GREHYS, 1996;Burn et al., 1997;Castellarin et al., 2001;Skaugen and Vaeringstad, 2005;Büchele et al., 2006), precipitation extremes (e.g.Buishand, 1991;Cong et al., 1993), daily SWE (e.g.Skaugen, 1999;Bocchiola et al., 2007) and three day snow depth for avalanche hazard mapping (e.g.Barbolini et al., 2002;Barbolini et al., 2003;Bocchiola et al., 2006), because it provides a considerable net gain in accuracy.Here, it will be shown that proposed regional approach results in a decrease of σ H 72i T for the reference return period T =300 years up to 300% or so with respect to the single site analysis.Therefore, it seems reasonable to preliminarily neglect inter site dependence.
In Table A4 the parent distributions parameters are reported.Three goodness of fit tests are carried out, namely Anderson-Darling (AD), Kolmogorov-Smirnov (KS) and χ 2 (see e.g.Kottegoda and Rosso, 1997), reported in Table A.4.For Switzerland, the proposed GEV distribution shows p-val>5% for AD and KS and p-val=1% for χ 2 .The dimensionless sample for the whole Switzerland can be reasonably well accommodated by a GEV distribution for the purpose of analyzing the extreme values of H 72 , but its use at single sites is more questionable (Wiltshire test in Sect. 3.7).For all the regions, p-val>5% are obtained.Only for region 4E, p-val=1% is obtained in the χ 2 test.
In region 1 some issues were found in stations 1IN and 5LQ, indicated in Sect.3.5 for high values of the D i statis-tic.Preliminary analysis based on Wiltshire test and plotting position showed that the distribution of H * 72i therein is considerably different as compared with the regional one (shown in Bianchi and Gorni, 2007), maybe due to low altitude of the two stations (1IN, A=580 m a.s.l., 5LQ, A=520 m a.s.l.).Thence, 1IN and 5LQ stations were not considered for parent distribution evaluation.In region 2W, a similar issue was found in station 4VI, showing a high value of D i .Again here an analysis based on Wiltshire test and plotting position showed that the distribution of H * 72i is different as compared with the regional one, maybe due to the low altitude of the station (4VI, A=660 m a.s.l.).The station is then withheld for regional distribution evaluation.In region 4W the station 6BE presented a similar issue.As reported, therein some non stationarity is detected.Analysis based on Wiltshire test and plotting position showed that also therein the single site distribution is different from the regional one.In view of its low altitude (6BE, A=230 m a.s.l.) and of the results here reported, the station is withheld for parent distribution evaluation.In Table A2, the withheld stations for evaluation of the parent distribution in each region are reported.Notice that these are all below 700 m a.s.l. or so (see Table A1).

Dependence of µ H 72 on altitude
In Fig. 4, µ H 72 is plotted against altitude.The results of weighted (on sample size Y i ) regression analysis are shown in Table A5, also for Switzerland as a whole.Significant (at a reference confidence level α=5%) increase of µ H 72 against A is observed in all regions but 2W and 5 (Fig. 4a, b).Notice that the observed dependence is of course representative of an average situation.For instance, in region 1 here reverse gradients have been observed in particular winters, as a result of certain weather situations and local topography (Rohrer, M., reviewer's communication, June 2008).However, these single occurrences are usually overridden by the viceversa case.For region 5, the findings here are consistent with e.g.Laternser (2002, p. 90).He reports on a study about the gradient of HN against altitude based on 2388 events in nine snow stations in the Davos region (in region 5 here), where no significant either positive or negative gradient exists.Concerning region 2, notice the considerable difference between the regions 2E and 2W (Fig. 4b), giving further evidence for their splitting.This independence (2W and 5) on altitude might be explained observing that the majority of the stations in these regions is above 1300 m a.s.l. in the inner-alpine dry region, thus not being affected by strong orographic precipitation.
Notice that the Sp indicates a rate of increase with altitude of H 72 for T =300 years of c Sp (300)=5 cm/100 m.Here, taking the average rate of increase of µ H 72 for Switzerland in Table A5, namely c=2.01 cm/100 m, multiplied by the 300 years quantile of the GEV distribution for Switzerland in Table A4, namely H * 72 (300)=2.63,one obtains c(300)=5.28cm/100 m.Thus, the rate of increase by the Sp is valid on average, but different rates should be considered in different regions.This is consistent e.g. with Laternser (2002, p. 91), when he states that it is ". . .not appropriate to calculate HN -altitude gradients across large areas, such as the entire Swiss Alps. . .".
During the analysis, stations 7CA and 7MA, in region 4E, showed unexpected behavior.The values of µ H 72 calculated therein are considerably higher than those expected from the general behavior of region 4E (compare Fig. 4c for region 4E and Table A1 for µ H 72 and A in 7MA and 7CA), and more similar to those observed in region 5, more constant against altitude (compare Fig. 4a for region 5).For 7CA, coupled analysis of D statistic (see Table A2) and dependence of µ H 72 vs. A, considering the surrounding stations (7BR, 7PO, 7PV) seems not to evidence significant heterogeneity with respect to the other stations in region 4E.Also for 7MA, coupled analysis of D and µ H 72 vs.A for the surrounding stations (7CO', 5JU, 5BI) seems not to display strong evidence that 7MA should be moved in region 5.Eventually, one might state that the definition of the boundary between region 4E and 5 feature some uncertainty in the area between stations 7MA and 5JU.However, the two stations 7CA and 7MA are not used here to evaluate the linear dependence of µ H 72 against A for region 4E, because they are not indicative of the behavior of µ H 72 in that region.Concerning region 4W, station 6BE was considered to evaluate the relationship between µ H 72 and A. Albeit non stationarity was found therein and the station was not used in parent distribution fitting, one can still assume that the long term mean µ H 72 is representative of the snow depth to altitude relationship.Because this conjecture seems well respected (see Fig. 4c), the mean of H 72 in 6BE is held as representative.The same reasoning applies to region 1, where stations 1IN and 5LQ were used for depth to altitude relation evaluation, as well as for region 2W, where station 4VI was similarly held.In Table A5 the standard deviation of µ H 72 when estimated from altitude is also reported, namely σE[H 72 ] .

Final remarks about the chosen regions
Comparison of the results here with the former study about clustering of snow stations by Laternser (2002) shows substantial agreement.Also, the regions proposed here can be compared to those by Baeriswyl and Rebetez (1997), based on monthly precipitation measured in 47 stations in Switzerland for the period 1961-1980.Using Cluster analysis and D. Bocchiola et al.: Regional three day snow depth in Switzerland PCA they found a number of seven climatological regions (Baeriswyl and Rebetez, 1997, p. 37, Fig. 1).Region I therein includes roughly the regions here named 3 and 4 (E and W).Region II therein coincides roughly with region 5 here, while region IV therein coincides with region 2W here.More difference is seen concerning the north western area, region 1 here, split in three parts by Baeriswyl and Rebetez (regions, V in the S-W, region VI in the centre and region III in the N-E).However, for this area the comparison is less suitable, because therein rainfall is of interest but avalanche hazard is not.Also, visual comparison with the seven traditional snow climatological regions of the Swiss Alps (see e.g. in Laternser and Schneebeli, 2003) shows reasonable consistence with those here proposed.
In region 1, a GEV distribution holds, featuring a positive value of k=0.12.Usually, either negative or null values of the k parameter are adopted (the latter indicating EV1 type one, or Gumbel distribution) for evaluation of the distribution of annual maxima of precipitation (e.g.Menabde and Sivapalan, 2000;Bocchiola et al., 2006).However, here poor fitting is obtained with a null value of k (i.e. using EV1, shown in Bianchi and Gorni, 2007).Regions 2E and 2W are separated from each other and from region 3 also in former studies.In all 2E and 2W, Gumbel distribution is well suited.Region 3 shows only 6 snow stations, thus probably requiring deeper investigation in the future.Region 4 is very likely to be divided in 4E and 4W, both featuring an EV1 parent distribution.Region 4E shows less homogeneity.However, no station shows considerable heterogeneity according to the D i test as well as to the depth to altitude relationship, and it seems reasonable to assume homogeneity of the region.Region 5 displays a good level of homogeneity and acceptable accommodation using a EV1 distribution.Also, the very weak dependence of H 72 against altitude seems consistent with previous works.
4 Accuracy of regional frequency estimation of H 72 4.1 Single site average µ H 72 A major deficiency for single site T -years quantiles evaluation is the estimation of the single site sample average µ H 72i , because it suffers from uncertainty due to scarce sample size.One is interested in knowing the least necessary number of samples, i.e. of sampled years, Y l , for reliable assessment of µ H 72i (e.g.Bocchiola et al., 2006).One can define the percentage level of accuracy of the estimate µ H 72i by its scaled mean square error of estimation σ * µi as is bigger than in the regional case (reg, use of n tot ), the final standard deviation, σ * T i is also bigger.For instance, considering 50 sampled years, close to the average here, the single site value of σ * T i is greater than the regional value by a factor ranging from 310% in the worst case of region 5 to 200% or so in the best case of region 3.As reported in Sect.3.8, the estimated values of σ * T i might be slightly underestimated due to inter site dependence, which is here preliminarily neglected.However, this seems reasonable in view of the considerable gain in accuracy, as reported in Fig. 5.

Discussion
The proposed results are valuable in the field of avalanche hazard mapping.The values of H 72i (T ) for T =30 and T =300 at a given site can be used to estimate the corresponding snow depth at release.This is used as an input to avalanche dynamics models, that can calculate the avalanche volume, the length of the runout zone and the impact pressures (e.g.Salm et al., 1990;Bartelt et al., 1999;Barbolini et al., 2002;Ancey et al., 2004), used for the evaluation of the red and blue zones, according to the Sp.Further, the proposed approach is valuable for use of dynamics models including mass uptake (e.g.Sovilla and Bartelt, 2002;Naaim et al., 2003;Eglit and Demidov, 2005;Sovilla et al., 2007), because the estimated values of H 72i (T ), explicitly depending on altitude, can be used to set initial conditions for erodible snow cover along the avalanche track (e.g.Bianchi Janetti and Gorni, 2007;Bianchi Janetti et al., 2008), of primary importance for mass uptake calculation (e.g.Sovilla et al., 2006).The estimated extension of the red and blue zones using any model is considerably influenced by the several uncertainties in the procedure, including model parameters, definition of release area, and H 72 , among others.The currently adopted approaches for the sensitivity analysis of the avalanche hazard mapping exercise include manual assessment (e.g.Bartelt et al., 1999;Sovilla and Bartelt, 2002;Bianchi Janetti et al., 2008), inter-model comparison (e.g.Barbolini et al., 2000), Monte Carlo Simulation (e.g.Natale et al., 2001;Barbolini et al., 2002), and statistically driven model calibration (e.g.Ancey et al., 2003Ancey et al., , 2004;;Naaim et al., 2004).To evaluate the deal of uncertainty in hazard maps brought by evaluation of H 72 , one can use a probabilistic assessment of the distribution of the input value H 72i (T ) (e.g.Barbolini et al., 2002).
Here, the T -years quantile can be approximately taken as normally distributed (e.g.De Michele and Rosso, 2001) with average H 72i (T ) and standard deviation σ H 72i (T ).The higher the latter, the higher the contribute to uncertainty brought by H 72 (T ).The results here show that the regional approach is more accurate (i.e.σ H 72i is considerably smaller) than single site analysis.Because the present approach to hazard mapping according to Sp is based on single site analysis of H 72 , use of the regional approach may result in decreased uncertainty of the hazard maps according to the Sp.
The regional approach is valuable also for use in ungauged sites inside a homogenous region, where local sample mean is not available.Therein, one can estimate the single site average µ H 72i , together with its standard deviation σE[H 72 ] , using an indirect approach, resulting in estimation of H 72 (T ) with quantified accuracy.For instance, regionally valid distributed maps of µ H 72i can be made available (as e.g. in Carroll, 1995 for SWE), and regression equations could be found for µ H 72i against topography or climatic attributes like those based on altitude here proposed.

Conclusions
Albeit Switzerland has relatively long time series of snow data, these are still short compared to the reference return periods for avalanche hazard mapping, i.e. 300 years.An attempt to apply a regional approach seems therefore warranted.Here, Switzerland is divided into seven homogenous regions with respect to H 72 .It is shown that the distribution of its scaled values is acceptably homogeneous in each of these regions.The scale factor, or index value, is given by the single site sample average.The observed frequencies are well accommodated by GEV distributions, including EV1, i.e.Gumbel, and the standard deviation of the T -years quantile is explicitly evaluated.Because the regional distributions are based on greater samples than the single site distributions the former are more reliable than the latter, particularly for high return periods.The estimation of the index value can be carried out by snow depth observation for a number of years, according to the required accuracy of the final T -years quantile estimates.
The approach holds even when the index value is estimated using indirect methods, if a measure of its accuracy is given.Regionally valid regression equations against altitude are provided, where significant, useful in ungauged sites.The proposed approach seems valuable for Switzerland, because it provides a nation wide valid tool for H 72 estimation inclusive of well defined accuracy, necessary for sensitivity analysis of avalanche hazard maps, according to the reference methods used in the present literature.Some issues are highlighted concerning the lowest stations, i.e. below 700 m a.s.l.However, low stations are of less interest as far as avalanche hazard is concerned.Also, some uncertainty still holds at the boundary between regions, where some degree of heterogeneity is found.In such cases, more detailed studies should be carried out to define the proper strategy for snow depth design for avalanche hazard mapping.Further studies might be devoted to evaluate the possible effect of inter site dependence on the accuracy of the estimated quantiles and on the accurate assessment of the index value in ungauged sites, for instance using other predictors than altitude.

Fig. 1 .
Fig. 1.Case study area, gauging stations network and proposed regional clustering.
L-Coefficient of Variation (ratio between second and first order L-moments) L-KU L-Coefficient of Kurtosis (ratio between fourth and second order L-moments) L-SK L-Coefficient of Skewness (ratio between third and second order L-moments) N Number of stations featuring at least 10 years of measured data Regional three day snow depth in Switzerland the most complete series feature a number of 76 years, while some stations feature less than 10 years of data.Eventually, a number of n tot =4954 values of H 72 could be collected.Few stations experienced a considerable change in position (i.e. in altitude) during the observation period.Four worst cases were found, namely stations 2ME, −160 m in Because the data set dates back to 1931, the stations with www.nat-hazards-earth-syst-sci.net/8/685/2008/Nat.Hazards Earth Syst.Sci., 8, 685-705, 2008 688 D. Bocchiola et al.: www.nat-hazards-earth-syst-sci.net/8/685/2008/Nat.Hazards Earth Syst.Sci., 8, 685-705, 2008 D. Bocchiola et al.: Regional three day snow depth in Switzerland

Table A1 .
Available data base.Y i is the number of available years (winter season).Bold indicates a number of years Y i < 10.The sample average µ H 72i is also reported.Also, the p-val for Mann Kendall (MK) test is reported, where suitable (stations with at least 30 consecutive years of data).Bold indicates p-val<5%, i.e. here significant decrease.

Table A3 .
Homogeneity test H for L-coefficients.n sim =1000.V obs is the observed value of V .µ V is the expected average of V and σ V is its standard deviation.Sw is Switzerland taken as a whole.Bold indicates a value H ≥2. All fields are dimensionless.