2.1. Avalanche occurrence and runout altitude data in the French Alps

The ‘Enquête Permanente sur les Avalanches’ (EPA) describes avalanche events on ∼3900 paths in France from the beginning of the 20th century (Reference MouginMougin, 1922). The common use of EPA data is for risk assessment at the path scale (e.g. Reference Ancey, Gervasoni and MeunierAncey and others, 2004; Reference Eckert, Parent and RichardEckert and others, 2007a, Reference Eckert, Parent, Faug and Naaim2009a, Reference Eckert, Naaim and Parent2010c), but links between avalanches and snow and weather covariates (e.g. Reference JomelliJomelli and others, 2007) or with dendrogeomorphological reconstructions (e.g. Reference Corona, Rovéra, Lopez Saez, Stoffel and PerfettiniCorona and others, 2010) have also been investigated.

This study involves all the avalanches recorded in the Alpine part of the database over 64 ‘full winters’ from 1946 to 2009, i.e. 54 641 avalanches (Table 1). Following the French convention, the ‘full winter’ starts on 1 September of a given year and ends on 30 August of the following year. Although the French Alps are divided into 23 massifs for operational forecasting, here we examine larger spatial scales: the entire French Alps and two subregions, i.e. the northern and the southern French Alps (Fig. 1), with the northern French Alps representing ∼70% of the data.

Table 1.

Avalanche occurrence data: full and filtered sample

Fig. 1.

Study area. The French Alps are divided into 23 massifs in operational forecasting. Here, in addition to the entire French Alps, only two groups of massifs are considered, the northern and southern French Alps, represented in light and dark grey respectively.

For detecting time trends, EPA’s major advantage is that it contains long data series from a sample of paths for which all avalanches are theoretically recorded, instead of trying to collect all major events everywhere (e.g. in an avalanche atlas). Although the protocol has seen several changes, including a major update in 2002 (Reference Burnet and AnceyBurnet, 2006), its philosophy has remained sufficiently the same to ensure a certain continuity in the data series, at least at scales sufficiently large to smooth discrepancies. Furthermore, in order to record mainly natural and unperturbed avalanche activity, the proportion of artificial or accidental triggers is very low on EPA paths, and they are little affected by the construction of recent countermeasures. For example, Reference Eckert, Baya and DeschatresEckert and others (2010b) found that similar conclusions were found if the few perturbed paths were included in or excluded from analysis.

Following Reference Eckert, Parent, Bélanger and GarciaEckert and others (2007b, Reference Eckert, Parent, Kies and Baya2010a), aggregation of occurrence data has been performed at the township scale. Hence, we define a_jt as the number of avalanches in the township j during the winter t, where j ∊ [1,M], t ∊ [t ₀, t ₀ + T _obs − 1]. M is the total number of townships in the region studied, T _obs the length of the observation period and t ₀ the first winter considered. We also define c_j as the number of surveyed paths in the township j, which is taken as constant over the entire considered period, and equal to the present number fixed during the last data collection protocol update.

For occurrences, the predominant source of remaining error is missing events. A simple test procedure (Cemagref ETNA, 2008) has been implemented to discard the township/winter couplets where the observed avalanche count is statistically ‘abnormally’ low owing to undercounting by local avalanche observers. The test is based on the comparison of each annual value to the local 20 year running mean, and while this discards ∼40% of all township/winter couplets, it retains most of the observed avalanche events (90–93% depending on the considered region; Table 1).

For safety reasons, rangers do not actually measure runout altitudes, but estimate them from a distant observation point, and then plot these estimations on a map. As a consequence, recorded runout altitudes are more uncertain than avalanche counts, and may be missed because of bad visibility. Following Reference Eckert, Baya and DeschatresEckert and others (2010b), rather than altitude we use the Runout Altitude Index (RAI):

where e = exp (1), denotes the runout altitude of the avalanches i ∊ [1, N_kt ] recorded in the avalanche path k ∊ [1, L] during the winter t, and is the minimal runout altitude possible in the path k (often the valley floor). By definition, RAI = 1 if is reached. If not, it is a continuous and decreasing function of the runout altitude belonging to [0,1]. Note that, in case of climbing the opposite side, is deemed to have been reached.

As it is a dimensionless scaled variable, RAI permits the comparison of runouts between avalanche paths of different size, aspect, altitude, etc. From this perspective, it bears a relation to the runout ratio index used in avalanche engineering to evaluate extreme avalanches (Reference McClung and LiedMcClung and Lied, 1987), although without the use of the ‘beta point’ position. On the other hand, RAI gives more weight to paths where the runout altitudes reached are far above the reference value . To limit this bias, minimal altitudes as realistic as possible were chosen, and, as in Reference Eckert, Baya and DeschatresEckert and others (2010b), data quality checks were performed using several deterministic and statistical procedures to discard paths with aberrant values from the study. The retained RAI data represent 35% of the original sample, on about 2600 paths, ∼1650 in the north and ∼950 in the south subregions. This loss of information is deemed necessary to ensure we obtain robust results.

2.2. Flow regime, and snow and weather covariates

For avalanche forecasting, Météo-France employs two numerical models, SAFRAN (Reference Durand, Brun, Mérindol, Guyomarc’h, Lesaffre and MartinDurand and others, 1993) and Crocus (Reference Brun, Martin, Simon, Gendre and ColéouBrun and others, 1989, Reference Brun, David, Sudul and Brunot1992), to assimilate all available snow and weather information and to then simulate meteorological parameters, snow and cover stratigraphy at various altitudes, aspects and slopes according to these data and physical rules. The models have been used for retrospective climate investigations for a period starting in winter 1958 (Reference Durand, Laternser, Giraud, Etchevers, Lesaffre and MérindolDurand and others, 2009a,Reference Durand, Giraud, Laternser, Etchevers, Mérindol and Lesaffreb). At the winter (15 December–15 June) scale, snow and weather covariates from these analyses have been successfully related to simple avalanche activity indices using regression models that represent trends and high/low peaks well (Reference Castebrunet, Eckert and GiraudCastebrunet and others, 2012). In this paper, we relate our avalanche data to two synthetic variables derived from this work: the SAFRAN mean winter temperature and Crocus mean winter snow depth at 2400 m averaged over the four slope expositions (north, south, east, west). These are denoted by Temp _t and Snow _t , respectively. They represent mean behaviour at large spatial scales, such as the whole French Alps and the north and south subregions, better than a single snow and weather point observation series, whose selection over others introduces difficulties.

