Three decades of landslide activity in western Nepal: new insights into trends and climate drivers

In recent decades, landslide disasters in the Himalayas, as in other mountain regions, are widely reported to have increased. While some studies have suggested a link to increasing heavy rainfall under a warmer climate, others pointed to anthropogenic influences on slope stability, and increasing exposure of people and assets located in harm’s way. A lack of sufficiently high-resolution regional landslide inventories, both spatially and temporally, has prevented any robust consensus so far. Focusing on Far-Western Nepal, we draw on remote sensing techniques to create a regional inventory of 26,350 single landslide events, of which 8778 date to the period 1992–2018. These events serve as a basis for the analyses of landslide frequency relationships and trends in relation to precipitation and temperature datasets. Results show a strong correlation between the annual number of shallow landslides and the accumulated monsoon precipitation (r = 0.74). Furthermore, warm and dry monsoons followed by especially rainy monsoons produce the highest incidence of shallow landslides (r = 0.77). However, we find strong spatial variability in the strength of these relationships, which is linked to recent demographic development in the region. This highlights the role of anthropogenic drivers, and in particular road cutting and land-use change, in amplifying the seasonal monsoon influence on slope stability. In parallel, the absence of any long-term trends in landslide activity, despite widely reported increase in landslide disasters, points strongly to increasing exposure of people and infrastructure as the main driver of landslide disasters in this region of Nepal. By contrast, no climate change signal is evident from the data.


Introduction
Landslides are a major cause of destruction to lives and livelihoods across the globe, with notable hotspots of landslide risk located along the great mountain chains of the world, particularly the Himalaya (Nadim et al. 2006;Petley 2012;Froude and Petley 2018). Here, and elsewhere, heavy rainfall is often associated with landslide triggering (Kirschbaum et al. 2015), with events predominating during the peak monsoon months of July-September. This clear link has led to increasing focus on climate change as a driver of past and future landslide activity (Gariano and Guzzetti 2016), particularly given expectations for increased rainfall extremes over much of the globe (Orlowsky and Seneviratne 2011), including the anticipated intensification of the South Asian Monsoon (Kitoh 2017).
However, despite reported increases in landslide fatalities and disasters (Stäubli et al. 2018), studies have been cautious to distinguish trends in actual landslide frequency, and attribution of any trends to climate change has hardly been possible. Firstly, the large-scale disaster databases upon which most landslide trends have been established may suffer from recording biases and have been shown to suffer from under-reporting of up to 2000% (Petley 2012;Kirschbaum et al. 2015). Secondly, even if disaster trends are considered robust, the difficulty comes from needing to distinguish any change in landslide activity (event frequency or magnitude) from other drivers of risk such as population growth and expansion of infrastructure into exposed areas. For example, in Peru, an increase in people affected by landslides has been reported, despite no clear trend in landslide frequency, strongly suggesting that enhanced vulnerability and exposure are driving landslide risk (Huggel et al. 2015). Conversely, in Italy, a decrease in large landslide disasters has been linked to reduced vulnerability and exposure, owing to improved information and warning systems (Gariano and Guzzetti 2016). Thirdly, even where significant trends in landslide frequency are identified, other non-climatic drivers of landslide conditioning and triggering may dominate. In Nepal, where a comprehensive database of fatal landslides has been maintained since the 1980s, an overall upward trend in fatal landslides is considered robust (Petley et al. 2007;Petley 2010), corroborating with the disaster inventories (Stäubli et al. 2018). While interannual variability in landsliding appeared to be linked to monsoon activity, Petley et al. (2007) concluded that poorly planned road building was likely a key driver of this trend, highlighting broader issues around regional development in Nepal and the other mountain regions (Hearn and Shakya 2017;Sudmeier-Rieux et al. 2019).
A general lack of robust evidence on past trends in landslide activity and associated climate drivers is reflected in the assessments of the Intergovernmental Panel on Climate Change (IPCC). In terms of detection and attribution, i.e., the formal process of linking climate change (and eventually anthropogenic influences on the climate) to physical impacts, the IPCC has remained essentially silent on landslides, concluding that there is no clear evidence that their frequency or magnitude has changed over the past decades (Cramer et al. 2014). This lack of evidence is crucial, as detection and attribution underpins issues of responsibility and compensation for climate change impacts (Huggel et al. 2020) and helps firmly establish those impacts that are most likely to become key risks under future warming.
In response to a clear need to better understand trends and drivers of landslide activity in monsoon-affected regions of the Himalaya, our study draws on a comprehensive inventory of 26,350 landslides in Western Nepal, and ERA-5 climate reanalysis for the period 1992-2018. Specific objectives of the study are to detect trends in landslide activity, both spatially and temporally, and to determine whether these trends can be statistically linked to climate (or other) drivers. Ultimately, results from this study should strengthen our understanding of how landslide risk will evolve under a future climate.

