Climate warming enhances snow avalanche risk in the Western Himalayas.

Ongoing climate warming has been demonstrated to impact the cryosphere in the Indian Himalayas, with substantial consequences for the risk of disasters, human well-being, and terrestrial ecosystems. Here, we present evidence that the warming observed in recent decades has been accompanied by increased snow avalanche frequency in the Western Indian Himalayas. Using dendrogeomorphic techniques, we reconstruct the longest time series (150 y) of the occurrence and runout distances of snow avalanches that is currently available for the Himalayas. We apply a generalized linear autoregressive moving average model to demonstrate linkages between climate warming and the observed increase in the incidence of snow avalanches. Warming air temperatures in winter and early spring have indeed favored the wetting of snow and the formation of wet snow avalanches, which are now able to reach down to subalpine slopes, where they have high potential to cause damage. These findings contradict the intuitive notion that warming results in less snow, and thus lower avalanche activity, and have major implications for the Western Himalayan region, an area where human pressure is constantly increasing. Specifically, increasing traffic on a steadily expanding road network is calling for an immediate design of risk mitigation strategies and disaster risk policies to enhance climate change adaption in the wider study region.

Climate warming at the Third Pole and the ongoing melting of snow and ice (1) are anticipated to change the magnitude and frequency of cryospheric hazards in ways that will permanently change high mountain landscapes and associated socioeconomic systems (2). In the Himalayas, ice and snow avalanches have left strong environmental footprints, but both past and present activity, and the associated risks, remain poorly documented (3, 4). In the Western Himalayan region, snow avalanches often block critical transport corridors (5, 6) and cause human, property, livestock, and/or infrastructure losses (7, 8), which, in turn, disrupt the status quo of communities and may put the future welfare of people living in mountain valleys at risk.

Snow avalanche activity is controlled by the variability of snow and weather conditions and their interactions with topography. Avalanche release is often spontaneous, but it is sometimes caused by an external loading (e.g., by humans or animals) (9). In addition, avalanche activity is intrinsically variable in terms of the frequency, magnitude, seasonality, and typology of events, which are all dependent upon snowpack characteristics, such as snow type, thickness, and stratigraphy (10). For instance, powder snow avalanches tend to occur after intense snow precipitation during cold winter conditions, whereas wet and dense flows often coincide with warm spells, typically toward the end of the winter and early spring (11, 12). As a consequence, changing climatic conditions may modify avalanche activity. Land cover changes, such as afforestation and deforestation, are also likely to play a role (13). In recent decades, several studies have demonstrated decreasing trends in snowfall and snow cover duration in low- to middle-elevation regions of mountainous areas in both North America and Europe (14–17). These changes have led to an upslope retreat of the areas affected by large avalanches (18) and to the occurrence of wet snow avalanches at places where avalanches were mostly dry before (19). At higher altitudes, especially in the world’s highest regions, relationships between climate and snow avalanching remain unclear due to the general lack of long-term observations. The ongoing warming and related changes in precipitation, in combination with the threshold effect of the freezing level, will add further complexity to changes in snowpack and snow avalanching at critical elevations (20–22). One can thus assume that more intense heavy snowfall (22) and/or increasing winter temperature variability (23) will likely lead to an increase in high-altitude snow avalanche activity (24, 25).

Here, we reconstruct a snow avalanche history for the Indian Himalayas and investigate whether ongoing climate warming has had an impact on the frequency and magnitude of snow avalanching in the recent past. We hypothesize that the observed positive trends in air temperature during winter and spring impact snow avalanche occurrence on subalpine slopes in the Himalayas. The study focuses on a site in the Western Indian Himalayas (Fig. 1 and details are provided in ), where damage has been reported repeatedly along the only reliable transport corridor linking the lowlands of the Indian subcontinent to the remote mountain regions of Lahaul-Spiti and, ultimately, Leh (Ladakh). Due to the key strategic role of the Rothang Pass and its closure in winter, the Government of India is currently finalizing construction of the longest tunnel (26) in the Hindukush-Karakoram-Himalayan (HKH) region to allow year-round access to Ladakh.

Fig. 1.

Avalanche slope in the Western Himalayas used for the reconstruction of changes in avalanche frequency. Red dots indicate the locations of sampled trees. Potential release areas are indicated with semitransparent white surfaces and have been detected using the approach suggested by Bühler et al. (26). The access road to the new Rothang tunnel crosses the lower part of the slope.