Since 1973, the flow regime has been recorded in the EPA, and, as an additional covariate to explain the annual fluctuations of avalanche activity, we consider the annual proportion of powder- and mixed-snow avalanches computed on the filtered runout altitude sample, PSA _t . The

annual proportion of purely dense snow avalanches is simply 1 – PSA _t . Temp _t and Snow _t , and the flow regime proportion PSA _t , have been processed with simple intermediate- (5 year) and low (15 year)-frequency running mean filters to highlight structured patterns at different frequencies to be compared to those inferred from the avalanche data as detailed in Section 3. The 1973–2009 flow regime proportion PSA _t has also been adjusted with a simple linear regression, as this has been found to be supported by the data.

3.1. Hierarchical modelling versus empirical estimation

Empirical estimates for the mean annual avalanche number per path, , the mean annual RAI, , the mean annual runout altitude, , and the annual proportion of high-magnitude avalanches reaching the valley floor, , can be obtained using

with the ‘hat’ indicating an estimated quantity, in contrast to the generally unknown, true value. Computing involves taking into account the proportion percobs _t of non-missing township/winter couplets for the winter _t , with the underlying assumption that the missing township/winter couplets behave like the observed ones. This may be critical when percobs _t is too low to consider that the regional behaviour is well captured in the available data. is derived from using , the mean altitude of the valley floor in the region studied. Finally, the indicator function if the minimal runout altitude is reached and 0 if it is not.

A Spearman’s rank correlation test between the chronological order and the magnitude of the empirical estimates indicates that, for the 5% significance level, the hypothesis of no significant correlation is only rejected for runout altitudes in the southern Alps. Hence, major non-stationarities do not exist in most of the empirical, filtered series. For instance, all filtered occurrences series are declared stationary by the test, whereas this is not the case for two non-filtered occurrence series out of three. This highlights the homogenization effect of the filtering procedure and indicates that refined methods have to be employed to extract significant temporal patterns from these data. As stated in the introduction, to obtain the common effect from a sample of paths and to depict associated trends and change-points, instead of working with empirical estimates only, we therefore perform time-series analysis within a parametrical hierarchical Bayesian modelling approach to:

In contrast to a simple empirical approach, a hierarchical Bayesian modelling approach is richer, allowing consistent inference of quantities of interest such as trends and change-points in short time series (in our case 65 years long). The approach compensates for time by dependent repetitions across space (paths/townships) and by assumptions regarding data distributions, form of the investigated trends, types of change, etc. On the other hand, these modelling assumptions may influence inference. Figure 2 shows that, in our case, they are not too constraining. Indeed, our model estimates are very close to the empirical ones, and generally not distinguishable as soon as the associated standard error (measured by the 95% credibility interval in Fig. 2) is considered. This indicates that the parametric framework used remains flexible enough to infer annual patterns in the data. The exception concerns very low annual occurrences where empirical estimates are strongly affected by undercounting by observers, while the model uses the spatiotemporal structure of the full dataset instead of only the annual percentage of missing township/winter couplets, providing more robust estimates. Therefore, the analyses made in the rest of the paper use model estimates instead of empirical estimates. The next subsections detail how these model estimates are obtained, and how major temporal patterns are isolated in the different series.

Fig. 2.

Empirical estimates vs model estimates (model M₀), entire French Alps. (a) Number of avalanches per path, . (b) Mean runout altitude, .

3.2. Extracting the mean avalanche number per winter and path

Following Reference Eckert, Parent, Kies and BayaEckert and others (2010a), the annual avalanche counts a_jt are modelled with a non-homogeneous Poisson process inspired by spatial epidemiology (Reference Elliott, Wakefield, Best and BriggsElliott and others, 2000), with parameters λ_jt , j ∊ [1, M], t ∊ [t ₀, t ₀ + T _obs − 1] summarizing the local annual avalanche activity:

A standardization by the number of avalanches e_j , j ∊ [1,M] expected in each township j under the hypothesis of space–time homogeneity is used to isolate RR_jt ∊ ]0, + ∞[, the ‘relative risk’. It indicates if the observed number of avalanches per path is significantly greater or lower than that for a mean township during a mean winter in township j during winter t:

where e_j is evaluated by weighting the mean annual number of observed avalanches per path by the number of paths c_j undre survey in the considered township:

Further decomposition of RR_jt into spatial and temporal effects is undertaken assuming full separability between space and time is realizable:

The locally unstructured term, v_j , takes into account any strong local excess or deficit in the local relative risk, whereas the structured spatial component, u_j , models the inter-township smooth signal. These spatial terms are not considered further in this paper. The annual term, gt, represents the interannual fluctuations in the relative risks that similarly affect all the townships. The model estimate of the mean avalanche number per path and winter in the considered region is then

3.3. Extracting the mean annual runout altitude

Following Reference Eckert, Baya and DeschatresEckert and others (2010b), the RAI is modelled using a mixture of two independent distributions. RAI _{ikt 1} is a Bernoulli variable taking the value 1 if avalanche i occurring during winter t on path k reached the minimum altitude, and 0 if not. Hence, p_t ∊[0, 1] is the annual probability of reaching the minimum altitude. RAI _{ikt 2} models all the smaller events by a beta distribution, with an annual parameter pair (α_t , β_t ), α_t > 0, β_t > 0:

The RAI model is simpler than the occurrence model in that it does not take into account spatial effects, but more complex because a mixture model is necessary to fit the data structure. From the linearity of mathematical expectancy, m_t , the annual mean RAI is

which gives

Consequently, the triplet (p_t , m_t , β_t ) fully characterizes the RAI annual distribution and may be compared to g_t in the occurrences model. The model estimate of the mean runout altitude in winter t, , derives from the model estimate of the RAI, :

3.4. Time trend modelling

We model g_t , m_t and p_t as latent variables, i.e. as model unknowns that behave as parameters with regard to the data, but whose distributions are indexed by (hyper-) parameters, i.e. within a hierarchical framework (e.g. Reference WikleWikle, 2003; Reference Banerjee, Carlin and GelfandBanerjee and others, 2004). Note that, on the other hand, the are taken as exchangeable parameters so that their possible smooth trend is not modelled. It has been checked after inference that their interannual variability is low, allowing the part of the smooth signal they capture to be neglected, even for the evaluation of trends in high-magnitude avalanches (see Section 3.6). Furthermore, we consider two different time trend models for g_t , m_t and p_t , so as to distinguish different changes affecting the signal at different frequencies.

3.5. Bayesian inference and shrinkage effect

Inference was implemented using Bayesian Markov chain Monte Carlo (MCMC) methods (Reference BrooksBrooks, 1998; Reference Gilks, Richardson and SpiegelhalterGilks and others, 2001) which are quite convenient for the complex models used but require careful handling (e.g. when ensuring convergence). Hence, for each analysed series, inference robustness has been checked using tests based on starting different simulation runs at different points of the parameter space (Reference Brooks and GelmanBrooks and Gelman, 1998). For all parameters except the number of jumps in model M₂ (a too high number of jumps has been penalized), poorly informative priors were used. From the joint posterior distribution of all parameters, latent variables and missing values we retained point estimates (the posterior mean), posterior standard deviations and 95% credibility intervals.

A great advantage of this framework is that the posterior distribution of any latent time series is likely to fit complex temporal patterns even with a relatively simple parametric model. For example, Figure 3 shows that the number of avalanches per winter and per path at the whole Alps scale is clearly captured with models M₁ and M₂, the level of agreement between model and empirical estimates being very good. This justifies the statement made in Section 3.1 of a limited influence of modelling assumptions on the inferred annual patterns. Furthermore, model M₁ captures linear trends before and after a nearly 10 year transition period 1978–88 whose flat shape reflects the uncertainty of the change-point date, τ. Model M₂ captures within its trend less regular behaviours such as the ‘bulge’ between 1950 and 1954, and the recent 2006–09 increase, justifying the ‘intermediate-frequency trend’ label. Uncertainties about the trends provided by the two models are similar, except when M₂ detects patterns not seen using M₁, i.e. in the 1950–54 and 2006–09 periods.

Fig. 3.

Hierarchical decomposition of the number of avalanches per winter and per path, , at the entire French Alps scale: annual signal and underlying trends. Annual estimates provided by the different models are indistinguishable, with only the green line for M₀ visible. Associated credibility intervals for the annual estimates are not shown, for reasons of clarity. Exceptional winters are detected with different thresholds.

The annual model estimates for models M₀, M₁ and M₂ are indistinguishable for occurrences in the whole Alps (Fig. 3). However, runout altitudes in the northern Alps differ, with model estimates provided by M₁ closer to the low-frequency trend (Fig. 4a). Shrinkage is spectacular for runout altitudes in the southern Alps at the beginning of the study period because of the small number of data available in this region at this time (Fig. 4b). Annual estimates provided by M₁ are then extremely close to the low-frequency trend, so the interannual variability is underestimated. Over more recent winters, the difference between the two models is reduced because the number of data is much greater. Hence, since M₀ retains greater variability, all further analyses regarding annual estimates are based on M₀ estimates, and, for the different variables, we compute the fluctuations (high-frequency signal) by subtracting M₁’s trend from M₀’s annual estimates.

Fig. 4.

Shrinkage effect: mean runout altitude . (a) Northern French Alps. (b) Southern French Alps. Empirical estimates and credibility intervals for annual estimates are not shown, for simplicity. Minimal runout altitudes are 1170 and 1448 m, respectively.

Table 2 quantitatively assesses these statements, showing the excellent correlation between model and empirical estimate and significant correlations (R = 0.43–0.85) between the annual estimates and the estimated low- and intermediate-frequency trends (the value for runout altitudes in the southern Alps is very high due to shrinkage). Fluctuations also remain strongly correlated with the annual estimates, which is not surprising as M₁ captures only the predominant low-frequency pattern.

Table 2.

Empirical correlation between annual estimates provided by model M₀ and the other terms for the entire and southern French Alps. Fluctuations (high-frequency signal) are obtained by subtracting M₁’s low-frequency trend from model M₀’s annual estimates. Correlations are evaluated for the 1946–2009 study period, except for runout altitudes in the southern French Alps for which the 1949–2009 subperiod is considered. All values are nonzero at the 5% significance level

3.6. Evaluating time trends in high return period avalanches

We may combine the different estimates to evaluate the temporal fluctuations of high-magnitude avalanches. The annual return period reaching the valley floor is

The associated low- and intermediate-frequency trends are obtained by considering instead of in Eqn (22). However, subsumes genuine change and improved precision of runout altitude records. Consequently, we have attempted to find a less biased indicator for the annual occurrence of high-magnitude avalanches.

The modelled annual distribution of the RAI can be explored by simulating a large sample (50 000 values were necessary) given , taking the percentiles of interest and using them in Eqn (16). Figure 5a shows the evolution of the runout altitudes corresponding to annual non-exceedence probabilities of 0.75, 0.84 and 0.90. If the exceedence probability of interest is higher than , the valley floor is reached. This is nearly always for a nonexceedence probability of 0.9, and ∼50% of the time for a non-exceedence probability of 0.84, which corresponds to the interannual mean of in the French Alps. A similar approach can be used to obtain the low- and intermediate-frequency trends for these percentiles, using instead of to simulate each annual distribution. Averaging over the must be done because, as noted earlier, their possible trend is not modelled. Figure 5b shows that this simplification is not too strong since a reasonable representation of the intermediate- and low-frequency trends for the runout altitude corresponding to an annual non-exceedence probability of 0.75 is obtained.

Fig. 5.

Simulation of runout altitude quantiles (entire French Alps). (a) Quantiles q₇₅, q₈₄ and q₉₀. (b) The simulated quantile q₇₅ and associated low- and intermediate-frequency trends.

Finally, the simulated annual percentiles can be combined with the annual avalanche occurrences to extract the runout altitude corresponding to a given return period. Indeed, if the return period of interest is T, then taking the percentile of the simulated RAI annual distribution and using it in Eqn (16) gives the runout altitude corresponding to the return period T. This approach has been used to study the evolution of runout altitudes and corresponding to return periods of 10 and 20 years. For higher return periods, the minimal altitude is always obtained so that little can be said about the runout behaviour of the most extreme events. Note also that the associated uncertainty levels could not be obtained fully rigorously for the underlying trends because of the approximation made while simulating given average. Finally, empirical estimates and have also been derived, combining empirical RAI annual percentiles with the from Eqn (2).

4.1. Mean avalanche occurrence

According to the M₀ estimates for , the mean annual avalanche number on a mean French Alpine path is close to 0.32 (0.32 in the north and 0.33 in the south subregions; Table 3). Interannual variability is strong, with an empirical standard deviation of annual estimates close to 0.1 avalanches path⁻¹ winter⁻¹ at the entire Alps scale, ranging up to 0.135 in the south subregion. Hence, there are considerable variations from one winter to another, and the trends at low and intermediate frequencies identified by M₁ and M₂ do not indicate marked systematic changes, capturing frac.struc = 20–25% and 24–27%, respectively, of the signal only.

Table 3.

Descriptive statistics for annual estimates, model M₀. Mean, standard deviation, minimum and maximum over the study period

A threshold of ±1.5 standard deviations highlights the winters 1950 (in fact 1950/51), 1977, 1985, 1994 and 1998 as high-activity winters, and 1947, 1948, 1955, 1963, 1972 and 2006 as low-activity winters at the entire Alps scale (Fig. 3). In both the north and south subregions, 1977 and 1985 are detected as high-activity winters, while 1987, 1994 and 1998 are detected as high-activity winters in the north subregion only. In contrast, 1950, 2008 and 2009 are detected as high-activity winters in the south subregion only, although 1950 is just below the threshold in the north subregion (Fig. 6).

Fig. 6.

Number of avalanches per winter and per path, , in (a) the northern French Alps and (b) the southern French Alps. For the southern Alps, model M₁ is fitted on the full study period and the 1946–2007 subperiod to highlight the ‘window effect’ on the low-frequency trend.

Reference Durand, Giraud, Laternser, Etchevers, Mérindol and LesaffreDurand and others (2009b) established that low avalanche activity in 1963 was due to an extremely weak snow cover. The famous avalanche cycle of February 1999, which included a major avalanche in Montroc (Reference Rousselot, Durand, Giraud, Merindol and DanielRousselot and others, 2010) and also caused widespread damage in Europe (SLF Reference DavosDavos, 2000), occurred within the highlighted 1998 winter. Similarly, the December 2008 avalanche cycle caused considerable traffic disturbances and damaged equipment and buildings in the eastern part of the southern French Alps (Reference Eckert, Coléou, Castebrunet, Giraud, Deschatres and GaumeEckert and others, 2010d). Hence, although our approach smooths the signal by cumulating avalanche counts, high/low values represent the observed fluctuations of avalanche occurrences well. More detailed analyses of the relations between high-/low-activity winters and their climatic drivers in the different considered regions are provided by Reference Castebrunet, Eckert and GiraudCastebrunet and others (2012).

Except for a concentration of low values at the beginning of the study period, which could be, despite our efforts to filter out such phenomena, a database effect, it is difficult to detect a change in the number or distribution of winters with low/high activity at the scale of the entire Alps or within the two subregions. However, the low-frequency from M₁ shows that, at the whole Alps scale, the mean number of avalanches per winter and path has increased in the first half of the study period from 0.24 in 1946 to more than 0.37 in ∼1980, and has decreased during the second half of the study period to 0.3 avalanches path⁻¹ winter⁻¹ in 2009 (Fig. 3). Both trends (b ₁/b ₂ parameters) are not fully significant at the 95% credibility level, but have a relatively high posterior probability of being positive/negative, respectively (Table 4). Transition occurs during the period 1976–85, when activity was stronger than during the rest of the study period, at ∼0.35 avalanches per winter and path. The rather smooth transition reflects the uncertainty regarding the date of change, with the best posterior estimate being , but with a relatively large posterior standard deviation of 5 years (Table 4). There is no marked difference between the variability around the trend before and after the transition period (σ ₁/σ ₂ parameters, Table 4), so that is not a change-point in variance.

Table 4.

Posterior estimates, model M₁. Mean, standard deviation and 95% credibility interval. (b ₁, b ₂) and (σ ₁, σ ₂) are the trends/standard deviations before/after the change-point τ, respectively

The transition occurs earlier ( for the north–south; Table 4) and is more marked (posterior standard deviation of 3–5 years) in the northern Alps than in the southern Alps. Before the change-point, the increase is very strong in the southern Alps (nonzero at the 95% credibility interval) and very weak in the northern Alps. After the change-point, the decrease is similar to the overall Alpine behaviour in the northern Alps (just nonzero at the 95% credibility interval), whereas a slight increasing trend is detected in the southern Alps. This latter, surprising result is driven by the 2008 and 2009 high-activity winters in this subregion, since a slight decrease agreeing with the overall result is obtained if these two winters are excluded from analysis (Fig. 6).

For the entire Alps, the shifting level intermediate-frequency trend, , from M₂ detects low activity followed by a bulge in the early 1950s corresponding to well-documented harsh avalanche winters in Europe (Reference VoellmyVoellmy, 1955). There follow three long, flat segments, one between roughly 1975 and 1988 corresponding to the period of high activity discussed above, and the other two, before and afterwards, being quite close to the average interannual activity. Finally, the model identifies a recent (since 2006) strong rise (Fig. 3). Overall, the total series is segmented into six subperiods: three that correspond to the low-frequency signal and three much shorter ones that cannot be detected with M₁. The strength of the different change-points is quantified by the posterior probability of a level shift . This highlights 1949 as a very strong change-point , while the beginning and end of the 1975–88 high-activity period are less strong, but noticeable, change-points .

In terms of north/south differences, model M₂ highlights high activity in the early 1950s in the north region only: while the northern Alps experienced a succession of harsh winters, only 1951 was severe in the southern Alps. Similarly, the beginning and end of the 1975–88 high-activity period are more visible in the northern Alps than in the southern Alps . On the other hand, in the southern Alps, 1959 is a noticeable breakpoint in the increasing trend over the first half of the study period, and the effect of the last two high-activity winters is much more visible (Fig. 6).

Hence, M₂’s results for the northern Alps are logically very similar to those obtained in Reference Eckert, Parent, Kies and BayaEckert and others (2010a) for the same region with the same model over the 1946– 2005 period. That study concluded that there has been no recent systematic evolution of the occurrence process in the northern French Alps. The current work has not only extended the spatio-temporal extent of the analysis but has permitted, with model M₁, the detection of a change-point and of a slight low-frequency trend at the entire French Alps scale and in the two north/south regions, which were hidden in the pseudo-cyclic variations highlighted in the previous study.

The northern Alps contribute ∼70% of the data and contain more homogeneous massifs than the southern Alps. Hence, their response is closer to that of the French Alps as a whole (empirical correlation coefficient between annual estimates R = 0.92) than the behaviour of the southern Alps (R = 0.64). However, significant correlations exist between the annual estimates in the northern and southern Alps (R = 0.4), reaching a maximum for the low-frequency trend (R = 0.71). The correlation is lower but remains significant for intermediate-frequency trends and for fluctuations (R = 0.46 and R = 0.34, respectively). Thus, in terms of trends, high-activity winters and the position of change-points, there is a partially coupled response between the northern and southern Alps.

In more detail, the centred standardized difference between the annual estimates in the two regions shows numerous winters with low difference in terms of relative activity, but also strong outliers (e.g. 1994 and 1998 with strong excesses in the northern Alps, and 2008 and 2009 with strong excesses in the southern Alps; Fig. 7a). The relative activity is much higher in the northern Alps before 1958 and less strongly greater in the southern Alps between 1958 and 1990. Since then, there has been a period of very variable relative activity where most of the ‘outliers’ appear. This may indicate that the north/south coupling is less strong than before.

Fig. 7.

Centred standardized north/south difference in annual estimates. (a) Number of avalanches per winter and per path. (b) Mean runout altitude. For each variable, the centred standardized difference is evaluated as , where x_. and denote the annual estimate and its interannual mean, subscripts N and S refer to the north and south subregions, respectively, and σ _d is the standard deviation of the centred difference .

4.2. Mean runout altitude

According to M₀, the interannual mean runout altitude on a mean French Alpine path is close to 1430 m (Table 3) and nearly 200 m higher in the southern Alps (1565 m) than in the northern Alps (1369 m). In the different regions, the empirical standard deviation of the annual estimates is close to 30 m, which is relatively low. A significant fraction of the temporal signal (frac.struc = 30–80%) is therefore captured by low- and intermediate-frequency trends, except for the intermediate-frequency trend in the northern Alps (16%), as discussed below.

The winters during which runout altitude was low on average are highly concentrated in the middle of the study period. Using a threshold of −1.5 standard deviations, 1970, 1971, 1976, 1977 and 1985 are identified at the whole French Alps scale (Fig. 8), compared to 1977, 1980, 1985 and 1987 in the northern Alps, and 1967, 1980 and 1981 in the southern Alps (Fig. 4). The winters when avalanches remained on average at higher altitudes are 1966, 1975, 2001, 2002 and 2004 at the whole Alps scale, and the same except 2004 in the northern Alps. For the southern Alps, 1994 and all winters since 2005 are above the +1.5 standard deviation threshold.

Fig. 8.

Mean runout altitude at the entire French Alps scale. Minimal runout altitude is 1246 m.

At the whole Alps scale, the low-frequency trend shows a clear change-point in the late 1970s, with a best posterior estimate, , similar to that obtained for occurrences and a posterior standard deviation of 4 years (Table 4). Before the change-point, the mean annual runout altitude decreased by 55 m from 1946 to ∼1980 (Fig. 8). Since then, avalanches have retreated again, reaching more or less the 1946 state in 2009. There is no clear difference in interannual variability before and after the change-point (σ ₁/σ ₂ parameters in Table 4). The increasing trend after the change-point is fully significant at the 95% credibility level (the b ₂ parameter is negative because the RAI is a decreasing function of the runout altitude), whereas the decreasing trend before the change-point is close to the 95% credibility level (Table 4). Hence, trends are well supported by data. Over the 1946–2006 period, these results are very similar to those found by Reference Eckert, Baya and DeschatresEckert and others (2010b) for the merged Alps and Pyrenees data, which is logical since the number of Pyrenean data is small compared to that from the Alps.

Similar to occurrences, Figure 4 shows that change in runout altitudes occurs earlier and is stronger (posterior standard deviation is 0.6; Table 4) in the northern Alps than in the southern Alps (, with a posterior standard deviation of 3.7 years), where it is smoother. Hence, in the northern Alps, the low-frequency trend shows a strong shift between two slightly marked increasing trends, whereas in the southern Alps the transition between a significant decreasing trend and an even stronger and significant increasing trend is more gradual. This confirms rather different behaviors at low frequency in the two regions.

At the entire Alps scale (Fig. 8), before 1990, M₂ identifies two long segments separated by a transition period lasting a few winters around 1970. Hence, there is a high-runout regime at the beginning of the study period (mean runout altitude 1435–1440 m) and a low-runout regime from ∼1972 to ∼1990 (mean runout altitude 1415–1420 m). From 1990 to the early 2000s, an increasing trend is visible, but then it stops, with avalanches again reaching lower runout altitudes during recent winters . This recent termination of the upslope retreat of large avalanches could not be demonstrated in Reference Eckert, Baya and DeschatresEckert and others (2010b) because model M₂ was not used in that study and only runout altitudes recorded up to 2006 were considered. In the northern Alps, the intermediate-frequency trend is almost ‘flat’, explaining the low fraction of the signal captured by M₂ in this region (Fig. 4a). For instance, no significant rise is inferred, whereas, in the southern Alps the recent increase is very strong (Fig. 4b), occurring through two successive levels .

The correlation between the annual runout altitude estimates in the northern and southern Alps is low (R = 0.07). Even for the low-frequency trend, there is no significant correlation between the northern and southern Alps. There is therefore a nearly fully decoupled behaviour of runout altitudes between the northern and southern Alps. Hence, the centred standardized difference between the two regions is often high (Fig. 7b), with mostly positive values between 1960 and 1992, and negative values since 1993, due to the rapid decrease of runout altitudes in the southern Alps over the recent period.

4.3. Inter-variable correlations

The winters during which avalanche occurrences were higher on average correspond quite well to the winters where avalanche runout altitudes were lower on average (R = −0.39). For instance, 1977 and 1985 are detected as abnormally harsh winters for both occurrences and runout altitudes using the 1.5 standard deviation threshold. This similarity is enhanced if one looks at trends. For example, R = −0.82 for low-frequency trends, the ‘V’-shaped evolution of mean runout altitudes with a minimum around 1980 corresponding well to the ‘flat inverted V’-shaped pattern in avalanche occurrences with a 1975–88 period of high activity. Correlation is also strong at intermediate frequency (−0.46), and low but still significant for fluctuations (R = −0.29). Hence, runout altitudes and occurrences are not independent processes at the annual timescale, an important result for hazard assessment which could not be demonstrated in the preliminary approaches of Reference Eckert, Anderssen, Braddock and NewhamEckert (2009) and Reference Eckert, Parent, Naaim, Schweizer and Van HerwijnenEckert and others (2009b) because the datasets then considered were not fully coherent in terms of spatiotemporal scales.

For avalanche occurrences, the interannual variability around the low-frequency trend is greater, the uncertainty around the change-point is higher and the trends before and after the change-point are less marked and less significant than for runout altitudes. This explains the weaker clustering of winters with a high number of avalanches around 1980 compared to the winters with lower runout altitudes. These results may be partially related to the different observation models used for the two variables (non-homogeneous spatio-temporal Poisson process versus beta-binomial temporal mixture model). However, they highlight one important difference between the temporal evolution of the two variables: the structured low- and intermediate-frequency signal is more pronounced for runout altitudes.

Differences exist in the strength of this correlation between the northern and southern Alps. The correlation is stronger in the northern Alps (R = −0.91 at low frequency, −0.76 at intermediate frequency and −0.43 at the annual scale), while in the southern Alps, only low-frequency trends are significantly correlated at the 95% level (R = −0.37). For instance, the last two winters have seen exceptionally high avalanche numbers but very few low runouts in the southern Alps (Figs 4b and 6b).

4.4. Probability of reaching the valley floor and the associated return period

An interesting output from the runout altitude model is the annual probability of reaching the valley floor (Fig. 9a). The low-frequency trend increases slightly until the 1978 change-point, since when it has decreased markedly and almost continuously until today. As already pointed out in Reference Eckert, Baya and DeschatresEckert and others (2010b), the variability about the trend is stronger after the change-point than before. This explains why, for some recent winters such as the catastrophic 1998/99 season, remained relatively high. This must be kept in mind when considering the recent very sharp retreat of large avalanches from a hazard-zoning perspective. Although the trend implies reduced risk, the increased variability makes winters with a high proportion of very large avalanches still possible. The intermediate-frequency trend clearly shows the concentration of high values around 1980, corresponding to the period of abnormally harsh winters discussed above. In addition, in contrast to the continual drop shown from M₁, M₂ shows quasi-constant values from the mid-1980s until the late 1990s followed by an extremely strong decrease between 2000 and 2006.

Fig. 9.

(a) Annual probability of reaching the valley floor and (b) associated return period at the entire French Alps scale.

Equation (22) gives the return period for reaching the valley floor. Given that occurrences and runout altitudes are modelled independently, the excellent agreement between empirical estimates, annual estimates and the two modelled trends (Fig. 9b) is remarkable. interannual mean is about 20–30 years, confirming that this variable quantifies the temporal evolution of large avalanches, but not of extreme ones. Its main patterns are a direct consequence of the behaviour of and . First comes a slight decrease due to continuously increasing values of and , perturbed at intermediate frequency by the short period of even higher values in the early 1950s. At ∼1980 there is a concentration of winters with a lot of major avalanches due to concomitant maximal values for both and , leading to minimal values of close to 10 years. Finally, between ∼1980 and 2010, a slight decrease in , combined with a very strong decrease in , leads to a dramatic rise in , increasing to nearly 50 years. M₂ suggests that this has occurred more precisely between 2000 and 2006 due to the surprising intermediate-frequency trend in over this period, and has been interrupted in recent winters because of increased values of and stabilized values of since 2006.

Although the trend in since ∼1980 is in agreement with observations for and , its magnitude is too strong to merely reflect physical reality, especially from 2000 to 2006. The recent improvement in the precision of runout altitude records following the latest EPA protocol review (Reference Bélanger, Cassayre and ElderBélanger and Cassayre, 2004) is probably partly responsible. Better maps and topographical descriptions mean that rangers now register runout altitudes more precisely. For example, a runout altitude of 1005 m on a path with a valley floor altitude of 1000 m is now recorded if that is what actually occurred, whereas previously the runout altitude would have been considered to be 1000 m, artificially inflating the proportion of avalanches that have reached their minimal possible altitude.

4.5. Runout altitudes corresponding to high return periods

The runout altitudes and correspond to return periods of 10 and 20 years, respectively. They reduce the bias in quantifying the evolution of high-magnitude avalanches, although all temporal patterns remain consequences of inferences on and , enhanced by their partial correlation. Hence, because of the high interannual variability of discussed above, there is a strong interannual variability in and (Figs 10 and 11). For these two variables, similar patterns are observed, with the difference that patterns for are more likely to be truncated by events reaching the valley floor. Indeed, and as soon as respectively.

Fig. 10.

Runout altitude corresponding to a 10 year return period : (a) entire French Alps, (b) northern French Alps and (c) southern French Alps.

Fig. 11.

Runout altitude corresponding to a 20 year return period: (a) entire French Alps, (b) northern French Alps and (c) southern French Alps. In (c) the interannual mean is .

At the entire Alps scale, for , one can detect a marked V-shaped low-frequency pattern (Fig. 10a) that parallels the one in (Fig. 8), but with a greater maximal amplitude: nearly 90 m between the beginning/end of the study period and several winters around 1980 for which . Low- and intermediate-frequency trends are slightly higher at this time, close to 1260 m, leading to a difference of ∼80 m between ∼1980 and the beginning/end of the study period. This makes a horizontal runout distance difference as high as ∼450 m on a typical 10° runout slope. Finally, there is a departure between M₁ and M₂ in the early part of the study period, with M₂ reflecting the higher values in this period. As for (Fig. 9b), the retreat of the 10 year return period runout altitude is interrupted since ∼2006, due to the slightly higher values for since ∼2006 (Fig. 3) and the slightly lower runout altitudes since ∼2000 (Fig. 8).

For , the minimal altitude possible is attained for annual estimates but also for both low- and intermediate-frequency trends for many winters in the middle of the study period (Fig. 11a). More generally, the low-frequency pattern looks more like that inferred for (Fig. 9b) than that inferred for the mean runout altitude z_t (Fig. 10a) with, for instance, an increase over the second half of the study period higher than the decrease over the first half of the study period. Thus, both M₁ and M₂ give elevations of ∼1260 m at the beginning of the study period, whereas, for the ten last winters of the study period, , which is 20–30 m above . This latter difference is greater than the recent gain in precision in the EPA runout altitude survey. Hence, even if the return period for reaching the valley floor is a partially biased indicator, its recent very important increase corresponds, at least partially, to a significant retreat of large avalanches over the last 30 winters, or at least over the ∼1980/85–2000/05 period if one takes into account the recent inflexion.

In terms of north/south differences, it should be noted that the interannual mean of is higher in the south than in the north. Hence, is attained during more winters in the south for any return period. Nevertheless, for (Fig. 10b and c) and (Fig. 11b and c) and in both regions, the major result is a large increase since a marked change-point in ∼1980. This confirms the general retreat of large avalanches since this time all over the French Alps, but with the change occurring earlier and more dramatically in the north, according to change-point dates for and (Table 4).

In detail, in the northern Alps, before the change-point, mean and low-frequency trends remained almost constant, ∼20–30 m above the interannual mean for , and at the interannual mean for . At intermediate frequency, marked low values are noticeable in the early 1950s due to higher values of . In 1976, both and trends fall very sharply before beginning a fairly steady rise from ∼1980 to the early 2000s, later for because is attained for nearly all winters between 1976 and 1990. At intermediate frequency, the ‘end of large avalanche retreat’ occurs earlier than at the entire Alps scale, i.e. just after 2000, but corresponds to a plateau rather than to a decreasing trend, the nearly continuous decrease in avalanche occurrences (Fig. 6a) being just compensated by avalanches reaching slightly lower runout altitudes since ∼2000 (Fig. 4a).

In the southern Alps, it is the decrease until about 1980/85 that is more continuous, the decreasing trend in (Fig. 4b) being enhanced by the increasing trend in (Fig. 6b). A sharp rise occurs after 1985, but is visible only for (for , is attained for nearly all winters between 1960 and 1990). The rise at low frequency then continues due to the increasing trend in , but less strongly than in the northern Alps because of the concomitant slight increasing trend in . At intermediate frequency, high-magnitude avalanches have begun to reach lower altitudes again only since ∼2006, but, unlike what happens in the northern Alps, this is because of the strong increase in avalanche counts over recent winters (Fig. 6b), which more than compensates the continuous increasing trend in (Fig. 4b).

4.6. Correlation with synthetic climatic covariates

At the entire Alps scale, the major low-frequency pattern in the SAFRAN winter temperature at 2400 m, Temp _t , is a smooth increase of >1 °C between 1980/85 and ∼2000 (Fig. 12d), which characterizes the well-documented and accelerated climate warming in the entire Alpine space over this period. Also noticeable are the nearly constant values (with a very slight decrease) before 1980, and the inflexion through colder winters again since ∼2000. Regarding mean Crocus snow depths at 2400 m, Snow _t (Fig. 12a), there is a sharp increase in the 5 year running mean in 1976, followed by a 10 year period of snowier winters and then a drop around 1990. The low-frequency pattern is flatter, but with a noticeable decreasing trend between ∼1980 and 2000, and a slight increase since then.

Fig. 12.

Synthetic snow and weather covariates. (a–c) Modelled Crocus mean winter snow depth: (a) entire French Alps, (b) northern French Alps and (c) southern French Alps. (d–f) Modelled SAFRAN mean winter temperature: (d) entire French Alps, (e) northern French Alps and (f) southern French Alps. Considered altitude is 2400 m.

In terms of north/south differences, the most remarkable features are the higher interannual mean snow depth and lower interannual mean temperature in the northern Alps which explain why mean and high-magnitude runouts are, in mean, lower in this region (Fig. 12b, c, e and f), and various differences in the Crocus winter snow depth series: a more marked ‘bulge’ of snowier winters around 1980, a higher interannual variability over the first winters of study in the northern Alps, and a higher interannual variability over the most recent winters in the southern Alps. Note also that the 1980–2000 low-frequency snow-cover decrease is more marked in the northern Alps, whereas a net increase in the low-frequency snow-cover pattern is visible since 1999 in the southern Alps, mostly because of high values during the last two winters of the study period (Fig. 12b and c). Except for a higher interannual variability in the northern Alps over recent winters, the patterns in SAFRAN temperatures are quite similar in the two regions, highlighting the larger spatial scale of temperature changes compared to changes in precipitation and snow cover (Fig. 12e and f).

Reference Eckert, Baya and DeschatresEckert and others (2010b) showed that runout altitude fluctuations at the entire French scale are well correlated with temperature and snow-depth measurements and other climate proxies at mid- and high altitude. This is even truer for the two synthetic climatic series considered here (Tables 5 and 6). For , correlations are positive with Snow _t and negative with Temp _t , whereas for and they are negative with Snow _t and positive with Temp _t . This indicates more avalanches and lower mean and high-magnitude runouts during snowier and colder winters. For , correlations with temperatures and snow depths are higher for fluctuations and annual values than for trends, while for , correlations are generally enhanced for low-frequency trends. As a synthetic variable combining and , correlations remain high and significant at all frequencies for , suggesting a mixed low- and high-frequency climate control of high-magnitude avalanches by temperature and snow depth.

Table 5.

Empirical correlation with mean snow depth at 2400 m at different frequencies. Considered subperiod is 1958–2008. For avalanche variables, annual estimates (M₀), low-frequency trend (M₁), intermediate-frequency trend (M₂) and fluctuations (M₀-M₁) are considered. For the snow depth data, annual values, 15 year running means (low frequency), 5 year running means (intermediate frequency) and fluctuations (annual–low frequency) are considered. Emboldened values are nonzero at the 5% significance level

Table 6.

Empirical correlation with mean winter temperature at 2400 m at different frequencies. For avalanche variables, annual estimates (M₀), low-frequency trend (M₁), intermediate-frequency trend (M₂) and fluctuations (M₀–M₁) are considered. For the temperature data, annual values, 15 year running means (low frequency), 5 year running means (intermediate frequency) and fluctuations (annual–low frequency) are considered. Emboldened values are nonzero at the 5% significance level

Although there are issues regarding the quality and consistency of the EPA protocol, these results, when taken together, constitute convincing evidence for a climatic explanation of the temporal fluctuations of our different avalanche indices. First, the 10 year period of higher avalanche numbers and lower runouts around 1980 is consistent with snowier and slightly colder winters, especially in the northern Alps. Second, the decreasing trend in avalanche numbers, coupled with the increasing trend in avalanche mean and high-magnitude runout elevations between 1980/85 and 2000/05, corresponds well to the period of marked warming, and to slightly decreasing snow covers. Third, the very recent ‘end of large avalanche retreat’ corresponds well to winter temperatures again becoming slightly lower in both regions, whereas the two last winters of high avalanche occurrences in the southern Alps are related to important snow-cover excesses.

At the entire Alps scale, Snow _t and Temp _t have roughly the same explanatory power for the different avalanche indices. However, snow depth seems to have a stronger influence in the northern Alps (Table 5), whereas correlations are better with temperature in the southern Alps (Table 6). This is particularly true for low- and intermediate-frequency trends. Hence, the dramatic ∼1977 change-point in occurrences and runout altitudes in the northern Alps corresponds closely to the shift in winter snow depths in this region, whereas the later and more gradual ∼1979–84 change-point in the southern Alps is similar to what is observed for temperatures.

Discussing in detail the impact of climate change on the physical processes that control avalanche release and flow (snow accumulation, snowpack transformation, snow transformation during flow, etc.) is beyond the scope of this paper. However, simple physical explanations of our results can be postulated. Avalanche release is climatically controlled through the amount, stratigraphy, humidity, etc., of the available snow. This forms the basis of existing forecasting methods and models (e.g. Gassner and Brabec, 2002) on short daily scales, and makes intuitive sense on longer timescales. Similarly, runout distance is generally positively correlated with the volume of flowing material (Reference Dade and HuppertDade and Huppert, 1998; Reference Bartelt, Bühler, Buser, Christen and MeierBartelt and others, 2012), which causes a direct control of the runout process from the amount of snow precipitation and an indirect control by temperature through higher snowmelt and/or a higher proportion of rain-on-snow events. Furthermore, snow quality (density, humidity, grain size, etc.) also influences friction during the flow. For instance, higher temperatures lead to higher basal friction close to rest (Reference Casassa, Narita and MaenoCasassa and others, 1989). Thus, there is greater drag when wet snow is involved, providing another connection between winter temperature and runout of high-magnitude events.

4.7. Links to flow regime type

Trends in the proportion of avalanches with a powder part, PSA _t , were analysed for the merged French Alps and Pyrenees data between 1946 and 2006 in Reference Eckert, Baya and DeschatresEckert and others (2010b), and showed a significant decrease. Focusing on the Alps and adding the last winters into the analysis leads to even more significant results because of the recent low values of PSA _t . At the scale of the entire French Alps, Figure 13a shows a negative linear trend of −0.3% winter⁻¹, from 25% in 1973 to around 13% in 2009, mainly because of the strong trend in the southern Alps shown in Figure 13c (−0.4% winter⁻¹, from 23% in 1973 to 7% in 2009). This trend also exists in the northern Alps, but is not nonzero at the 5% significance level over the full 1973–2009 period. However, it is significant over the 1977–2009 period, i.e. starting at the preferred date of change previously highlighted in this region (Fig. 13b).

Fig. 13.

Proportion of powder-snow avalanches: (a) entire French Alps, (b) northern French Alps and (c) southern French Alps. The linear fit is made on the full 1973–2009 subperiod (i.e. without considering a possible change-point).

At the scale of the entire French Alps, because of the overall decreasing trend, high annual proportions are concentrated around 1980, with three values above 30%, but a sharp peak corresponding to the 1997 and 1998 winters is detected by the 5 year running mean filter. At the northern Alps scale, there are even four annual proportions above 30% around 1980, and the 1997–98 peak is >40%. In the southern Alps, things are quite different, with the strong overall decreasing trend in PSAt mostly driven by the (very) low annual proportions recorded since ∼2000, but low values around 1980 and a long period of rather high values between 1985 and 1998 are also noticeable.

These patterns are in agreement with our avalanche indices and their climatic controls. Indeed, the development of powder clouds during avalanche flow generally requires harsh winter conditions with a ready supply of cold dry snow, and long runouts often correspond to powder-snow or mixed avalanches. It is therefore logical to have positive correlations between PSA _t and , and negative correlations between PSA _t and both and , indicating more avalanches with a powder part during winters with more avalanches, and with lower mean and high-magnitude runouts. This also explains well the concentration of winters with high proportions of powder-snow avalanches around 1980, and the predominant decreasing trend in PSA _t concomitant with the warming, with especially low values since 2000 corresponding to winters where the probability of reaching low runouts is lowest over the study period. Finally, north/south differences also agree with regional differences in the evolution of the main climatic drivers. In the northern Alps, the decreasing pattern in PSA _t is less marked than at the entire Alps scale, interrupted in 1997/98, starting abruptly in ∼1977 and hence closely related to the shift in Crocus snow depth at this time. This is similar to what is highlighted by model M₁ for runout altitudes, and, although to a lesser extent, for avalanche occurrences in this region (Figs 4a and 6a). By contrast, the stronger decreasing trend in the southern Alps is better related to the low-frequency pattern in SAFRAN temperatures, but with a later and smoother inflexion just before the mid-1980s, similar to what is highlighted by M₁ for the avalanche data in this region (Figs 4b and 6b).

At all frequencies, correlations are slightly weaker between and PSA _t than between PSA _t and and/or (Table 7). This is presumably because flow regime really controls the runout process whereas there is only an indirect relation between the number of events and the flow regime through the amount, nature and repartition of snow. Furthermore, correlations are very strong for the low-frequency pattern, and remain rather strong for the intermediate-frequency pattern, at least with and , whereas they are non-significant between fluctuations, and significant for the annual values only for the southern Alps and for and . Hence, the flow regime, avalanche occurrences and runout altitude indices may be linked by a long-term, joint climate control rather than by a common response to year-to-year variability, possibly explaining why the recent movement towards an increase in and a decrease in is not visible in Figure 13. Other possible explanations are the lower quality of the PSA _t data and that we are examining proportions only and not the number of events.

Table 7.

Empirical correlation with the annual proportion of avalanches with a powder part at different frequencies. Considered subperiod is 1973–2009. For avalanche occurrences and runout altitudes, annual estimates (M₀), low-frequency trend (M₁), intermediate-frequency trend (M₂) and fluctuations (M₀ – M₁) are considered. For the proportion of powder-snow avalanches, annual values, 15 year running means (low frequency), 5 year running means (intermediate frequency) and fluctuations (annual–low frequency) are considered. Emboldened values are nonzero at the 5% significance level