Study site
The study area is entirely located in the westernmost Sudurpashchim Province, formerly known as the Far-Western Development region of Nepal (Fig. 1). The area of interest (AOI) partially covers the Seti sub-basin and has an extent of 6460 km 2 , with elevations ranging from 235 to 4280 m a.s.l. The AOI comprises three physiographic/geological regions (Fig. 1A) aligned approximately WNW-ESE and gaining in elevation from South to North: (i) southernmost, the Sub-Himalayan geological zone of the Siwaliks (200-2000 m a.s.l.) is composed of very thick (4000 to 6000 m) friable sedimentary successions of mudstone, sandstone, and conglomerate with upward coarsening sequences (Gansser 1964) that are severely affected by tectonics and weathering processes, owing to frequent flash floods and debris flows (Dhital 2015); (ii) the Middle Mountains region (600-3000 m a.s.l.), comprising the Proterozoic Lesser Himalaya geological sequence, formed by a range of sedimentary, low-grade metamorphic, and crystalline rocks; and (iii) northernmost, the High Mountain region with the Greater Himalayan Crystalline Complex reaching the highest elevations and steepest topographies. Here, facies comprise granitic intrusions and metamorphic and crystalline rocks, frequently serving as sliding surfaces of the uppermost layers of fragile, fractured, and weathered sedimentary lithologies (Hasegawa et al. 2008).
Climate is of the monsoon type, strongly influenced by the airflow coming from the Bay of Bengal (Pokharel et al. 2019). The rainy season is concentrated from June to September with sparse rains during the dry winter season (Dahal 2012). Total annual precipitation averages 1870 mm, of which 80% is concentrated during the monsoon season, although local anomalies exist (Thyer 1985). Temperature shows predictable variability across the study area, with changes controlled not only by geography but also largely by differences in elevation. Temperatures range from warm-humid in the tropical southern regions, moderate temperatures in the Mid Mountains, to cold temperatures in the northern High Mountains (Shrestha and Aryal 2010).
The steep topography exerts strong control on societal development in the region. In the Far-Western Region, 45.6% of the local population lives in poverty (ADB 2017) with limited access to basic services (CBS 2016). Traditional cultural systems, closely linked to religion and traditions, greatly limit development (UNFCO 2012). Likewise, demography in the study area has shown one of the largest relative increases in Nepal over the last two decades (Worldpop 2020). This development has led to considerable expansion of farmlands on steep slopes and the development of infrastructures. In particular, the expansion of the road network stands out, often being associated with substandard engineering practices (Sudmeier-Rieux et al. 2019).

Methods
To achieve our goal in a standardized and replicable framework, we make use of remote sensing techniques and statistical methods in a GIS environment using a process that involves the following: (i) generation of a regional landslide inventory by manual satellite imagery interpretation; (ii) dating of the inventory by automated remote sensing techniques; (iii) time-series correlation of landslides with climatic variables; and (iv) analysis of trends, both for landslides and for climatic variables.