To overcome the striking lack of data and to put recent snow avalanching into perspective, we used growth-ring records of trees on a forested avalanche slope to develop the longest avalanche reconstruction available for the HKH region. Tree-ring–based snow avalanche reconstructions have been shown to be a highly accurate means to retrospectively date the occurrence of snow avalanches in time and space and with annual-to-seasonal precision (27, 28). The use of tree rings has a long tradition in snow avalanche reconstructions in North America and Europe (29, 30), but information is not yet available for the Indian Himalayan range (3). To analyze snow avalanche/climate linkages, we use data from climate reanalyses and an innovative modeling strategy that takes into account year-to-year memory, namely, a generalized linear autoregressive moving average model [GLARMA (31)] ().

Long-Term Snow Avalanche Reconstruction in the Western Himalayas

Tree-ring analyses of 144 trees growing on the avalanche slope allowed detection of 521 growth anomalies induced by past snow avalanches (Fig. 2 and ) and 38 y since 1855 CE in which snow avalanches occurred. We find a long-term average occurrence ratio of 0.24 avalanches per year (Fig. 3). Our snow avalanche reconstruction reveals evident contrasts in activity over time; trend analyses indeed point to phases with virtually no activity between the 1940s and 1960s. By contrast, we observe very high activity between 1970 and 1977 and from 1989 to 2003 (with occurrence rates >0.875 y−1). Spatial analyses of the areas affected by snow avalanches point to large spatial footprints during the 1970s and 1990s, and especially in the years 1991, 1974, and 2006. By contrast, snow avalanches left smaller spatial footprints before the 1970s. The minimum width of the area affected by snow avalanches in any given year is 132 m, which is much larger than the minimum distance between analyzed trees ().

Fig. 2.

(A) General view of the investigated avalanche slope. (B) Trees disturbed by snow avalanches with characteristic avalanche scars on their stem surface. (C) Sampling of avalanche-damaged trees with an increment borer.

Fig. 3.

Snow avalanche reconstruction based on tree-ring records from 144 broadleaved and conifer trees. (Upper) Chronology of 38 reconstructed snow avalanches. N, no; Y, yes. (Lower) W index (blue columns), the number of GDs (gray columns) observed each year, the sample depth (i.e., number of trees available for analysis in any given year), and the related W and number of GD thresholds (blue and red lines, respectively) used to distinguish avalanche signal from noise.

Hence, the tree-ring–based snow avalanche reconstruction reveals quite substantial changes in process activity since the 1970s, both in terms of the area affected and the frequency of snow avalanching. The robustness of these patterns of change is clearly supported by the increase in the number of intermediate and strong growth anomalies in the tree-ring records since the 1970s (). These changes in reconstructed snow avalanching can be explained neither by changes in sample size (i.e., trees available for analysis in any given year) nor by the distribution of trees sampled on the slope (32, 33). In fact, the snow avalanche reconstruction is based on a substantial number of trees covering the entire 20th century (>85 disturbed trees) and on a homogeneous ratio between conifer and broadleaved trees (∼1.26; ). A possible masking of old scars () after several decades (34, 35) cannot explain the detected trends in activity either, as scars occurred in around 60% of all cases since the 1870s, indicating that our sampling strategy minimized the risk of missing old, overgrown scars. Thus, and even if we cannot rule out the possibility that some avalanches may have passed at the site unnoticed (36), especially during the first half of the 20th century, the reconstructed change in snow avalanching is indeed reflecting changes in process behavior.

More generally, the recent intense snow avalanche activity reconstructed from tree rings is fully supported by direct observations of snow avalanching at adjacent sites, even if these records are too diffuse in space and time to allow sound time series analysis (5). Indeed, during the 1990s, avalanches were recorded in most years in the wider study region (Kullu district), and during the 2000s, more than 15 events were noted, with the largest ones on record in March 2002, March 2003, January 2006, and January 2008 (37). The high snow avalanche frequency reconstructed over the past four decades also matches with 21st century avalanche records from the nearby Lahaul (∼40 km) (38), Chamba, and Kinnaur regions (37). By contrast, to avalanche time series reconstructed from tree-ring evidence, such direct observations provide additional information regarding the timing and nature of recent events. Notably, many of the reported avalanches occurred in late winter and/or early spring, and wet snow avalanche deposits predominated in recent observations (39). The above core of evidence supports our findings that a substantial shift in snow avalanche activity occurred over recent decades, not only on the studied slope but also on similar slopes in the surrounding region. Recent climate warming is therefore the most plausible explanation for the drastic increase in process activity at the bottom of slopes, where events are nowadays being felt and hence recorded.

