Temporal trends in avalanche activity in the French Alps and subregions: from occurrences and runout altitudes to unsteady return periods

Abstract

We present an analysis of temporal trends in ∼55 000 avalanches recorded between 1946 and 2010 in the French Alps and two north/south subregions. First, Bayesian hierarchical modelling is used to isolate low-, intermediate- and high-frequency trends in the mean avalanche occurrence and runout altitude per year/winter. Variables are then combined to investigate their correlation and the recent evolution of large avalanches. Comparisons are also made to climatic and flow regime covariates. The results are important for risk assessment, and the development of new high-altitude climate proxies. At the entire French Alps scale, a major change-point exists in ∼1978 at the heart of a 10 year period of high occurrences and low runout altitudes corresponding to colder and snowier winters. The differences between this change-point and the beginning/end of the study period are 0.1 avalanche occurrences per winter and per path and 55 m in runout altitude. Trends before/after are well correlated, leading to enhanced minimal altitudes for large avalanches at this time. A marked upslope retreat (80 m for the 10 year return period runout altitude) accompanied by a 12% decrease in the proportion of powder snow avalanches has occurred since then, interrupted from about 2000. The snow-depth and temperature control on these patterns seems significant (R = 0.4–0.6), but is stronger at high frequencies for occurrences, and at lower frequencies for runout altitudes. Occurrences between the northern and southern French Alps are partially coupled (R∼0.4, higher at low frequencies). In the north, the main change-point was an earlier shift in ∼1977, and winter snow depth seems to be the main control parameter. In the south, the main change-point occurred later, ∼1979–84, was more gradual, and trends are more strongly correlated with winter temperature.

1. Introduction

Indirect avalanche data from dendrochronology (Jomelli and Pech, 2004) and lichenometry (McCarroll and others, 1995) indicate that major avalanches of the type that occurred during the Little Ice Age have not been encountered in recent decades. Models of snowpack evolution following climate-change scenarios also suggest that changes in triggering mechanisms are already in progress (Martin and others, 2001), and that this trend may persist during the 21 st century (Lazar and Williams, 2008), especially at low and mid-altitudes (López-Moreno and others, 2009). Hence, for hazard mitigation, the assumption of stationarity of high-magnitude avalanches, nearly always made when deriving reference scenarios from a sample of past observations (e.g. Keylock and others, 1999), may be questionable.

The problem of assessing temporal trends in avalanche data has received relatively little attention in the literature. Indeed, past work has tried to correlate avalanche activity to climatic factors, rather than to analyse avalanche time series directly (Keylock, 2003; García-Sellés and others, 2010), primarily because most available avalanche data series are short, incomplete and inhomogeneous. In addition, while possible changes in avalanche activity are likely to be related to climate fluctuations, historical records are also affected by the construction of countermeasures. This makes standard statistical methodologies for trend detection such as stationarity tests (e.g. Burn and Hag Elnur, 2002) hard to implement, precluding firm conclusions despite increasing knowledge regarding recent changes of mountain climate and snow cover (e.g. Beniston, 1997, 2003; Falarz, 2004; Valt and Cianfarra, 2010). For example, Laternser and Scheneebeli (2002) found no changes in avalanche activity over the 1950–2000 period in Switzerland, and Schneebeli and others (1997) found no modifications in the number of catastrophic avalanches around Davos, Switzerland, during the 20th century.

Recently, Eckert and others (2010a,b) introduced a model-based approach for extracting the predominant temporal patterns common to a set of local avalanche series within a hierarchical Bayesian framework. The idea is that averaging the record over a large number of paths should be relatively free from local artefacts and may therefore be more confidently linked to regional forcing such as climate change than a single series. Furthermore, with regard to more empirical approaches, Bayesian hierarchical modelling permits refined underlying trends and significant patterns such as change-points to be extracted and studied, with the different sources of uncertainty treated rigorously (e.g. taking into account missing values and the uncertainty regarding annual estimates when inferring the temporal patterns of interest). Application to avalanche occurrences and runout altitudes from the exceptionally detailed French avalanche chronicle has given promising first results. For instance, Eckert and others (2010a) implemented different autoregressive and shifting level models to highlight abrupt changes in avalanche occurrences in the northern French Alps over the 1946–2005 period, while Eckert and others (2010b) used a single change-point model to highlight a clearer temporal pattern to changes in avalanche runout altitudes well correlated with a few direct and indirect climate data at the scale of the whole French avalanche database, including the Alps and the Pyrenees.