Landslide inventory visual interpretation
Interpretation of the landslide inventory for the Seti basin was based on the visual and systematic high-resolution and freely accessible satellite imagery from Google Earth (CNSS-Airbus and Pléiades). Overall, 93 different sources for 79 different time slices between 2002 and 2018 have been used for interpretation. We discarded overlapping imagery with excessive cloud coverage, blurring effects or darkness for a better interpretation of landslide-related features. Additionally, GIS tools have been used for topology correction processes, as well as for the postprocessing of the features by spatial intersections, spatial joins, identification, area calculation, and summary statistics, among others, so as to create a comprehensive landslide inventory.
A mesh of 3500 × 2000 m was used to identify past landslide activity in each of the 1023 quadrants. Events were then digitized in the form of landslide polygons and added to the inventory. For the remotely based interpretation process, we employed a maximum scale of interpretation of 1:1000, whereas for digitalization, the scale was maintained between 1:800 and 1:2000, resulting in a scale of 1:1000 for the visualization of the output inventory. Owing to the The areas in the inset figure of Nepal refer to the administrative divisions of Nepal: 1, Province No. 1; 2, Province No. 2; 3, Bagmati Province; 4, Gandaki Province; 5, Lumbini Province; 6, Karnali Province; and 7, Sudurpashchim Province, formerly known as the Far-Western Region and where our study area is entirely located scale of interpretation and intrinsic difficulties of the area for landslide identification, which is densely vegetated, we identified, for the vast majority of cases, the top, lateral, and, approximately, lower boundaries of landslides. Therefore, each of the final landslide polygons bound completely the "zone of depletion" and partially the "zone of accumulation" as the areas which, together, include all other morphological features (IAEG 1990) (Fig. 2). Hence, these digitized zones include the surface of rupture in its full extent and more or less partially the depositional zone of a landslide. Due to this limitation regarding boundary area, in which the lower limit of the rupture area cannot be properly distinguished, we do not include landslide area in analyses. Instead, we use point locations for each landslide, which roughly corresponds to the depletion zone centroid. This location has been calculated automatically for all landslides and corresponds to the central part of the upper ¾ of the digitized area.
Landslides can be classified as shallow or deep-seated landslides, depending on their maximum vertical depth from the original ground surface to the surface of rupture (Hungr et al. 2013). In shallow landslides, the slip surface is situated within the soil mantle or at the contact between a soil and a shallow weathered rock. This sliding surface can be situated at depths ranging from less than 1 m to a few meters. Meanwhile, in deep-seated landslides failing as translational and rotational slides, the sliding surface is often at depths ranging from a few to tens of meters. Shallow soil failures are generally rapid and initiated primarily by periods of persistent rainfalls or short-lived yet intense storms (Montrasio and Valentino 2008). By contrast, the slower, deepseated landslide triggering is favored by the succession of more complex processes (Yerro et al. 2015), often associated with gravitational phenomena in tectonically active areas (Dramis and Sorriso-Valvo 1994). Accordingly, prior to assessing correlations between landslides and climatic factors, one has to distinguish shallow from deep-seated landslides. We developed an automated procedure to this end (Online Resource 1) so as to extract approximate landslide depth using an ALOS PALSAR DEM with a 12.5-m pixel resolution (ASF 2020) and calculated maximum depths of the depletion zone (h) at a cross section perpendicular to the sliding axis and passing through the calculated starting point of each landslide (Fig. 2). Finally, based on the resulting data and crosschecking with visual interpretation, we assigned a conservative threshold of 10 m to differentiate between shallow and deepseated types. Besides, for the choice of this threshold, due to the lack of field measurements, we used the rational criterion that most deep-seated landslides occur within the limits of the regolith layer, while deep-seated ones may involve the bedrock. The selected depth threshold, up to 10 m, coincides with the depth of the regolith below which most regolith slides occur in other areas of the Himalayas (e.g., Basharat et al. 2018) with a similar regolith depth distribution to our study area (Pelletier et al. 2016). Furthermore, this depth also coincides with that used by Dou et al. (2015) where an automatic DEM method is also used for landslide depth classification.

