The NH-7 ascends from 400 m at Rishikesh to approximately 2000 m at Joshimath, crossing steep terrain with soil mantled slopes. Mean annual rainfall (1970–2019) varies from 1500 to 2000 mm around Rishikesh to 1000–1200 mm in Joshimath, with 80 %–86 % and 60 %–70 % delivered by the Indian summer monsoon (June to September), respectively (Pai et al., 2014; Swarnkar et al., 2021). Air temperature in Rishikesh is always above freezing and ranges between 4 and 40 °C, whereas the temperature in Joshimath varies between −10 and 20 °C. This climatic gradient is evident in the gradual transition of vegetation. In the lower-lying subtropical region, dense deciduous forests dominate. As elevation increases, these forests give way to temperate broadleaf mixed forests and temperate shrub and grassland communities.

The geological framework of the study area is largely determined by the ongoing Indo–Asian collision that causes crustal thickening and exhumation along large-scale detachment zones and thrust faults. Most of the study area lies within the Lesser Himalaya, between the Main Boundary Thrust in the south and the Main Central Thrust (MCT) in the north, which are both splays of the root detachment, the Main Himalayan Thrust (Fig. 1). As the present-day India–Eurasia convergence is on the order of 36–40 mm yr⁻¹ (e.g., Wang et al., 2001) and approximately half of this is accommodated within the Himalayas (e.g., Lavé and Avouac, 2000), the region is seismically active and bears the potential for large earthquakes (e.g., Kayal et al., 2003; Bollinger et al., 2014; Rajendran et al., 2017).

The highway runs perpendicular to the strike of the orogen and crosses rocks of the Lesser Himalayan Sequence (LHS) and the High Himalayan Crystalline (HHC). These units represent the ancient passive Indian margin and are separated by the MCT. The LHS is mainly composed of sedimentary and low-grade metasedimentary rocks, quartzite, shale, phyllites, and slate with occasional limestone and dolomite, whereas the HHC is characterized by high-grade schist, gneiss and quartzite. These rocks feature a high density of discontinuities like faults, fractures and joints that are important seepage pathways. In locations where the road cuts through weathered rocks and intersects with major faults, the hillslopes are particularly fragile (Prasad and Siddique, 2020).

During the week preceding our survey, large parts of north India experienced a period of strong rainfall. The state of Uttarakhand registered a departure of 1040 % from the long-term (1961–2010) average. The districts of Tehri Garhwal, Pauri Garhwal, Rudraprayag and Chamoli, through which the NH-7 runs, recorded a weekly surplus of 419 %, 679 %, 218 % and 1855 %, respectively (Table S2). Given that the average rainfall amount in October is around 35 mm, approximately 3 times the monthly average rainfall occurred in only 1 week, which is close to the rate that prevails during the monsoon.

The NH-7 is a lifeline for socioeconomic development, which is mainly based on agriculture, trade, tourism, mining and hydropower. Furthermore, the highway is vital for the Indian military to transport personnel and equipment to their outposts along the Indian and Chinese–Tibetan border. Between May and October, more than 1 million pilgrims travel the highway to visit the holy shrines of Badrinath and Kedarnath. Moreover, the highway follows the course of the Ganges and its tributary Alaknanda and thus passes the river confluences known as the five Prayags, namely Devprayag, Rudraprayag, Karnaprayag, Nandaprayag and Vishnuprayag (Fig. 1). In Hinduism, these confluences are considered sacred, attracting pilgrims who bathe in the waters before worshiping the rivers. Finally, there are 10 hydroelectric power plants within 20 km distance to the road, and more projects are planned or under construction, highlighting the road's significance for energy security and economic value.

Traveling to fieldwork in the Chamoli area (Uttarakhand) on 15 October 2022, we recognized numerous partially road-blocking landslides along the road and decided to investigate their occurrence using the following approach.