Warming Temperatures Increase Probability of Snow Avalanching

A GLARMA model was used to investigate relationships between monthly climate covariates (40) and the occurrence of snow avalanches since the beginning of the 20th century. Unlike previous work, which relied mostly on logistic regression of classification trees (41–43), we applied a model that explicitly includes an autoregressive component to explain the nonlinear nature of predisposing factors of snow avalanches, such that an avalanche occurrence in 1 y can influence the likelihood of subsequent avalanches via, for example, changes in vegetation patterns (44, 45). Principal component analysis (PCA; details provided in ) was applied to the set of available standardized climate variables and supports a predominant contribution of cumulated precipitation during February and March (PCA axis1 = 24.4%; ). The second component of the PCA is represented by warming air temperatures between December and March (PCA axis2 = 16%). The third component of the PCA describes the variability of air temperatures in January (PCA axis3 = 14.8%). It is noteworthy that if combined, the three principal axes of the PCA are able to explain 55.2% of the total variance of climate variables (relative contributions are provided in ).

A closer look at the GLARMA model outcomes suggests that snow avalanche probabilities are highest if warmer temperatures persist during these months. Thus, the only significant influence on the onset probability of snow avalanche is given by PCA axis2, as shown by the z-ratio tests on parameter estimates (Table 1). The Akaike information criterion (AIC; ) supports the statistical model based exclusively on this first-order autoregressive parameter (AIC = 106.1) against two- or even three-order parameter models (AIC = 116.2 and 115.4, respectively). This finding is consistent with existing significant influences from previous-year conditions on present-year probabilities of snow avalanche occurrence. Interpretation of a GLARMA process is also supported by both the likelihood-ratio test (LRT) and Wald test (). The fact that this model only points to a first-order term is explained by the intense snow avalanche activity at the study site, which could thus be strong enough to erase longer time-scale memory effects.

Table 1.

Model terms for the period 1900–2010

Model terms	Estimate	SE	Z-ratio	Pr(>|z|)
(Intercept)	−1.46975	0.75416	−1.949	0.0513
PCA1	0.28107	0.26049	1.079	0.2806
PCA2	0.64892	0.29461	2.203	0.0276*
PCA3	−0.02726	0.23955	−0.114	0.9094
ɸ1	0.93906	0.04106	22.87	<2e-16***

Null deviance is 137.61 on 109° df, residual deviance is 15.57 on 105° df, and the AIC is 106.1154. Pr(>|z|) is the probability of finding the observed Z-ratio in the normal distribution of Z with a critical point of |z|. *P = 0.05; ***P = 0.001.

Interestingly, model predictions point to an increased probability of snow avalanches toward the second half of the 20th century, which coincides with warming trends observed between December and March, as well as with an increase in accumulated precipitation totals in January and February (Fig. 4). Both of these tendencies clearly increase the probability of snow avalanching. Nevertheless, any increase in snow avalanche probability above a given threshold of 0.5 will also require notable departures from mean temperatures and cumulated precipitation totals, as shown in Fig. 4. It is worth noting that all of the assumptions of the GLARMA model are met, as shown by the graphical analyses of residuals (). In addition, the fitted GLARMA model shows a high goodness of fit with a deviance-based R2 of 0.88. These outcomes are also consistent with the performance of the GLARMA model when applied only to the period 1950–2010 (), therefore underlining the robustness of conclusions, even if, hypothetically, our sampling strategy would have added uncertainties or biases to the reconstruction during the first half of the 20th century.

Fig. 4.

(A) Snow avalanche predictions based on the retained generalized linear mixed model (gray line) vs. the retained GLARMA model (colored surfaces) for the 20th and early 21st centuries using the same covariates. Circles represent information gathered from tree-ring records. Red dots indicate years with snow avalanche activity according to our reconstruction, and blue dots mark years without snow avalanche activity according to our reconstruction. Similarly, the red-colored areas given by the GLARMA represent avalanche probabilities exceeding 0.5. (B) Predictions for the probability of snow avalanches to occur as a function of PCA2. Principal climatic variables in terms of their individual contribution with the PCA2 are included below the plot. Min, minimum; Prec., precipitation; Temp., temperature.