Automated landslide dating
The identified landslides were then seasonally dated using differences in the Normalized Difference Vegetation Index (NDVI) on consecutive Landsat 5, 7, and 8 satellite imagery between 1992 and 2018, using our automated, GIS-based approach (Online Resource 2). With this dating procedure, we can deduce the age of the events not for a precise moment, but as occurring within a monsoon season of each year. To evaluate landslides at the intra-seasonal scale, we compared imagery taken before the start and after the end of the monsoon season. To avoid possible influences from clouds and snow, Landsat datasets were selected from April 20 to May 30 for the period just before the monsoon, maximizing the azimuth of the sun and thereby minimizing errors related to topographic shadows. Noteworthy, we also performed a preprocessing of each Landsat image for radiometric calibration, topographic correction, and masking of clouds, cloud shadows, and snow areas.
The approach computes differences in the NDVI for each pixel over consecutive images, so as to detect evidence for the potential removal of vegetation by landslide processes (Montrasio and Valentino 2016). Thereby only pixel changes falling inside the interpreted landslide polygons are considered. In this way, if a sudden reduction in vegetation is found between the pre-and post-monsoon Landsat images, landslide initiation or reactivation is attributed to a given monsoon period (Fig. 3). In addition, we discarded large changes in NDVI in (i) flat areas, as they are more indicative of vegetation changes at valley bottoms, crops, or riverbanks; (ii) the case of possible changes that are directly related with a clearing of vegetation for road construction; or (iii) water differences stemming from variations in secondary or tertiary drainage on steeper slopes. Due to the Landsat pre-processing Fig. 2 Example from the inventory of a digitized landslide boundary. Identified parts of the interpreted landslides correspond to the following: 1, main scarp; 2, top; 3, tip; 4, foot; 5, toe; 6, zone of depletion; 8, lateral flanks; 9, accumulation zone; and 10, toe of surface of rupture. Note that the limits of the depositional zones are less clear than for the surface of rupture masking, the area considered for landslide dating was reduced to a portion of the total study area, such that the analyzed area in a given year results from the intersection between the area of interpretation in each image and the automated Landsat imagery available in a given year.
The performance of the automated method and the choice of the NDVI threshold were validated through a random preliminary selection of 500 equally distributed (both spatially and temporally), dated landslide events. To this end, Landsat Pan-sharpen images (1999-2018) with a 15-m resolution were compared visually before and after the monsoon to derive true positives (TP) and false positives (FP) ( Table 1). Best results were obtained with a NDVI difference percentile threshold >0.95 with positive predictive values (PPV = TP/(TP + FP)) equal to 0.900. Although lower thresholds (0.925) would provide a larger number of TP, they also identify an overly large number of FP, lowering the PPV value. A high number of FP are significant for tests with low specificity, and therefore there is a lower probability that a positive result is indeed a true positive. The opposite happens with higher thresholds, where too many landslide events would have been missed.

Climate variables
We use the ERA5 (ERA5 2020) precipitation and temperature climate datasets for the landslide-climate analysis due to the consistency of the time series in mountain environments (Scherrer 2020). Precipitation and 2-m air temperatures are part of the ERA5 2D surface analysis products derived from the 4D-Var data assimilation in CY41R2 of ECMWF's Integrated Forecast System (ECMWF 2016). The datasets cover the period 1979-present with a cell-size resolution of~30 km (0.25 × 0.25°). The hourly mean total precipitation rate-given as in kg m -2 s -1 -was converted to daily accumulated total precipitation (in mm) so as to allow calculation of monsoon (June-September) accumulated precipitation. Moreover, maximum accumulated monsoon precipitation has been calculated for 1, 3, 5, 7, and 10 consecutive days with the aim of testing the incidence of short-lived yet intense precipitation events on landslide triggering. Besides, we also converted monthly averaged 2-m air temperature (provided in°K) to monsoon average temperatures (in°C). These average and cumulative values were calculated for each cell size of the original data of our AOI, as well as for the entire Nepal, and also given as a weighted mean for both areas.