First, we started to record these landslides using handheld GPS devices, as road workers had already begun to clean the road of the debris, thus rapidly removing evidence of smaller landslides that had detached in close vicinity to the road. We registered landslides along the road on our way from Rishikesh to Joshimath (Fig. 1) on 15 and 16 October, as well as on our way back on 18 October 2022. We only recorded landslides with runouts affecting the road, thus partially or fully blocking it (Fig. 2). We classified the landslides as partially blocking if the emplaced deposits either substantially narrowed the road or debris ran over marginal road markings. Very small landslides with an area of less than ∼10 m² were not considered. Secondly, we cross-checked each recorded landslide with Google Earth using the most recent and historic high-resolution imagery. The images enabled us to identify the landslides that occurred between January and October 2022, encompassing the entire annual monsoon season as well as the heavy rainfall period at the beginning of October 2022. We initially tried to attribute the landslides to the latter period by matching our data with the acquisition dates of the images, but the temporal resolution for many road sections was insufficient. Moreover, historic imagery has become unavailable for India (https://www.thehindu.com/news/national/google-historical-satellite-imagery-disappears-for-india/article66834033.ece, last access: 25 June 2024). We also identified reactivated landslides, where a slip surface and a scar were visible in the imagery prior to January 2022. This was not always straightforward since landslide scars cannot always be clearly distinguished from bare, engineered slopes and road widening. Landslides confidently determined to have occurred before January 2022, as well as those obscured by shadows in Google Earth images, were omitted from the analysis.

Figure 2Examples of partially road-blocking landslides along the highway NH-7. Panels (a) and (c) show locations where hillslopes parallel the foliation and/or bedding. Panels (e) and (f) show a landslide that was reactivated in the observation period (Google, © 2024 Airbus).

Thirdly, we analyzed landslides along the road as a network-attached spatial point process. A spatial point process is a stochastic mechanism that controls the spatial distribution of events or occurrences (Baddeley et al., 2015). As our mapped landslides are events that occur along the road – and only these have been mapped – these events are constrained to lie on a network of lines (Baddeley et al., 2021; Okabe et al., 2006). In our case, the network is rather simple as it consists of only one polyline, but the approach can be equally applied to more complicated network topologies.

We used the open-source MATLAB extension TopoToolbox (Schwanghart and Scherler, 2014) and its numerical object PPS (Schwanghart et al., 2021) to analyze, visualize and model the density of landslides along the road. PPS means point patterns on stream networks, but it can be used with any type of dendritic, undirected network. The numerical approach consists of a fine-pixel approximation, which is controlled by the geometry of the digital elevation model (DEM) from which the data are retrieved. This means that the vector shape of the road is pixelated with the same geometry as the DEM (Schwanghart et al., 2021). We model landslide densities with an inhomogeneous Poisson process model, which is defined by the intensity function λ(u), with u being the horizontal distance along the road. A common parametric model of the intensity is the log-linear model:

where X is a matrix of predictor variables (covariates) that vary along the road, β₀ is an intercept and β is a vector of model parameters. Our approach uses a spatial logistic regression, where the relation between presence probabilities p and explanatory variables X is controlled by the form of the logistic link function logit $p=\ln (p/(\mathrm{1}-p\left)\right)$. As pixel size tends to zero, we have p→0 and $\ln (p/(\mathrm{1}-p\left)\right)\to \ln \left(p\right)$. The limiting Poisson intensity is thus a log-linear function of the covariates (Baddeley et al., 2010). A key property of the model is that the events are independent from each other; i.e., the probability of landslide occurrence is independent from landslides nearby. Spatial dependence of events can occur in different ways, leading to clustering, i.e., points tend to occur close to other points, or inhibition, i.e., there is a characteristic distance or regularity in the spacing between the points. Spatial clustering of landslide events has previously been addressed by Lombardo et al. (2018, 2019) using a Cox process model to emulate the latent spatial effects of unobserved variables, whereas inhibition can be observed, for example, in data where areal non-overlapping phenomena are represented as points (Evans, 2012; Schwanghart et al., 2021). At this stage, we will not include these potential second-order effects on the density of landslides in our model, but we will investigate their possibility using the inhomogeneous K function defined by Ang et al. (2012) once we have modeled first-order effects.

We used the following candidate predictor variables in the log-linear model introduced above. First, we hypothesized that steep hillslope gradients next to the road are more susceptible to mass wasting events. Therefore, we calculated surface gradients based on the Copernicus 30 m DEM (European Space Agency, 2021). We identified areas that lie to the right or left of the road and that are higher than the road itself within a buffer zone of ∼210 m (or 7 pixels). We then selected the nearest DEM pixels and mapped the mean slopes of these nearest pixels to the road network. We chose 210 m (7× the spatial resolution of the DEM) to have sufficient pixels to obtain robust estimates of the slope. These values still vary greatly over short distances along the road and thus we smoothed them using the algorithms (with smoothness penalty parameter K=4) described by Schwanghart and Scherler (2017).

Rainfall patterns exert a strong influence on the spatial occurrence of landslides (e.g., Ching et al., 2011; Joshi and Kumar, 2006; Gariano and Guzzetti, 2016; Ozturk et al., 2021). Thus, we chose five rainfall/precipitation products (see Table 1) to characterize spatial patterns of accumulated rainfall between 1 January and 10 October 2022. Among the five rainfall/precipitation products, IMD1 relies solely on the interpolation of gauge-based measurements from a network of stations provided by the Indian Meteorological Department (Pai et al., 2014). In addition to gauge measurements, IMD2 is enhanced with estimates from the Integrated Multi-satellite Retrievals for Global Precipitation Measurement (IMERG) (Mitra et al., 2009). MSWEP v2 is another merged product that incorporates reanalysis-based, gauge and satellite-derived rainfall estimates (Beck et al., 2017). CHIRPS v2 provides a high-resolution record by combining gauge and satellite data from the NOAA (National Oceanic and Atmospheric Administration; Funk et al., 2015), and the IMERG late run provides gridded multi-satellite precipitation estimates with quasi-Lagrangian time interpolation from NASA (National Aeronautics and Space Administration; Huffman et al., 2019). We projected the geographic coordinates of the rainfall/precipitation grids to UTM zone 44N and bilinearly resampled the data to the resolution of the DEM before extracting the values of each road pixel.

Table 1Overview of the rainfall products.

Download Print Version | Download XLSX

Next, we incorporated information about road widening by creating a logical mask that identifies the widened segments of the highway. Due to the absence of more-detailed data, we rely on the map by Sati et al. (2011), which shows the locations of road widening completed before the 2010 monsoon season. Our envisioned strategy of compiling a more recent database of road widening using historic imagery from Google Earth was rendered impossible after most of the historic images were removed from the public archive in April 2023.

Finally, we obtained a digitized version of the lithological map of Uttarakhand (map scale 1:2 000 000) from the Geological Survey of India (2022). The data contain both stratigraphic and lithological information. Accordingly, the NH-7 crosses 34 different lithologies between Rishikesh and Joshimath. To reduce the number of categories, we summarized and aggregated the lithological information into lithozones with a lesser focus on the stratigraphic context. This aggregation resulted in five lithozones that are dominated by carbonate rocks (1), phyllite and shale (2), quartzite (3), quartzite and igneous rocks (4), and crystalline high-grade metamorphic rocks (5). The reclassification is shown in Table 2. Again, we gridded these data and assigned corresponding lithozones to each road pixel.

Table 2Definition of lithozones.

^* ID refers to the UID given in the original data (Geological Survey of India, 2022).

Download Print Version | Download XLSX

We adopt a Bayesian strategy to infer and identify predictor variables using the function bayesloglinear of the PPS numerical class (Schwanghart et al., 2021). The function provides an interface to bayesreg (Makalic and Schmidt, 2016), a MATLAB toolbox that enables efficient Bayesian modeling and regularization of high-dimensional data. We use a Bayesian lasso estimator with Laplace prior distributions for the regression coefficients (Park and Casella, 2008). This prior results in posterior-mode estimates that are similar to estimates obtained with the lasso penalty (van Erp et al., 2019). Before calculating 1000 samples of the posterior parameter distributions, we calculated 1000 burn-in samples. These are calculated to ensure that the Gibbs sampler converges to the target distribution. In addition, we used a level of thinning of five samples. This means that only every fifth sample was retained in the generated sequence to reduce autocorrelation between the samples and to obtain more independent and representative Bayesian posterior distributions. To this end, we find that 1000 samples are sufficient to characterize the posterior distributions.

Finally, we evaluate the model based on the receiver-operating characteristics (ROC) area-under-the-curve (AUC) approach (e.g., Hanley and McNeil, 1982). We visualize the predictions and inspect and analyze spatial densities obtained from random realizations of the fitted inhomogeneous Poisson process model. Moreover, we evaluate the predictive performance of the model using a 5-fold cross-validation. In this approach, we subdivide the road into 15 km segments, which are then randomly partitioned into five groups. The first group is used as test data while the remaining groups are used to train the model. This step is then repeated for each group and the performance summarized with the AUC. In addition, we test whether additional covariates provide opportunities for further improving the model. The selected attributes include terrain roughness and total curvature as well as land cover derived from the Copernicus Global Land Operations (CGLOPS-1; moderate dynamic land cover 100 m, version 3; Buchhorn et al., 2020), which we reclassify according to Table 3. We investigate these models using a frequentist modeling approach (see the PPS function fitloglinear) and compare models with additional covariates with the Akaike information criterion (AIC, Akaike, 1974).

Table 3Aggregation of land cover classes derived from CGLOPS-1 (Buchhorn et al., 2020). The remaining map codes in the original data were not present along the NH-7.

Download Print Version | Download XLSX

We recorded 309 fully or partially road-blocking landslides along the 247 km long road between Rishikesh and Joshimath, which amounts to an average landslide density of 1.25 landslides per kilometer. A two-sample Kolmogorov–Smirnoff test between the road distance (uniform distribution between start and end of the surveyed road) and the road distances measured at the landslides rejects the null hypothesis, with p≈0 that landslide locations follow a completely random spatial distribution. Visually inspecting the landslide locations using Google Earth reveals that 80 % of the recorded landslides with road blockages occurred after January 2022 (Fig. 1). Of these 250 landslides, 30 % were most likely reactivated because they could not be identified as road blocking before the rainy season (Fig. 2e and f).

The spatial distribution and amounts of accumulated rainfall between 1 January and 10 October differ between the rainfall products (Fig. 3). Since independent measurements based on rain gauges are unavailable, we investigate the performance of the rainfall products to explain the spatial distribution of landslides. The approach uses the AIC to iteratively evaluate five models that include hillslope gradient, lithozones and road widening but each using a different rainfall product. AIC values vary between 1876 and 1913, with CHIRPS v2 having the lowest AIC. CHIRPS v2 correlates positively with landslide density, whereas all other rainfall products show negative or no significant correlation (see Table 1 and Figs. S1–S5). We thus use CHIRPS v2 in the development of the subsequent models. We emphasize that including different rainfall products in the model has no strong effect on the remaining model parameters that determine the influence of slope and lithozones and road widening (see Table S1). In other words, while determining the overall performance, the choice of rainfall product does not affect our results and conclusions about the topographic and lithologic controls on landslide occurrence.

Figure 3Accumulated rainfall/precipitation amounts from different gridded products between 1 January and 10 October 2023. The sources of these data products are listed in Table 2. The black line indicates the road between Rishikesh and Joshimath (see Fig. 1).

Bayesian log-linear modeling of the landslide density (Figs. 4 and 5a and b) reveals a credible influence of the covariates rainfall (Fig. 5c), slope (Fig. 5d), lithozones (Fig. 5e) and road widening (Fig. 5f; see Fig. 4 for posterior means and 95 % highest density intervals and Fig. 6 for individual effects). A Bayesian feature rank algorithm based on the absolute magnitude of the parameters in each posterior sample (Makalic and Schmidt, 2011) ranks slope as the top covariate in terms of explanatory power together with rainfall, followed by road widening and lithozones (Table 4). Among the lithozones, zones 2 and 4 stand out as categories improving the explanatory power of the model. The individual effects of the covariates reveal a positive influence of rainfall and slope on landslides (Fig. 6a and b). Predictions of landslide densities in lithozone 4 are credibly lower than in lithozone 2 (Fig. 6c). Also, our model suggests that road widening affects the density of landslides. Corrected for the influence of other parameters, average densities are twice as high in widened road sections (Fig. 6d). The spatial pattern of predicted landslide density (Fig. 5g) is consistent with observed spatial density variations, but the high small-scale variability reflects the importance of slope as a predictor variable.

Figure 4Posterior parameter samples of the log-linear model of landslide occurrence along the NH-7. CHIRPS v2 is the gridded rainfall product, gradient determines the slope within 210 m of the road and lithozones were aggregated from a geological map. Note that lithozone 1 is missing since the parameter is encapsulated in the intercept.

Download

Figure 5Predictor and response variables used in the model. Panel (a) shows the occurrences of the fully and partially road-blocking landslides together with the elevation profile of the road. (b) Landslide density along the road measured using histogram binning (bin width is 9.8 km). (c) Accumulated rainfall between 1 January 1 and 10 October 2022 from CHIRPS v2. (d) Mean upslope gradient within a distance of 210 m of the road. (e) Lithozones along the road (see also Fig. 1). (f) Widened road segments. (g) Predicted landslide density using a model involving rainfall, slope, road widening and lithozones as covariates.

Download

Figure 6Main effects of predictors in the log-linear model of road-blocking landslides along the NH-7. (a) Effect of accumulated rainfall amount between 1 January and 10 October 2022 on the occurrence of landslides, averaging out the effects of the other predictors. The orange lines indicate the occurrences of landslides. (b) Effect of hillslope gradient, (c) of lithozone and (d) of road widening.

Download

Table 4Summary of posterior distributions of parameters fitted by the Bayesian log-linear point process model. Ranks in the model were calculated according to Makalic and Schmidt (2011).

Download Print Version | Download XLSX

The AUC is an aggregated metric for a point pattern model across thresholds and ranges between 0.5 and 1, with values close to 1 indicating high performance of the model. Our log-linear model has an AUC value of 0.79 (Fig. 7a). The inhomogeneous K function shown in Fig. 7b quantifies the expected number of points as a function of distance from each point, adjusted for the modeled inhomogeneous intensity of the point pattern. Distances between individual landslides are calculated as the shortest-path distance along the road rather than the direct Euclidean distance. Acceptance intervals around the theoretical K function were derived from repeated simulations of the inhomogeneous Bayesian log-linear model. The actual point pattern's K function is within these acceptance intervals, suggesting that the landslide locations do not exhibit clustering. A comparison of 100 randomly simulated and actual point densities (Fig. 7c) shows that the modeled and observed spatial landslide densities are consistent, although the second, smaller peak of landslide densities close to Joshimath is not captured by the model well. Moreover, 5-fold cross-validation based on 15 km long road segments reveals AUC values between 0.77 and 0.78.

Figure 7Evaluation of the log-linear model including rainfall, slope, lithozones and road widening. (a) Receiver operating characteristics (ROC) curve. The area-under-the-curve (AUC) metric is 0.79 (0.76–0.82, 95 % bootstrap confidence intervals). (b) The inhomogeneous K function corrects for the influence of an inhomogeneous Poisson point process model and tests for second-order effects (e.g., spatial clustering). Acceptance intervals of a theoretical model with independent point independence are shown in gray and were computed from 100 bootstrap samples. The red line is the empirical inhomogeneous K function. As the curve remains within the bootstrapped acceptance intervals, we retain the null hypothesis of inter-point independence. (c) Comparison of observed landslide densities (black line) with densities obtained from 100 random realizations (gray lines) from the model.

Download

Can the model be improved by incorporating more explanatory covariates? Our impression in the field was that landslides detach independently of planform hillslope geometry, as they occurred on spurs and in convex hollows. Nevertheless, we calculated total curvature and topographic roughness as potential predictor candidates. In addition, we used land cover (Table 3) and distance to faults (Fig. 1), as they are commonly used in susceptibility studies (e.g., Stanley and Kirschbaum, 2017; Li et al., 2020; Ozturk et al., 2021) and potentially control the density of roadside slope failures. Yet including these metrics in the model barely contributes to improving the model fit, and their incorporation in the model would, according to the constancy of the AIC, lead to overfitting (Fig. 8).

Figure 8Forward stepwise selection of additional explanatory covariates in the log-linear model of road-blocking landslides.

Download

We recorded more than one landslide per road kilometer along the NH-7 between Rishikesh and Joshimath. The fact that this road is strongly affected by landslides has been previously described and attributed to the region's fragile slopes, intense rainfall and frequent seismicity (Sati et al., 2011). In addition to the environmental conditions, road construction and widening have contributed to the occurrence of new landslides, which are often shallow and small, but nevertheless lead to fatalities, inflict severe damage on infrastructure and cause traffic disruption (Sati et al., 2011). We conducted a systematic survey of roadside landslides and derived a statistical model that quantifies landslide susceptibility along the NH-7 at a high spatial resolution.

Our analysis relied on a GPS-based survey of landslides while traveling from Rishikesh to Joshimath shortly after a period of anomalously high rainfall. Using this approach, we mapped landslides irrespective of cloud cover and without acquiring high-resolution satellite imagery, which is usually needed to reliably detect small landslides. A drawback, however, is that we may have missed landslides where debris had already been removed by road work. Also, detailed mapping of the areal extent of the landslides was not possible during drive-by, so we did not quantify the size of landslides. Thus, our analysis treats all landslides the same, irrespective of their areal extent and volume. To this end, however, this enables us to adopt a modeling approach that conceptualizes landslides as an unmarked network-attached point pattern (Baddeley et al., 2021; Okabe and Sugihara, 2012). Representing landslides as network-attached points and not as areal features has advantages and disadvantages. Constraining landslide locations to lie on roads demands that all predictor variables also need to be mapped to the road network, which entails some generalization and additional degrees of freedom about the choice and aggregation of 2D variables. For landslide susceptibility analysis, for example, this means that spatial variables characterizing the source area (e.g., hillslope gradient) are projected onto the road. At the same time, model development and fitting benefits from smaller sample sizes so that data amounts are moderate and computational demands during Bayesian posterior sampling remain manageable.

We detected a profound difference between rainfall products, a detailed analysis of which is outside the scope of this study. Several studies along the Himalayan orographic front have previously detected these differences. They have been attributed to the sparse network of rain gauges or, in cases where rainfall estimates are based on remote sensing data, to irregular acquisition times, which can result in missing individual rainfall events (Andermann et al., 2011; Hu et al., 2016). To this end, these uncertainties are detrimental to accurately capturing the spatial patterns of landslides (Ozturk et al., 2021). We found that CHIRPS v2 performed best in predicting the spatial landslide patterns along NH-7, but the search strategy employed must be viewed critically as the reverse conclusion, that landslides are controlled by these patterns, is not necessarily true. Indeed, studies come to different conclusions about the performance of CHIRPS v2 and other gridded rainfall products. For example, Kumar et al. (2021) studied different gridded rainfall products in the eastern Himalayas and found that CHIRPS v2 overestimated the monsoon but underestimated annual precipitation. Based on the analysis of 18 extreme precipitation events during 2014–2016 in the northwest Himalayas (including our study site), however, Jena et al. (2020) concluded that CHIRPS v2 provides the most reliable precipitation estimates. Indeed, the two-peaked rainfall pattern of CHIRPS v2 is most consistent with long-term rainfall patterns obtained from the interpolation of 44 rainfall gauge records averaged over the period from 1901 to 1950, which show highest values along the Himalayan front and along the physiographic transition to the Higher Himalaya (Basistha et al., 2008; Bookhagen, 2010).

Lithozones derived from a geological map contributed to the explanatory power of the model, thus highlighting the role of rock properties in modulating landslide susceptibility. As we did not measure the actual geotechnical and geomechanical properties or, e.g., bedding and foliation along our route, we can only provide a first-order reasoning of the prediction capacity of the lithozones. The high landslide density in lithozone 2 is likely related to the pronounced fissility and cleavage of the dominating shales and phyllites. Moreover, material softening, percolation and weathering cause a general decrease in rock strength. Tectonic activity adds to a general decrease in rock strength by creating shear surfaces with low friction angles (Stead, 2016). In addition, road segments where the adjoining hillslopes parallel bedding, joints or foliation planes are particularly vulnerable (e.g., Bartarya and Valdiya, 1989). Conversely, lithozone 4 is characterized by quartzite and igneous rocks. These have undergone low- to high-grade metamorphism and are generally harder and have a more irregular fabric that restrains the formation of planar slide surfaces. Moreover, these rocks tend to develop more-stable regolith mantles (Gerrard, 1994) and are thus less susceptible to landsliding. Our model shows that under the same topographic conditions and rainfall amounts, the rock types of lithozone 2 are 2–6 times more susceptible to landslides than those in lithozone 4 (Fig. 6c). The remaining lithozones are not credibly different from each other. In part, the lack of differences can be attributed to the low number of samples, as only a small fraction of the road traverses these lithozones (Table 2). Notwithstanding, a general trend towards lower landslide susceptibility from lithozones 1 to 5 is consistent with a previous review study about lithological controls on the occurrence of mass movements in the Himalayas (Gerrard, 1994).

Previous studies indicate that human activities have played a crucial role in predisposing slopes to failure (Li et al., 2020; van Westen et al., 2008; Tanyaş et al., 2022). Our analysis underscores these findings and quantifies the occurrence of landslides as being twice as high along widened sections of NH-7 if other predictor variables are held constant. The process of cutting into mountain sides to create wider pathways often creates unstable slopes that are prone to failure (Barnard et al., 2001; Haigh and Rawat, 2011; Sati et al., 2011; Li et al., 2020). Clearing trees and vegetation for road widening eliminates their stabilizing influence on slopes, thereby increasing landslide hazard (e.g., Guthrie, 2002). Additionally, widened roads can lead to increased surface runoff during heavy rainfall or snowmelt, saturating the soil and making it more susceptible to landslides (e.g., Wadhawan et al., 2020). The rock blasting required during the road widening process can lead to the fracturing and weakening of rock masses, creating potential landslide-prone zones along the road corridors (Sati et al., 2011). Moreover, road widening alters natural drainage patterns and potentially redirects water flow to adjacent slopes, thus causing water saturation, erosion and instability. In fact, these disturbances have previously led to frequent landslides along the NH-7: Sati et al. (2011) also report about ∼300 landslides occurring along the road more than 10 years ago. Our data indicate that 30 % of the recorded landslides are reactivated slope failures, which highlights the fact that slopes are recurrently unstable during periods of intense rainfall (Joshi and Kumar, 2006). During mapping, we also noticed that some slopes were engineered during the last years with retaining walls, yet many of these also failed.

Clearly, our model may miss important predictor variables that control the occurrence of landslides. We included variables that characterize environmental conditions and found that slope, rainfall, road widening and lithology largely explain the variability in landslide density. Variables such as land use or topographic derivatives do not improve the performance of the model as measured by the AIC, at least at the spatial scale at which these variables were available. However, small-scale topographic changes due to construction or land use changes (e.g., abandonment and degradation of agricultural fields and terrace systems) may exacerbate roadside slope failures (Jaquet et al., 2015; Mauri et al., 2022). A more-up-to-date DEM with higher resolution may indeed help to improve the spatial prediction of landslide, although higher digital terrain model resolutions were shown to not necessarily improve model performance, in particular along roads (Penna et al., 2014). Moreover, the propensity for landslides along widened segments of the road may change over time. Road widening and slope undercutting can be viewed as a perturbation to which slopes will adjust. Timescales of this adjustment may vary and depend on several factors. Likewise, remedial measures such as slope enforcement, artificial drainage and retention walls will affect slope stability. Including these activities of landslide mitigation was not possible in our study but could support planning of structural measures.

Our model hindcasts the spatial pattern of road-blocking landslides, and we posit that it can be used as a time-predictive model as well. Rainfall is one of the main covariates in the model and is also the one with the largest uncertainties, as shown by the discrepancies between the gridded rainfall products. A denser network of rain gauges and better availability of these data would likely contribute, together with weather forecasting, to more accurate estimates of landslide occurrences, which ultimately would facilitate more efficient allocation of resources for road maintenance. Also, recurrent slope failures should be monitored more closely to direct efforts for slope reinforcement. The land cover data, which we included in the model, is too coarse to capture the widespread lack of vegetation along the road. As many of the landslides were shallow, revegetating slopes may contribute to their stabilization (Vergani et al., 2017).

The NH-7 is a key arterial road, and landslides make transport of goods and people difficult, thus causing serious economic disruption. Moreover, slope failures along the road have led to fatalities in places where roads were widened, and recent heavy rainfall in July 2023 has caused several severe landslides (https://www.hindustantimes.com/india-news/ landslide-on-char-dham-road-leaves-pilgrims-stranded-heavy-rains-cause-damage-and-blockages-in-uttarakhand-101690225215947.html, last access: 26 June 2024). Such damage and commensurate fatalities may become even more frequent in the future. The entire upper Ganga Basin is susceptible to extreme rainfall events (Joshi and Kumar, 2006), and climate change projections – although subject to high uncertainties – indicate a trend towards more-frequent extreme events due to elevation-dependent warming and a likely increase in summer monsoon precipitation by 4 %–25 % (Krishnan et al., 2019). In addition, exposure to landslides is likely to increase in Uttarakhand as the Char Dham National Highway project is implemented (Chouhan et al., 2023). Road construction and increased traffic volumes attract more people, who will strive for new economic opportunities associated with sites along roads (Fort et al., 2010; Chouhan et al., 2023). These sites are often more susceptible to landslides as construction often implies vegetation removal and slope destabilization (Petley et al., 2007; Li et al., 2020). A reduction in traffic may disrupt the cycle of increasing hazard and exposure. The commissioning of the currently constructed 125 km long broad-gauge railway between Rishikesh and Karnprayag (Azad et al., 2022) might be a major step towards this goal.

The code and the data to run the analysis are available in the Supplement.

The supplement related to this article is available online at: https://doi.org/10.5194/nhess-24-3207-2024-supplement.

All authors conceptualized the study and conducted mapping landslides. RKG retrieved and analyzed the rainfall data. JM and AP compiled, analyzed and visualized the data, and WS conducted the statistical analysis. All authors wrote the paper.

The contact author has declared that none of the authors has any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

This article is part of the special issue “Estimating and predicting natural hazards and vulnerabilities in the Himalayan region”. It is not associated with a conference.

This research has been supported by the Deutscher Akademischer Austauschdienst (grant no. 57553291) and the Prime Minister's Research Fellowship, government of India (grant no. PM-31-22-695-414), awarded to Ravi Kumar Guntu.

This paper was edited by Olivier Dewitte and reviewed by two anonymous referees.