Based on this work, the objectives of this paper are:

to apply the two best-adapted models previously tested to all occurrence and runout altitude data available in the French Alps over the 1946–2010 period. These models are aimed at detecting complementary patterns rather than searching for the one that is optimally adapted to each analysed series. Hence, they quantify the mean evolution as precisely as possible, as well as the presence of underlying trends or change-points in low- and intermediate-frequency signals and in annual fluctuations (e.g. at different timescales). Here we expand their application to fully coherent datasets in terms of spatio-temporal scales, which facilitates inferences regarding the correlation between occurrences and runout altitudes at different frequencies. This allows us to quantify the extent to which winters with many avalanches correspond to winters where average runout altitudes are low;

to combine the occurrence and runout altitude variables to extract major patterns at different frequencies for high-magnitude avalanches in the French Alps. These results, especially for low-frequency trends, are even more crucial for quantifying possible changes in risk and are the first of their kind in the avalanche field;

to quantify the correlations with synthetic climatic covariates resulting from the assimilation of all available snow and weather data, and evaluate how this connects to changes in avalanche flow regimes. This analysis is necessary to investigate whether the changes we highlight in our avalanche data series are driven by climate rather than artefacts;

to consider two subregions so as to infer possible deviations around the mean French alpine effect. This is motivated by different predominant atmospheric patterns in the northern and southern French Alps: mostly Atlantic flows, and mixed Atlantic/Mediterranean flows, respectively. Hence, inference of the predominant climatic drivers in each region becomes possible, and the level of coupling between the two regions can be quantified at the different considered frequencies.

The paper is organized as follows: Section 2 describes the data used. Section 3 briefly presents the models used for extracting major temporal patterns at different frequencies in avalanche occurrences and runout altitudes. How these can be combined to evaluate the recent patterns of behaviour for high-magnitude avalanches is also detailed and the advantages of the chosen methodology for our problem are illustrated. Section 4 presents and discusses the results obtained for the different regions/variables studied, while Section 5 summarizes the main outcomes of the work and points out possible developments.