Statistical analysis between landslides and climate variables
We calculated Pearson correlation coefficients (r) to measure the strength of association and direction of linear relationship between landslide occurrences and climate variables (Pearson and Pearson 1922), to obtain values between +1 (perfect positive correlation) and −1 (perfect negative correlation), or a value of 0 in the case of non-correlation. Prior to the use of Pearson's correlation for our time series, we verified that the distribution of each pair of variables was normal by performing a Shapiro-Wilk test (Shapiro and Wilk 1965) and to remove possible outliers. Besides, we selected a 95% confidence level (alpha = 0.05) for all statistical analysis. Next, we searched for a monotonic upward (or positive) or downward (negative) trend in the data by using the nonparametric and distribution-free Mann-Kendall test (Gilbert 1997). The null hypothesis for this test is that there is no trend over the time period, while the alternative hypothesis was that an upward or downward trend exists in the one-sided test. The Mann-Kendall test statistics (S) is given by Eq. (1): where x j and x i are the time series and where sgn (or sign function) is expressed with Eq. (2): In cases where S > 0, later observations in the time series tend to be larger than those that appear earlier in the time series, whereas the opposite is true in cases where S < 0. Therefore, if the alternative hypothesis is true with the p value being less than our significance level value (alpha), our data will follow an upward (S > 0, p value <alpha) or downward (S < 0, p value <alpha) trend. The Mann-Kendall test serves the identification of statistically significant trends in time-series data, whereas the Sen slope estimator (Theil 1950) is used to measure the magnitude of the trend. The latter is nonparametric as it does not follow any probability distribution and is non-sensitive to outliers. Magnitude of the trend is given by the slope, which is calculated by choosing the median of all lines through all the pairs of sample points in the time-series data (Eq. (3)).
Sen ′ s slope = Median (3)where x j and x i are data values at time j and i, respectively.

Results and discussion
Multi-temporal landslide inventory Image analysis and landslide digitization yielded a final inventory with 26,350 different events between 1992 and 2018 (Fig. 4). We find that translational, rotational, compound, and flow slides as well as debris flows (Highland and Bobrowsky 2008) are the main landslide types found in the area. The total inventoried landslide area covers 139 km 2 (2.17% of the study area) with a density of 4.08 landslides km -2 . Such a density can be considered unusually high when compared with other recent regional landslide inventories of similar resolution in Nepal (e.g., Dahal 2013;Meena et al. 2019). The vast majority of events (19,661, 74.6%) are interpreted as shallow landslides, whereas 6689 (25.4%) are deep-seated landslides. A conservative width/length ratio (W/L) threshold of 0.23 (Mikoš et al. 2017) was then applied to distinguish channelized (<0.23) from non-channelized landslides (≥0.23). Only 10.8% of all landslide falls show a W/L ratio < 0.23 and can therefore be considered possible channelized events, likely in the form of flow slides or debris flows. The remaining 89.2% exhibit a larger W/L ratio with less confined paths. In addition, many of the nonchannelized landslides correspond to reactivations on previous landslide surfaces. One-third (33.9%) of all landslides analyzed are interpreted as complex events occurring on previously unstable slopes. In the Seti basin area, recurrent deep-seated landslides show an evident spatial correlation with the proximity of geological lineaments that are interpreted as possible active tectonic lineaments (Online Resource 3).
The automated landslide identification method developed in this study permits the annual dating of 8778 landslides (33.3% of the inventory), for 23 different time slices between 1992 and 2018 (Fig. 5). All landslides from the inventory that could not be dated have likely been active prior to 1992. Of all dated events, 6856 (78.1%) correspond to shallow landslides, whereas 1922 (21.9 %) are deep-seated events. Therefore, the relative percentages are roughly maintained with respect to the overall inventory. One can thus infer that there is no relative increase in a specific type of landslide type during the last decades or that the automated dating tool is equally efficient in identifying landslides of both types. Using area-normalized numbers (Fig. 5b), the years 2007 and 2018 stand out with the highest numbers of landslides; the lowest number of landslide incidences is found in 2001, 2004, 2005, 2006, 2014, and 2015. The general spatial distribution of landslides shows a higher density of events in the northern half of the study area, where population density is also highest (CBS 2016). For the subsequent statistical analyses, initiation points of dated landslides have been regrouped spatially to the ERA5 grid cell in which they are located (Fig. 5a).