Finally, our results confirm that ongoing climate warming has been responsible for a shift in snow avalanching on subalpine slopes in the Western Indian Himalayas in recent decades (9, 10, 46). According to direct observations, this influence concerns wet snow avalanche activity in particular. Process-based studies also highlight that increasing air temperature would lead to a rise in liquid water content of the snowpack and an increase in its shear deformation rate, which, in turn, favor increased strain at the interface of slab and/or weak layers, and thus, ultimately, the release of wet snow avalanches (12). The increased avalanching in the Western Indian Himalayas is thus likely to result from the preservation of a sufficiently thick snow cover in high-elevation avalanche release areas combined with more frequent crossing of the melting point of snow, which, in turn, results in more frequent, mostly wet, snow avalanches. The transformation of dry snow packs into wet snow packs is decisive for the release of snow avalanches in the region, as the reconstructed increase in process activity occurred even in the absence of increased snowfall. This conclusion for the Western Himalayas is in line with observed changes in the typology and timing of avalanche activity in the European Alps (19), as well as with existing, yet still scarce, data on any related elevation-dependent change in these patterns (24, 25).

Also, above a certain threshold, the increase in the liquid water content of snow in motion will tend to reduce friction, increasing avalanche runout distances (47), while conserving high-impact pressures even close to the point of rest (48), thereby smoothing the avalanche path (9, 12). The documented increase in the extent of avalanches and the severe damage observed in trees () further support the fact that the current late-winter and/or early-spring warming has enhanced the occurrence of wet snow avalanches that are able to reach the bottom of subalpine slopes in the Western Indian Himalayas. Given the potentially dramatic consequences of such high-magnitude events with larger extent and higher pressures, we conclude that ongoing climate warming is currently exacerbating avalanche risk in the study region. This finding therefore calls for well-coordinated actions aimed at improving risk management of snow avalanches and disaster risk policies to enhance climate change adaption in the wider study region.

Material and Methods

Sampling was conducted on a characteristic east-facing snow avalanche slope in the Western Himalayas, located between the villages of Solang and Dhundi, Kullu district, Himachal Pradesh, India (longitude/latitude: 32.33° N/77.14° E; Fig. 1). The slope includes different, rather unconfined, avalanche paths with the highest release areas (26) at 4,200 m above sea level (masl) and lowest runout locations at 2,600 masl. The avalanche path and runout zones are covered by a loose-forest stand growing on a rather homogeneous slope with an average angle of ∼35°. The transport zone is partially vegetated and has a mixed forest of maple and spruce, with an average slope of ∼20° and a maximum slope width of ∼1,500 m. The runout zone ends at the contact of the slope with the Beas River at ∼2,600 masl. Precipitation in winter is controlled primarily by western disturbances.