2. Data

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 (Mougin, 1922). The common use of EPA data is for risk assessment at the path scale (e.g. Ancey and others, 2004; Eckert and others, 2007a, 2009a, 2010c), but links between avalanches and snow and weather covariates (e.g. Jomelli and others, 2007) or with dendrogeomorphological reconstructions (e.g. Corona 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 (Burnet, 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, Eckert and others (2010b) found that similar conclusions were found if the few perturbed paths were included in or excluded from analysis.

Following Eckert and others (2007b, 2010a), 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 Eckert and others (2010b), rather than altitude we use the Runout Altitude Index (RAI):

(1)

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 (McClung 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 Eckert 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 (Durand and others, 1993) and Crocus (Brun and others, 1989, 1992), 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 (Durand and others, 2009a,b). 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 (Castebrunet 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. Extracting Temporal Patterns from Avalanche Time Series

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

(2)

(3)

(4)

(5)

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:

extract annual model estimates from the data;

separate possible systematic trends and change-points from the interannual fluctuations taking into account that the annual common effect is not observed and, hence, not known with certainty, ensuring the significance of possible temporal patterns is not overestimated.

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 Eckert and others (2010a), the annual avalanche counts a_jt are modelled with a non-homogeneous Poisson process inspired by spatial epidemiology (Elliott and others, 2000), with parameters λ_jt , j ∊ [1, M], t ∊ [t ₀, t ₀ + T _obs − 1] summarizing the local annual avalanche activity:

(6)

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:

(7)

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:

(8)

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

(9)

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

(10)

3.3. Extracting the mean annual runout altitude

Following Eckert 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:

(11)

(12)

(13)

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

(14)

which gives

(15)

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, :

(16)

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. Wikle, 2003; Banerjee 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.4.1. Low-frequency linear trend, M1

The low-frequency trend is extracted using model M₁, a single change-point model originally developed in hydrology (Perreault and others, 2000a,b) and successfully applied to various proxies since then (e.g. Eckert and others, 2011). Defining the winter of a possible change-point in behaviour as τ (the change is assumed to occur between τ and τ + 1), then, before and after the change-point, g_t , m_t or logit(p_t ) is separated into a random noise and a linear trend, trendx_t = a_. + b_.t, where x is replaced by the considered variable. The random noise, with variance , models the residual interannual fluctuation. The notation is combined with a subscript to indicate that the parameters can take different values before and after the change-point, for example, (b ₁, b ₂), respectively:

(17)

where N indicates a normal distribution. This model is relatively simple, but, depending on the continuity of trendx_t around τ and on (b ₁, b ₂)and (σ ₁, σ ₂) values, it can capture a monotonic trend and various types of change in mean and variance.

This model has not been previously applied to avalanche counts, but already to runout altitudes in Eckert and others (2010b). Here we just use an additional logit transformation for p_t to facilitate numerical inference. A logit transformation could also be used for m_t ∊ ]0,1 [, but data quantity is sufficiently large, even in the north and south subregions, to constrain its value strongly and avoid any numerical trap during inference.

We search for τ over the subperiod only to prevent ‘boundary effects’ at the beginning and end of the time series. The model imposes the same change-point for m_t and logit(p_t ), whereas a different winter can be selected for g_t since avalanche occurrences and runout altitudes are processed separately. Finally, to obtain , and the associated uncertainty, we substitute and for and and their credibility intervals in Eqns (10) and (16), respectively. Similarly, and the associated uncertainty are obtained by applying the inverse logit transformation to and its credibility interval.

3.4.2. Intermediate-frequency shifting level trend, M2

The intermediate-frequency trend is extracted using a shifting level model (M₂) developed by Salas and Boes (1980) and successfully applied to hydrological series by Fortin and others (2004). It considers any time-series variable x_t to be decomposable into a white-noise component and intermediate-frequency segments of constant trend. The parameter ζ quantifies the annual probability of a change (level shift) in the intermediate-frequency trend. If the Bernoulli variable B_t = 0, it remains identical. If not, a new regime is reached:

(18)

(19)

newmean _t is distributed as a white noise with a variance quantifying the inter-regime variability that is to be compared with the white-noise σ ² component of the x_t terms quantifying the intra-regime fluctuations around the trend:

(20)

Hence, level shifts break the autocorrelation structure, since newmean _t does not depend on newmean _{t −1}. is the interannual mean of the levels. It is set to zero for g_t which is a centred excess/deficit, and estimated for m_t and logit(p_t ). Finally, the balance between the inter-regime variability and σ ² is constrained for the model to be identifiable.

For avalanche counts, this model has been used in Eckert and others (2010a), but never before for runout altitudes. Since the multiple change-points detected are too different for m_t and logit(p_t ), the two series have been modelled independently contrary to what has been done with model M₁. Finally, the trends of interest and the associated uncertainty are obtained, as for model M₁, by applying Eqns (10) and (16) and the inverse logit transformation to the modelled latent variables.

3.4.3. Explained variance

To compare the respective contributions of the low-/intermediate-frequency trend and the random fluctuations, we define, for both models M₁ and M₂, the ratio of explained variance frac.struc:

(21)

3.4.4. M₀: a null model with no trend

A null model, M₀, to which the behaviour of the trend models, M₁ and M₂, may be compared, is formed by modelling the latent variables, g_t , m_t and and logit(p_t ), as white noises with no trends. Hence, for M₀, frac.struc is forced to zero. This model makes use of the shrinkage effect, whereby the temporal structure in M₁ and M₂ constrains the annual estimates (see below).

3.5. Bayesian inference and shrinkage effect

Inference was implemented using Bayesian Markov chain Monte Carlo (MCMC) methods (Brooks, 1998; Gilks 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 (Brooks 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

(22)

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. Results and Discussion

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.

Durand 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 (Rousselot and others, 2010) and also caused widespread damage in Europe (SLF Davos, 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 (Eckert 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 Castebrunet 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 (Voellmy, 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 Eckert 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 Eckert 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 Eckert 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 Eckert (2009) and Eckert 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 Eckert 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 (Bé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).

Eckert 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 (Dade and Huppert, 1998; Bartelt 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 (Casassa 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 Eckert 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

5. Conclusion

We have used an advanced statistical framework to extract temporal patterns from different avalanche data series from the French Alps. For both occurrences and runout altitudes, we separated the hidden temporal pattern common to the different local data series from spatial effect. The spatial effect was explicitly taken into account in the occurrence model, leading to a two-way variance decomposition. It was considered as already separated from the scaled RAI variable, leading to a one-way variance decomposition performed on a non-Gaussian and discontinuous variable. In addition, hierarchical modelling permitted low-, intermediate- and high-frequency signals to be extracted using two distinct time-series models, aimed at detecting complementary patterns, rather than searching for the model that is optimally adapted to each analysed series, allowing a model-based spectrum analysis to be performed.

After checking that the modelling assumptions made were not too strong to produce biased estimates, annual effects and the associated trends were systematically reworked, leading to the mean avalanche number and mean runout altitude per year/winter on a mean path at the whole French Alps scale and for the north/south subregions. This allowed expansion of previous results to datasets fully coherent in terms of spatio-temporal scales, study of the north/south coupling, and of the intervariable correlation in each region. Occurrences and runout altitudes were then combined to evaluate the temporal patterns in (relatively) high-magnitude avalanches rigorously, lowering the bias with regard to the probability of reaching the valley floor, which had been previously adopted. Finally, a correlation study with two synthetic climatic covariates and avalanche flow regime was performed, searching for similarities, so as to determine the main drivers of the highlighted evolutions. Our main results may be summarized as follows:

for occurrences, a partial coupling exists between the north/south regions (R = 0.4), especially at low frequencies (R = 0.71), but it has weakened in recent winters; by contrast, runout altitudes between the north/south regions are nearly decoupled;

the time series for occurrences is less structured than for runout altitudes, making it harder to distinguish low- and intermediate-frequency patterns from the interannual variability. However, for both variables, there is a major change-point ∼1978, with a difference of ∼ 0.1 avalanches per winter and per path in occurrences and ∼55 m in runout altitude between this change-point and the beginning/end of the study period. The change occurred slightly later in the southern Alps, the mean alpine behaviour being the north/south mean weighted by the number of data in the two subregions. The change was also more of a dramatic shift between two distinct levels in 1977 in the northern Alps and a more gradual 1979–84 transition in the southern Alps. In the northern Alps, there are coherent trends after the change-point, and nearly no trend before, except a short period of high activity in the early 1950s. In the southern Alps, significant trends exist before and after the change-point, although their coherence decreases after the change-point;

there is a significant correlation at the annual scale between occurrences and runout altitudes (R∼−0.4), except in the southern Alps, and it enhances temporal patterns in high return period avalanches. This correlation is also enhanced at low frequency (R∼−0.82), becoming significant even in the southern Alps. The concomitant high avalanche occurrences and low runout altitudes lead to minimum high return period runout altitudes around 1980;

a marked upslope retreat of high return period avalanches occurred over the 1980/85–2000/05 period, for instance ∼80 m for the 10 year return period runout altitude, which makes a horizontal runout distance difference as high as ∼450 m on a typical 10° runout slope. However, higher avalanche counts, largely in the southern Alps, since around 2005, and lower runout altitudes, generally in the northern Alps, since around 2000, have led to high return period avalanches again slightly lower in the most recent winters;

there has been a general decrease of ∼12% in the proportion of powder-snow avalanches since 1973, mostly consistent with the evolution of occurrences and mean and high-magnitude runouts;

all these patterns are highly correlated with two synthetic temperature and snow-depth covariates (R = 0.3–0.6), especially in terms of change-point dates, and of low-and intermediate-frequency trends (R up to 0.8), with a greater influence of snow depth in the northern Alps, and temperature in the southern Alps. The climate control seems stronger at high frequencies for avalanche occurrences and at low frequencies for runout altitudes and flow regime. This leads to a mixed control on high return period avalanches, but with a clear impact from warming on large avalanche retreat over 1980/85–2000/05.

Although filtering procedures have been used to exclude major error sources from analysis, the usual limits to avalanche data mean that all results should be interpreted with care. Hence, the discrepancies between the different variables and subregions that have been shown are possibly partially linked to data limitation such as fewer data in the southern Alps during the first part of the study period, and fewer avalanche paths and less homogeneous massifs in this region. However, features such as the ∼1978 change-point and the retreat of large avalanches over the 1980/85–2000/05 period are so clear in all datasets that they reflect reality. The strong and significant correlations with climatic drivers and flow regime proportions provide additional support for this assertion.

Hence, the detected patterns constitute new high-altitude winter climate proxies for the Alps and are potentially relevant for risk assessment considerations. They definitely challenge the assumption of stationarity generally made in long-term forecasting. For instance, the significant link between warming and the upslope retreat of large avalanches over the 1980/85–2000/05 period indicates that the already observed changes may be amplified in the upcoming winters due to ongoing climate change. However, the apparent decreasing exposure of French mountain communities to avalanche risk must be tempered for different reasons: first, because of the ‘end of large avalanche retreat’ observed over the most recent winters, and, more generally, the difficulty of making future predictions on the basis of time trends alone (see below); second, because of the higher variability of over the second half of the study period discussed in Section 4.4; third, because of the significant negative correlation at the winter scale between avalanche occurrences and runout altitudes, indicating that one must still be prepared to face a high number of potentially damaging avalanches simultaneously. This latter point implies that further work is required to undertake an explicit joint approach to the two variables that were here modelled independently.

A limitation of our approach is that it uses time (and space for avalanche occurrences) as covariates, rather than the true physical drivers. The highlighted patterns therefore depend on the time period considered (window effect). For example, in the southern Alps, as discussed in Section 4.1, the low-frequency trend after the change-point changed dramatically if the two last winters were included in the analysis (Fig. 6b). Having the two highest counts at the end of the study period makes future prediction difficult but also shows the complexity of forecasting avalanche time series. To make progress on this problem, in future work it will be necessary to replace time by time-indexed climate covariates, i.e. expand our preliminary use of synthetic climatic covariates by linking our approach to that used by Castebrunet and others (2012) to study the avalanche–climate relation.

Other outstanding questions are whether the concept of a mean temporal signal common to a sample of avalanche paths is appropriate, and what is the best scale to detect it. A partial answer to the first question was given in Eckert and others (2010a), showing that the common temporal signal represents a small but not negligible part of the total variability of avalanche counts, presumably because of similar responses in terms of release/non-release to regional snow and weather forcing (‘avalanche cycles’). For runout altitudes, this quantification remains to be done, since the modelling approach taken here ceases to consider the spatial variability once the valley floor scaling is completed.

Regarding the question of the best scale to detect a common signal, this study has shown that north/south differences exist, leading to regional patterns slightly different from the overall pattern at the entire French Alps scale, and better correlated with the regional evolution of climatic drivers. Hence, even smaller subregions could be studied in further work, with the advantage of presumably even more homogeneous avalanche activity. This may allow more than two significantly divergent temporal patterns to be inferred, but at the cost of a smaller number of data in each region.

Furthermore, we have chosen to fix the definition of the north and south regions based on climatic knowledge. This assumption is reconsidered for avalanche counts by Lavigne and others (2012) who included the classification problem in the modelling, showing distinct temporal patterns in different groups of townships that do not correspond fully to the north/south regions considered in this work. This divergence highlights model sensitivity and a strong altitudinal control on the temporal evolution of avalanche activity. It also suggests that further work is required to better discriminate spatio-temporal and altitudinal effects on avalanche variables before attempting future predictions.

Finally, we analysed only ‘full winter’ annual series, and shorter time periods may be worth considering in further work more linked to short-term forecasting. This would imply adapting our models to quantify the evolution of major avalanche cycles rather than annual behaviour. The two problems are not wholly disconnected since major avalanches generally occur during the strongest storms, which predominantly contribute to the high-activity winters highlighted in this work.