Correlations between landslide occurrences and climate variables
Linear Pearson correlations between the total number of landslide events per year and different weighted average ERA5 variables have been analyzed for the entire AOI (Fig. 6). Results show quite clearly that the climate variables considered in this study have a larger influence on rapid shallow than on deep-seated landslides, yielding high correlations of landslides with total cumulative precipitation. Total monsoon precipitation (TP) has a high correlation with shallow landslides (r = 0.74, p value = 3.99e−05) but no statistically significant relationship with deepseated landslides. These findings are in line with work evidencing that physically based models testing different rainfall thresholds perform reasonably well in the prediction of shallow landslides but lose effectiveness to predict deep-seated landslides (Guzzetti et al. 2007). This can be explained by the fact that deep-seated landslides typically require more complex prior conditioning factors than shallow landslides where soil water saturation alone is already able to trigger slope failure (Starkel 1996;Terlien 1997). Furthermore, numerical models have shown how high temperatures could increase the susceptibility to shallow landslides by increasing heat conduction in the soil and increasing the mobility of infiltrated water (Liu et al. 2019). Our results indeed demonstrate that warmer summers (T) have a certain influence on shallow landsliding (r = 0.41, p value = 0.045) but no significant influence on the occurrence of deep-seated landslides.
Climate-landslide relationships are also stronger for shallow events when considering the influence of precedent climatic conditions causing later failure. We find that dry monsoons (PTP) lead to an increasing landslide frequency during the next monsoon season (r = −0.41, p value = 0.049). Meanwhile, precedent temperature (PT) as an isolated factor does not correlate strongly with landslide incidence. Instead, if warm monsoons are followed by an especially humid monsoon (TPPT), the highest incidence of landslides occurs, reaching the highest correlation for both shallow (r = 0.77, p value = 4.19e−05) and deep-seated landslides (r = 0.42, p value = 0.047). In fact, this is the only statistically significant correlation found between climate and deep-seated landslides in this study. Therefore, precedent dry and warm conditions could facilitate soil water transport by enlarging macropores in certain soils and thereby increase landslide susceptibility under subsequent wet conditions (Tichavský et al. 2019). This result also suggests that projections of future landslide activity will need to consider changes in variability of extreme events, rather than focusing solely on changing rainfall magnitudes. Fig. 3 Automatic landslide dating tool: the landslide used in this example presents a complex evolution, a phenomenon that has been observed frequently in the area. This means that we do not identify a single landslide event; instead, three different consecutive events defined the current boundary of the landslide area. The automatic landslide tool has efficiently allowed for annual dating of such events