Due to the lack of reliable long-term climate data in the area, we used gridded temperature and precipitation records from the University of East Anglia, Climatic Research Unit Climatic Research Unit TS3.2 dataset (49) and the Global Precipitation Climatology Centre (gpcc.dwd.de). The selection of these particular gridded datasets was based on (i) their temporal coverage (i.e., 1900–2010) and (ii) significant positive correlations with available records of the Indian Meteorological Department (IMD) at the monthly scale. Reliability of datasets was validated with meteorological records from the Manali station and with NASA’s Tropical Rainfall Measuring Mission (https://trmm.gsfc.nasa.gov/) data (). We also used the IMD’s long-term records from the Shimla station to detect trends in different time windows ().

Dating of Snow Avalanche Events.

We used snow avalanche damage in tree-ring records to date past snow avalanches in this study with annual resolution (36). In the field, trees disturbed by past snow avalanche activity were screened along the avalanche path for (i) scars induced by material transported in snow avalanches, (ii) trees tilted by avalanche pressure, (iii) trees decapitated by avalanches, as well as (iv) surviving neighboring trees growing along the snow avalanche track (28, 32, 50). Disturbed trees record growth disturbances (GDs) in their tree-ring series in the form of (i) injuries and callus tissues, (ii) tangential rows of traumatic resin ducts, (iii) reaction wood, (iv) abrupt growth decreases, or (v) abrupt growth releases. In addition to the impacted trees, we sampled undisturbed trees from a nearby slope to build a reference chronology representing “normal” growth at the site. The reference series was used for visual cross-dating of disturbed samples and for the detection of avalanche-induced growth anomalies. Trees were sampled with increment borers and a handsaw. We recorded additional information, such as the typology of GDs, geographical location, and graphical information. In the laboratory, samples were analyzed following standard protocols (29, 36), which involved (i) sample preparation (polishing and sanding), (ii) tree-ring width measurement, (iii) cross-dating of series using so-called “pointer years” [i.e., years with remarkable growth responses at the stand level (51)] from the reference chronology, as well as (iv) identification of growth anomalies in the tree-ring series of “avalanche trees.”

Past snow avalanches were reconstructed on the basis of (i) the number of GDs observed for any particular year and (ii) the relative number of trees showing a GD (compared with all trees available for analysis in that year) and the intensity of GDs in that year [i.e., weighted index (36)]. To take account of the increasing number of samples available for analysis and the intensity of GDs, we computed the weighted, W, index as follows (52):

where the sum of trees with injuries (T) was multiplied by a factor 7, the sum of trees with a strong GD (T) was multiplied by a factor 5, and the sum of trees with an intermediate GD (T) was multiplied by a factor 3. R represents the number of trees showing GDs as a response to a snow avalanche in year t, and A represents the total number of trees alive in year t. The threshold used for the GDs and W index was defined dynamically according to the existing literature (32). A given year was accepted as an avalanche event year when both the GD and W indices passed their threshold. The location of disturbed trees was used to estimate the relative runout distance and magnitude of each avalanche. To this end, the width of the snow avalanche is given by the maximum orthogonal distance between affected trees. The relative runout during each snow avalanche year is given by the maximum orthogonal distance between the lowest tree showing avalanche damage and a reference line located at the level of the highest tree affected by the snow avalanche. The relative total affected area was defined by the bivariate assessment of the width and relative runout. Finally, to differentiate between large, moderate, and low size events, we performed a cluster analyses on the relative total affected area.

Statistical Analyses.

We identified trends in the climate data, both temperature and precipitation, using the Mann–Kendall test applied to moving arbitrary time windows of 10 y (53). Then, we looked at the statistical relationship between climate and snow avalanche events. For this analysis, we considered temperature and precipitation variables at a monthly scale, specifically, maximum and minimum average monthly temperature, as well as accumulated precipitation and anomalies for December, January, February, and March. The following equation was applied to estimate monthly anomalies in the series:

where µref and σref are the mean and SD for the reference period and Xij is the value for each specific month, respectively. The reference period considered here was 1980–2010.

For the statistical modeling of the avalanche/climate linkage, we first reduced redundant dimensionality in the set of climate variables by using PCA with all climatic variables centered and standardized by the norm (54). Once the climatic information was conveniently summarized, a GLARMA was fitted to analyze the relationship between the binary response variable (i.e., occurrence of snow avalanche events at an annual time scale from 1900 to 2010) and the first three axes of the PCA analysis as orthogonal predictor variables. GLARMA models are a class of observation-driven state-space models where the state vector consists of a linear regression component plus an observation-driven component consisting of an autoregressive-moving average filter of past predictive residuals. The model estimates GLARMA parameters using a Newton–Raphson iterative process. The linear prediction of the response is as follows:

where Z is the infinite moving average assessed using the autoregressive moving average recursions:

where ϕ is the autoregressive parameter and p is the order of the autoregressive term. To validate the inclusion of the autoregressive parameters in the model, we tested the null hypothesis that autoregressive parameters do not have any relevant contribution. Two independent tests were applied to test the aforementioned null hypothesis: (i) a likelihood ratio test and (ii) a Wald test. P values equal or lower than 0.05 suggest the rejection of the null hypothesis; thus, a GLARMA model is supported. We used the AIC with a correction for small sample sizes (AICc) to select the appropriate autoregressive order in the GLARMA model. AICc criteria combine a measure of goodness of fit with a penalty term based on the number of parameters (k) used in the model. The influence of fixed parameters (e.g., parameters associated with PCA components as explanative variables) in the linear model was tested using a z-ratio with a threshold probability of 0.05 for rejection of the null hypothesis that a given parameter is equal to zero. GLARMA models were conducted with package “glarma” in the R environment (31).