Original Paper
When comparing the correlations for maximum monsoon rainfall intensities over 1, 3, 5, 7, and 10 days, we obtain highest correlations for the longest time window of accumulated precipitation, i.e., 10 days (MAX10D). Shallow landslides show a remarkable correlation of 0.68 (p value = 0.00037), whereas deep-seated landslides exhibit very low correlation values, in  TP, monsoon accumulated precipitation; T, monsoon average temperature; TPPT, monsoon accumulated precipitation * (previous monsoon average temperature); MAX10D, monsoon maximum 10 days accumulated precipitation accordance with the differing hydrological characteristics between both landslides types. Loss of cohesion between surface layers in shallow translational landslides usually occurs due to rapid water infiltration, whereas in deep-seated landslides, the reduction of shear strength is linked to a steady rise of the groundwater level, which normally occurs after longer periods with abundant precipitation (Iverson 2000;Zêzere et al. 2005). Interestingly, correlations are not spatially homogeneous when analyzed at the 30-km scale corresponding to the ERA5 cell grid with important differences between the cells of the 18 different areas covering the AOI (Fig. 7). In general terms, for all the variables, higher, statistically significant (p value <0.05) correlations are observed in the northern areas 12 and 14 and to a lesser extent in area 13 (p values >0.05 and < 0.1). These three areas also exhibit the highest landslide densities for the entire study area. At the grid scale, as for the entire AOI, climatic variables exhibit strongest relationships with shallow landslides with less statistical significance in the case of deep-seated landslides. While there is a clear north-south contrast regarding the importance of precipitation for soil failure, we do not find equally clear limits for temperature as they show only a moderate positive correlation with landsliding across a large number of grid cells.
The observed spatial variability in the relationship between landslide activity and climate variables seem to be most strongly controlled by both demographic development and geological characteristics. Moreover, some influence of elevation can be found, particularly for deep-seated landslides, while slope inclination has the least influence (Table 2). Typically, at higher elevations, snow cover is more important and persists longer, increasing the depth of the moving regolith, and biasing the susceptibility of the terrain more towards deep-seated rather than shallow landslides (Anderson et al. 2009).
Various anthropogenic factors (Fig. 8) appear to enhance the relative influence of precipitation on landslide activity (r = 0.61). Those areas who experienced higher rates of demographic development between 2000 and2018 (i.e., areas 8, 12, 13, and14 in Fig. 8) inevitably become areas with greater changes in land use and further development of infrastructure, as obvious from the higher density of the road network (Meyer and Turner 1992). Similar factors have been identified widely in the literature as drivers of increased landslide susceptibility in mountain areas (e.g., Hung et al. 2015;Brenning et al. 2015). The particular geological characteristics of areas 12, 13, and 14 likely add in a decisive way to their susceptibility to landslide activity (Online Resource 4). Three geological formations are concentrated in these areas, accounting for 68% of the total shallow and 72% of the deep-seated landslides. These formations form an E-W corridor in which phyllites, slates, and schists are interleaved and are highly susceptible to soil failure, easily forming slip planes between the weathered upper materials and bedrock (Dhital 2015).

Fig. 8 Spatial distribution of anthropogenic factors controlling differential spatial distributions of landslide incidences
Trend analysis As a result of the positive correlations between some of the climate variables and landslide activity seen in the "Correlations between landslide occurrences and climate variables" section, we now look into possible longer-term trends of these variables. The idea here is whether a statistically significant positive or negative trend of one of the climatic variables could translate into a positive or negative trend in landslide frequency over time (Ávila et al. 2016).
Results across the entire study region do not show monotonic trend (Mann-Kendall) in the yearly total number of landslides, not even if shallow and deep-seated landslides are analyzed separately over the period 1992-2018, and even if p values exceeding 0.05 by far are applied (Fig. 9). As a consequence, results clearly do not reconcile with observations of global disaster catalogs (RED U 2020; Kirschbaum et al. 2010;Kirschbaum et al. 2015) for which an increase in landslide incidences is derived for the study area over the past two decades. This mismatch between the increase in landslide disasters, but not necessarily landslide events, provides robust evidence that the widely reported increasing trend in landslide disasters is driven primarily by changes in exposure-with the latter being ascribed to the demographic development in the region (Fig. 8). Trends for the same period for TP, TPPT, and MAX10D likewise do not show any statistically significant trends, consistent with previous studies that have been realized in the wider study region (Khatiwada et al. 2016). Meanwhile, mean temperature shows an average rise of almost 0.4°C, although the trend is not statistically significant (p value = 0.2).
Looking at spatial differences in trends using the ERA5 grid cells (Fig. 9b), we observe some positive trends in landslide incidence in some grid cells, but none of them match to those where higher correlations with the climatic variables were observed. It is thus likely that the localized increase in landslide incidences seen in these cases is associated with other, non-climatic drivers. While we are unaware of the reasons for the increase in landslide activity in areas 6, 9, and 11, it is reasonable to assume that the increase in areas 8 and 13 may have an anthropogenic component, as these are also the areas with the largest demographic growth (Fig. 8). What does seem clear is that a non-positive trend in climatic variables is Fig. 9 Mann-Kendall trends for the period 1992-2018: (a) landslide and climatic trends for the entire AOI and (b) Landslide and climatic trends per ERA5 grid cell seen in those areas with both the highest correlation between landslide activity and climate drivers, and the largest landslide frequency (areas 12 and 14). Hence, this could be an important reason for the resulting lack of a landslide trend observed throughout the study area. In other words, landslides are strongly related to interannual climate variability in this region, but to date, there is no clear climate change signal visible yet in landslide activity.
The extent to which these results from Western Nepal may be indicative for the rest of the country is unclear. Climatic trends calculated for the entire country for the period 1979-2018 show the contrast between our AOI and the rest of Nepal (Fig. 10). Mean values for the entire country present clear positive trends for total monsoon precipitation (p value = 0.0489) and average air temperature (p value = 4.778e −05). In total, average precipitation and temperature have increased by 190 mm and 0.6°C, respectively, over the past four decades. Noteworthy, the increase in monsoon precipitation and temperature has been most pronounced over the high mountain region-although not the most populated, this zone contains the sensitive cryosphere and numerous glacial lakes (Khadka et al. 2018); it also functions as an initiation zone for flash floods and cascading processes (e.g., landslide or glacier dammed lakes; Ruiz-Villanueva et al. 2016 Zheng et al. in press) that can threaten downstream communities (Allen et al. 2019). Furthermore, increasing trends have been strongest over the densely populated Central provinces, where Petley et al. (2007) have demonstrated a clear clustering of fatal landslide disasters. While we may reasonably speculate, given strong landslide-climate correlations observed in Western Nepal, that an increasing trend in monsoon rainfall over central regions of Nepal may at least partially be driving overall increasing landslide fatalities in Nepal (Petley et al. 2007;Froude and Petley 2018), the necessary robust multi-temporal data on landslide activity is lacking. Upscaling of the methodological approach from western Nepal to the entire country or broader Himalayan region should therefore be prioritized and could make use of citizen science approaches (Paul et al. 2017) to establish the initial high-resolution database of landslide polygons.

Conclusions
In this paper, we present new insights into landslide activity that was enabled through the creation of a high-resolution, multitemporal landslide inventory. Due to the high spatial variability observed, any analysis of linkages between landslides and climatic or other drivers should ideally be undertaken in such a systematic approach, applied at regional or larger scales. The combined use of manual and automated remote sensing interpretation techniques has also allowed us to create a unique historical inventory within one of the most landslide-prone areas of the world. Thereby, a total of 26,350 landslides have been identified, of which 8778 could be dated annually for the period 1992-2018 and used for statistical analysis in relation to climate variables. The dated shallow landslides show high correlation with both total monsoon precipitation as well as with shorter periods of 10 days of accumulated precipitation. Furthermore, we demonstrated how climatic conditions of the previous monsoon can have an important impact when warm monsoons are followed by particularly wet monsoons, leading to highest correlations with the observed occurrences of shallow landslides. Marked spatial variability in the strength of these correlations has been revealed, in which the geological characteristics, differences in elevation, but especially, the remarkable demographic development seem to act as key, magnifying the effects of seasonal climate variability on landslide frequency. Contrary to what is seen in disaster catalogs, and often reported in the media, we did not observe any increase in landslide events for the period 1992-2018. This strongly implies that the increase in reported landslide fatalities in western Nepal relates primarily to the notable demographic expansion in this area, placing more people and assets in harm's way. Hence, there would appear to be significant potential and need for community-focused disaster risk management strategies to reduce local exposure levels and minimize landslide losses. While no influence of climate change on landslide activity in the region can be demonstrated for the past decades, projected changes in extreme precipitation, if realized, will undoubtedly increase landslide frequency in the future. This will be of upmost concern in those areas where poorly engineered road cutting or other detrimental land-use practices continue unabated.