Characteristics of landslide path dependency revealed through multiple resolution landslide inventories in the Nepal Himalaya

14 Recent research in Umbria, Italy, has shown that landslide susceptibility is controlled by a process 15 called path-dependency, which describes how past landslides control the locations of future landslides. 16 To date, landslide path-dependency has only been characterised in Italy. This raises the question of 17 whether this process occurs in other geomorphic settings, and thus whether path-dependency should be 18 more universally included in landslide susceptibility assessments. Therefore, the main aim of this paper 19 is to investigate and quantify landslide path dependency in the Nepal Himalaya. This is achieved by 20 applying several path dependent metrics to three monsoon-triggered landslide inventories for the 21 central-eastern Nepal Himalaya. These inventories were developed at two different spatial and temporal 22 resolutions. As such, we aim not just to quantify landslide path-dependency, but also assess whether 23 path dependency characteristics are resolution-dependent. We find strong evidence that landslide path


Introduction
Landslides are a ubiquitous hazard in areas with mountainous terrain and frequent trigger events that can create significant risk when near populations and critical infrastructure.From January 2004 to December 2016, 55,997 people were killed globally in 4862 landslide events (Froude and Petley, 2018).
Furthermore, the impacts of landslides are typically greater in less economically developed countries owing to increased vulnerability and decreased resilience (Kjekstad and Highland, 2009).It is clear that landslides are a significant global hazard, and that further research is required to reduce the economic and societal toll of landslide occurrences.
Landslide mitigation typically relies upon assessments of landslide susceptibility, which quantifies the likely locations of future landslides (Guzzetti et al., 2006).Typical landslide susceptibility approaches are time-independent, meaning that landslide susceptibility is assumed to be stationary through time (Aleotti and Chowdhury, 1999;Samia et al., 2017;Varnes, 1978).However, this assumption is challenged by recent research showing that several processes impart time-dependent controls on landsliding.For example, it is now well-described that the damage signature of large magnitude earthquakes can transiently precondition landscapes such that they experience an increased likelihood of landslide occurrence for years to decades after the event (Marc et al., 2015;Parker et al., 2015).
Furthermore, recent analysis of a long-term landslide record in Italy suggests that past landslide events can transiently influence future landslide occurrence via a non-linear process defined by Samia et al. (2017) as path dependency (Samia et al., 2018(Samia et al., , 2017;;Temme et al., 2020).As established by Samia et al. (2017), path dependency describes how pre-existing landslides can impart a legacy effect that controls the occurrence and size of new landslides.Specifically, this legacy effect increases the propensity of new landslides to overlap or be contiguous with pre-existing landslides.Such pathdependent landslides are termed here as spatially associated landslides, whilst new landslides that do not overlap a pre-existing landslide are termed as spatially unassociated landslides.It is also observed that spatially associated landslides can have different geometric characteristics to spatially unassociated landslides, with the former being less round and larger on average (Samia et al., 2017).
Furthermore, it has since been shown that the inclusion of time-dependent path dependency in landslide susceptibility models, via various approaches, can significantly improve model performance (Samia et al., 2018;2020).However, the wider applicability of using path dependency in landslide susceptibility modelling remains uncertain, as path dependency is yet to be rigorously tested in other regions.As such, quantifying whether landslides exhibit path dependent behaviour in other geomorphic regions should be considered vital (Samia et al., 2017).One of the main reasons why this process has yet to be widely investigated is that it is thought to require large-scale multi-temporal landslide inventories to observe.
For example, the path dependency characteristics quantified by Samia et al. (2017) in Italy were derived from a 60-year multi-seasonal inventory across a 78.9 km 2 region.These landslides were mapped using stereoscopic aerial photographs ranging from 1:13,000 to 1:33,000 for the period 1941 -1997, and using field surveys and satellite imagery with better than 1 m resolution for subsequent years (Samia et al., 2017).Unfortunately, such regional multi-temporal inventories are uncommon as they are timeconsuming to develop, particularly over very large (> 100 km 2 ) regions (Reichenbach et al., 2018).As such, if path dependency is to be investigated more widely, it would be beneficial to have a greater understanding of what spatial and temporal resolution datasets are required to fully describe its behaviour.For example, over what time period should landslide data be collected and across how many time-slices?What size study region is required?At what resolution should landslides be mapped?A better understanding of these questions would allow other studies to target landslide data collection more efficiently for the purpose of investigating path dependency.
Therefore, this paper has two main aims.First, to investigate whether monsoon-triggered landslides exhibit path-dependency in the Nepal Himalaya and, if so, to quantify the characteristics of this process.Path dependency will be tested using several metrics that quantify the degree of overlap between new and pre-existing landslides following the approach of Samia et al. (2017).In addition, we hypothesise that, as observed in Italy (Samia et al., 2017), spatially associated landslides in the Nepal Himalaya will have different spatial and size characteristics to spatially unassociated landslides.Second, we investigate the spatial and temporal resolutions of landslide data that are required to observe landslide path-dependency.We do this by quantifying landslide path-dependency across three regions with two different resolution landslide inventories.One region with landslide data collected across a 45,000 km 2 region and 29-time slices, but with landslides mapped at a relatively low resolution of 15 -30 m.And two regions with landslide data collected across smaller (513 and 608 km 2 ) regions and just five time slices, but with landslides mapped at a higher (5 m) resolution.

Regional Setting
The main region selected for this study is a 45,000 km² area of the central-eastern Himalaya, Nepal (Fig. 1).Nepal was selected because it has an annual monsoon season that facilitates the systematic development of multi-temporal landslide inventories.Nepal also suffers frequently from other landslide trigger events such as earthquakes, storms and floods which, combined with its high relief, makes landslides and associated hazards pervasive (Cook et al., 2018;Dhital, 2003;Marc et al., 2015;Petley et al., 2007).The specific size and location of the main region was selected to ensure that our inventory included landsliding across a range of geomorphological and geological settings.The selected region includes catchments in all four of the main tectonic units that comprise the Nepal Himalaya: the Sub, Lesser, Greater and Tethyan Himalaya (DeCelles et al., 2001).These tectonic units have distinctly different geomorphologies and geologies across the study region.For example, the Tethyan Himalaya is dominated by granitic complexes (DeCelles et al., 1998) and has high average elevations of ~5000 m (Lavé and Avouac, 2001).Consequently, the Tethyan Himalaya are strongly influenced by glacial and periglacial processes, with diverse landforms such as active glaciers, moraines, alluvial fans, landslides and rock avalanches (e.g.Jones et al., 2020).In contrast, the bedrock geology of the Greater and Lesser Himalaya is composed of lower-grade metamorphic, metasedimentary and sedimentary rocks (Robinson et al., 2001;Upreti, 1999).Topographically the Greater and Lesser Himalaya are characterised by steep, fluvially incised valleys, including the Mahabharat Mountains and higher peaks to the north, which range from ~2000 -7000 m (Lavé and Avouac, 2001).Finally, the Sub-Himalaya are composed of Quaternary sandstones and conglomerates (Decelles et al., 2001;Robinson et al., 2001), with a topography that extends from the flat lowlands of the Indo-Gangetic Plain to the Siwalik Hills.Overall, by selecting a main region with a diverse set of landscapes, we can be confident that any observed path-dependency is not unique to any specific geomorphic setting, which, as described above, is a key aim of this paper.
Whilst the main study region allows us to investigate path-dependency across a range of geomorphic settings, we selected two-sub regions within the main region to assess how different inventory resolution impacts path-dependency results.These sub-regions are a 609 km² area following the Pasang Lhamu highway to Rasuawa Ghadhi (sub-region 1; Fig 1a and b) and a 513 km² area following the Araniko highway from Kathmandu to Kodari (sub-region 2; Fig 1a and c).These sub-regions were kept a similar size to allow the two to be easily compared, whilst their extents were defined within the main region to ensure that any differences between the main and sub regions are due to data resolution, rather than a change in landslide process.As with the main region, the specific locations of these two subregions were selected largely for their geomorphological and geological diversity, as well as their high landslide densities.The topography of both regions is extreme, with elevations ranging from 600 to 5800 m in sub-region 1 and 850 to 5900 m in sub-region 2, and hillslope angles in both regions ranging from 0 to 77° with an average of 32°.Across the two sub-regions there are seven lithologies: dolomites; granites; gneisses; phyllites; schists; quartzites, and shales.These rock types are characteristic of the Higher Himalaya Central Crystalline zone (Stöcklin, 1980).Both sub-regions contain large areas of granites and gneisses in the north, whilst the lithologies in the southern portions of both areas are more heterogeneous.This geomorphological and geological diversity allows investigation into how landslide distributions vary across a variety of spatially heterogeneous landslide control factors.

Data collection
To investigate landslide path dependency, the methods of Samia et al. (2017) were followed using three newly developed multi-temporal landslide inventories for each region described above.The following sections describe the main data collection and analysis methodologies used throughout this paper.

Landslide inventory mapping
Three monsoon-triggered landslide inventories were developed across central-eastern Nepal (Fig. 1).
The first, referred to as the "main region" inventory, was a multi-seasonal inventory developed for the ~45,000 km² main region described in section 2.0, using 15 -30 m Landsat 4/5/8 satellite imagery.This inventory included landslides mapped within 29 distinct time-slices across a 30 year period, where each time-slice encompassed a single monsoon season between 1988 and 2018, excluding 2011 and 2012 which could not be mapped due to Landsat 7 scanline errors (Alexandridis et al., 2013).In total, the main region inventory included 12,901 monsoon-triggered landslides (Supplementary Table 1).The second and third inventories, referred to as the "sub-region 1" and "sub-region 2" inventories, were developed across the two smaller sub-regions described in section 2.0.These sub-regions both fell within the main region, where sub-region 1 covers 609 km² along the Pasang Lhamu Highway, and subregion 2 covers 513 km² along the Araniko Highway (Fig. 1a -c).Both sub-region inventories were developed using 5 m resolution RapidEye satellite imagery and included landslides mapped within 5 distinct time-slices, where each time-slice encompassed two monsoon seasons (and the month's inbetween) from 2009 and 2019.This resulted in a total of 623 landslides in Sub-region 1 and 569 landslides in Sub-region 2 (Supplementary Tables 1 and 2).For all three inventories, landslides in a given time-slice were mapped through the visual comparison of pre-and post-time-slice satellite imagery.All used satellite images had < 15% cloud cover and were obtained during the Nepalese dry season (1st October -31st March).The specific satellite types and imagery dates used to map each individual landslides are shown in Supplementary Tables 1 -3.Landslides were mapped as polygons that included the scar, deposition and runout zones of each event, with care taken to avoid landslide amalgamation (Marc and Hovius, 2015) and mass-wasting due to anthropogenic activity such as cut-and-fill practices, land-clearing and road-tipping.For both resolution inventories, road-tipping events were identified as any mass-wasting event that occurred below and concurrently with a new road.We also ensured that our monsoon-triggered inventories did not include any of the coseismic landslides triggered by the 2015 Gorkha earthquake by discounting any landslides included in the comprehensive coseismic inventory of Roback et al., (2018).Within all three inventories, landslide spatial associations were classified into three categories (Fig. 2) using the ArcGIS functions: outside, partial, and inside.The outside function selected any landslide which at the time of its occurrence made no contact with an earlier landslide polygon.I.e., it defined all spatially unassociated landslides.The partial function selected any landslide which crossed the outline of an earlier landslide polygon, whilst the inside function selected any landslide which occurred fully within an earlier landslide polygon.I.e., these functions defined all spatially associated landslides, which were subdivided into spatially associated (partial) and spatially associated (inside) respectively.
Example landslides of each path-dependency type are shown in Figure 2.

Size-frequency and roundness analysis
Assessment of landslide size has important implications for the design and implementation of appropriate hazard mitigation and management strategies.In this case, it would be particularly useful to understand whether spatially associated landslides have similar size-frequency distributions to spatially unassociated landslides, and thus whether they need to be managed differently.As such, prior to analysing to landslide path dependency metrics, for our main inventory and combined sub-inventories we quantify the landslide area-frequency distributions for all landslide data and subsets of both the spatially unassociated and spatially associated (partial) landslide types.It should be noted that we did not calculate distributions for the spatially associated (inside) landslide types as there were very few cases of these within the inventories (46 in the main inventory and none in the sub-inventories).
Landslide size-distributions are typically characterised according to the probability density function of landslide areas p(AL), which is defined by Equation 1 (Malamud et al., 2004): Where NLT is the total number of landslides in the inventory, AL L is the number of landslides with areas between AL and AL L. Probability density functions of landslide area are commonly observed to exhibit a power-law decay and exponential roll-over at smaller areas that can be modelled by a three-parameter inverse-gamma distribution (Malamud et al., 2004), as defined by: (Equation 2) of the inverse-power law (i.e., the steepness of the right tail of the probability density function), rollover, which represents the size of the most frequent landslides, is typically used as a means of comparing the completeness of different inventories, as it indicates the landslide area below which landslide frequency-density values begin to decrease 2019).The exponent of the inverse power-law describes the rate at which the probability of getting larger landslides decreases.A larger exponent indicates that the probability of getting larger events is decreasing quickly, and thus that larger landslides are contributing less to each inventory.
Conversely, a smaller exponent indicates that the probability of getting larger events is decreasing more slowly, and thus that larger landslides are contributing more to each inventory (Borgomeo et al., 2014;Van Den Eeckhaut et al., 2007).
Here, we use the LANDSTAT software (version 10) to fit the three-parameter inverse gamma distribution (Malamud et al., 2004; Equation 2) to the probability density functions of the areafrequency data for three subsets (all landslides, spatially unassociated landslides, and spatially associated (partial) landslides) of our main inventory and combined sub-region inventories.This software uses Maximum Likelihood Estimation (MLE) to optimise the parameters of equations 1 and 2, and a bootstrapped (here with 1000 simulations) Kolmogorov-Smirnov (K-S) test to estimate parameter uncertainties and overall goodness of fit of the inventory data to the fitted distributions.
Finally, we also calculated the roundness of the landslides in both inventories.This was done using the methodology outlined by Samia et al. (2017), whereby roundness is defined as the ratio between the actual perimeter of a landslide, and the theoretical perimeter a landslide would have if it were a perfect circle with the same area.

Path dependency analyses
To quantify path dependency in the Himalaya, we apply the methods of Samia et al. (2017) to our three inventories.These methods are the overlap index, the unaffected area, and the number of overlaps.The overlap index quantifies the degree of overlap between landslides in two different inventory time slices according to equation 3 (Samia et al., 2017).

= ( )
(Equation 3) Where, t is the average time of the first time slice, t + n is the average time of the second time slice, and is the geometric intersection (overlapping area) between two time slices.Average times were taken as the middle of the monsoon period for the main inventory and the middle of the two-year time slice period for the sub-region inventories.By plotting the overlap index values against the respective time interval between time slices, the relationship between the overlap of earlier and new landslides and time since an earlier landslide can be established.For each inventory, the overlap index was calculated for every possible combination of time-slices, as well as between each year in each inventory and the 2015 Gorkha earthquake coseismic landslides of Roback et al. (2018).This gave 435 overlap index values for the main region inventory and 15 for each of the sub-region inventories.
The unaffected area method compares the actual area of land unaffected by landsliding to the theoretical area of land that would be unaffected by landsliding if no overlapping of landslides in different time slices occurred.The Actually Unaffected Area (AUA) is given by equation 4 (Samia et al., 2017): Where, t is the time slice number, ALTi is the total area of landsliding in time slice i, AS is the area of the study region, and U is the union of all landslides between time slice i =1 and time slice t.In other words, this metric calculates the total area of landsliding up to a given time slice, after accounting for overlaps, as a dimensionless fraction of the size of the study region.
Conversely, as shown in equation 5 (Samia et al., 2017) the Theoretically Unaffected Area (TUA) does not account for overlaps: = / (Equation 5) This metric essentially calculates the total area of landsliding up to a given time slice, without accounting for overlaps (i.e., assuming all landslides are spatially unassociated), as a dimensionless fraction of the size of the study region.As such, if spatially associated landslides are occurring, when the AUA and TUA are plotted against one another through time, we would expect the AUA to plot progressively lower than the TUA, with the distance between the two indicating the degree to which landslides are overlapping across the study region.These metrics were calculated at every time slice of each of our three inventories, i.e., calculated 29 times between 1988 and 2018 for the main region inventory, and 5 times between 2009 and 2019 for the two sub-region inventories.
The number of overlaps method compares for each inventory the actual area of landsliding at different degrees of overlap to the area of landsliding predicted at each degree of overlap by a random model.
The actual number of overlaps in each inventory was calculated by first converting the landslide polygons in each time slice into rasters using the ArcGIS 'Feature to Raster' tool.These time slices were rasterised at the same resolution as the landslide mapping (i.e., 30x30 m cells for the main region inventory and 5x5 m cells for the sub region inventories).All landslide presence cells were given a value of 1, whilst landslide absence cells were given values of zero.Then, by summing all time-slice rasters for each inventory, the total number of landslide presence cells (and thus area of landsliding) at different degrees of overlap were calculated for each inventory.The random model used a random number generator in Matlab to randomly assign the same number of landslide presence cells as existed in each time slice of the actual inventories to a grid with the same number of cells as each study region.
The randomly assigned time slices were then summed to count the total number of landslide presence cells (and thus area of landsliding) at different degrees of overlap for the random model.Each random model was run 50 times per inventory to get a mean number of overlaps.If landslides in each inventory are exhibiting path dependency, it would be expected that the actual data will show a greater area of landslides at higher degrees of overlap than is predicted by the random model.

Landslide control factor analysis
The second hypothesis to be tested in this paper is that spatially associated and spatially unassociated landslides have different spatial distributions and size characteristics.To investigate the spatial distributions, we first used standard ArcGIS tools to extract the slope, elevation, aspect, and lithology at the centroid of each landslide in each inventory, where all topographical data were derived from the ALOS World 3D-30m (©JAXA) product.We then split each of our inventories into four subsets: all data, just spatially unassociated landslides and just spatially associated (partial) landslides.Again, owing to the low number of cases, we do not include the spatially associated (inside) landslides in this analysis.We then used standard ArcGIS zonal statistics tools to calculate for each subset of the main inventory and the combined sub-region inventories, the percentage of landslides that occurred across different bins of each control factor.The percentages of the main inventory study region landscape and combined sub-region inventory landscapes that fell within each bin were also calculated.Finally, for the main inventory and combined sub-region inventories, the percentage distributions of each landslide subset were plotted against each other and the percentage distribution of the respective study region landscapes to compare how their distributions across each control factor varied.

Landslide inventory, area-frequency and roundness
The main region landslide inventory contained 12,901 landslides; 11,031 of which were spatially unassociated landslides, 1,824 of which were spatially associated (partial) landslides and 46 of which were spatially associated (inside) landslides.These landslides had a mean area of 12,616 m 2 , with a minimum area of 198 m 2 and a maximum area of 684,780 m 2 (Table 1).Sub-region 1 contained 623 landslides, of which 545 were spatially unassociated landslides and 78 were spatially associated (partial) landslides, with no inside cases.These landslide areas varied from 190 m² to 71,000 m² with a mean of 5,980 m² (Table 1).Sub-region 2 contained 569 landslides, of which 483 were spatially unassociated landslides and 86 were spatially associated (partial), with no inside cases.These landslide areas varied from 170 m² to 45,780 m², with a mean of 4,995 m² (Table 1).
To investigate whether spatially associated and spatially unassociated landslides have different size characteristics, three parameter inverse gamma distributions were fitted to the probability density functions of landslide area for all of the data, just the spatially unassociated landslides, and just the spatially associated (partial) landslides for both the main inventory and the combined sub-region inventories (Fig. 3).The scaling parameters and rollover values of these distributions are shown in Table 1 alongside the mean, minimum and maximum landslide areas for each case.This highlights several significant differences between both the two inventory types and the two landslide path-dependent types.In terms of the inventory types, we see that the main region distributions have higher scaling exponents ( ) of 1.62 -1.83 compared to the Sub-region distributions, which have exponents of 1.35 -1.45.Similarly, the main region distributions have higher mean landslide areas (12,000 -16,400 m 2 ) and rollover values (3260 -3440 m 2 ) compared to the sub-region inventory mean areas (5300 -6860 m 2 ) and rollovers (1090 -1170 m 2 ).In terms of spatially associated versus spatially unassociated landslides, we see that for both inventories the spatially associated (partial) landslides have higher mean areas and rollovers than the spatially unassociated landslides.In the case of the main region inventory, the spatially associated (partial) landslides have a mean area of 16,410 m 2 compared to 12,017 m 2 for the spatially unassociated landslides, and rollovers of 3444 m 2 compared to 3256 m 2 .Likewise, for the sub-region inventories, the spatially associated (partial) landslides have mean areas of 6861 m 2 compared to 5295 m 2 for the spatially unassociated landslides, and rollovers of 1171 m 2 compared to 1087 m 2 .The spatially associated (partial) landslides also have lower scaling exponents than the spatially unassociated landslides, with a value of 1.62 compared to 1.83 in the main region, and a value of 1.35 compared to 1.45 in the sub-regions.
Landslide roundness was also determined to compare differences between spatially associated and unassociated landslides, as well as between the main region and sub-regions to highlight inventory resolution differences (Fig. 4).The main region contains landslides with a higher average roundness than the sub-regions.The average landslide roundness is ~0.78 in the main region and ~0.65 in the sub-regions.In terms of spatial associations, spatially unassociated landslides have slightly larger average roundness than spatially associated landslides in both regions.This difference is most pronounced in the main region data, where spatially associated landslides have a mean roundness of 0.76 compared to 0.8 for those that are spatially unassociated.

Path dependency
Path dependency was investigated for all of our inventories using three metrics proposed by Samia et al. (2017); number of overlaps, overlap index, and unaffected areas.In terms of number of overlaps, all inventories are found to have a greater area of landsliding at higher degrees of overlap than is expected based on the random model (Fig. 5).This is particularly evident in the main region, where the random model never observes more than three degrees of overlap, compared to five degrees of overlap in the actual data.Furthermore, whilst the two sub-regions do not have more degrees of overlap than predicted by the random model, the actual data have 3-4 times more area of landsliding at three degrees of overlap than predicted, and approximately 1.5 times more landslide area at two degrees of overlap than predicted.
From the overlap index analysis, it is apparent that there is a weak negative coincidence between amount of overlap between landslides in two time-slices and the time period between those time-slices (Fig. 6).
In the main region, as time increases between time-slices the overlap index decreases gradually from 0.01 -0.02 at one year between time slices to 0 -0.005 at 30 years between time slices (Fig. 6a).The two sub regions show a similar trend, with a gradual decrease in overlap index from 0.01 -0.05 at one year between time slices to 0 -0.001 at eight years between time slices (Fig. 6b -c).Time-slices compared to the Roback (2018) Gorkha earthquake landslide inventory show a generally lower overlap index in the main region whilst maintaining a similar overlap index in the sub-regions (Fig. 6a -c).
Interestingly, whilst both the main region and sub-regions show the same overall trend, the characteristics of the trends appear to be resolution dependent, with the negative trend of the index values for sub-regions 1 and 2 significantly steeper than the trend in the main region (Fig. 6d).However, these observations are tentative, with significant variation in overlap index at different time intervals, and thus no conclusive R 2 values.
In terms of the unaffected area method, in the main region and sub-region 1, the actually unaffected area (AUA) gradually falls below the theoretically unaffected area (TUA) as time increases (Fig. 7ab).However, the divergence between the AUA and TUA in both regions is small, with a maximum difference of 3x10 -4 in the main region and 6x10 -5 in sub-region 1.Furthermore, in sub-region 2, the AUA and TUA remain similar throughout, with almost negligible divergence between the two metrics (Fig. 7c).It should also be noted that all three regions exhibit large perturbations in both the AUA and TUA that are coincident with cloud outburst storms in 1993 and 2002 (Fig. 7a) and the Gorkha earthquake in 2015 (Fig. 7a -c).

Landslide control factor distributions
To test whether spatially associated and spatially unassociated landslides have different spatial characteristics, the spatial distributions of the spatially associated and unassociated landslides in the main region and combined sub-region inventories were quantified with respect to four potential landslide controlling factors: elevation, hillslope angle, aspect, and geology (Fig. 8a -h).
It is apparent that the sub-regions have higher average elevations than the main region, with most frequent elevations of 2000 -2400 m compared to 800 -1200 m (Fig 7a -b).This is unsurprising given that the main region encompasses large parts of the Sub and Lesser Himalaya, whereas the two subregions are located towards the north of the main region in the more topographically extreme Greater Himalaya.Despite these differences, overall landslides are distributed similarly across both the main and sub-regions, with landslides typically occurring between 400 -3600 m.In the main region, this elevation distribution is similar for all landslides and both the spatially associated and unassociated landslide subsets (Fig. 8a).However, in the sub-regions, it is apparent that the spatially associated (partial) landslide types are 5 -10% more likely to occur at elevations of 2000 -2800 m than spatially unassociated landslides, and 2 -5% less likely to occur at higher elevations of 2800 -3600 m (Fig. 8b).
In terms of hillslope angle, both the main region and sub-regions have similar slope distributions that are normally distributed around 30 -35 o (Fig. 8c -d).In both the main and sub-regions, landslides also show a normal distribution with respect to slope, though shifted 5 -10 o higher than the study region distribution, such that landslides are 5 -10% more likely to occur at slope angles of 35 -55 o and less likely at lower slope angles.Spatially unassociated and spatially associated (partial) landslides occur on similar slopes in the main region (Fig. 8c).However, in the sub-regions, spatially associated (partial) landslides are 2 -4% less likely to occur at slope angles of 0 -20° than spatially unassociated landslides, but 2 -5% more likely to occur at slope angles of 25 -30 o and at slope angles > 50° than spatially unassociated landslides (Fig. 8d).
Hillslope aspect across both the main and sub-region study areas is evenly distributed across all orientations, with almost negligible amounts of flat terrain (Fig. 8e -f).However, landslides in both regions are 15 -20% more likely to occur on S and SE aspects.In the sub-regions (Fig. 8f), spatially associated (partial) landslides are ~5% less likely than spatially unassociated landslides on SE aspects, but up to 10% more likely on S -SW aspects.
Lithological variability is largely similar between the main region and sub-region study areas, with granites and phyllites dominating, though the sub-regions do not contain significant occurrences of marble, schist, quartzites, shales, or Quaternary sandstones and conglomerates (Fig. 8g -h).Generally, landslides in both regions are 10 -20% more likely to occur in phyllites, but 20 -30% less likely to occur in granites.In the sub-regions, spatially associated and unassociated landslides are similarly distributed across the different lithologies (Fig. 8h).Spatially associated (partial) landslides are slightly less likely than other types in Quaternary conglomerates and sandstone, but up to 5% more likely in granites and phyllites (Fig. 8g).

Does landslide path dependency occur in the Nepal Himalaya?
The first aim of this paper was to investigate whether landslides exhibit path-dependency in the Nepal Himalayas.To complete this aim we quantified the three path dependent metrics used by Samia et al.
(2017) (number of overlaps, overlap index and unaffected area) for three multi-temporal landslide inventories that covered the central-eastern Nepal Himalaya.For all three landslide inventories, these metrics suggest that there is strong evidence of path dependency in this region.
In terms of number of overlaps, all three of our inventories for central-eastern Nepal are found to have larger areas of landsliding at higher degrees of overlap than is predicted by a random model (Fig. 5ac).In other words, this shows that new landslides overlap with earlier landslides to a greater extent than would be expected if the spatial distributions of landslides through time were random.This result is indicative of path dependency and follows a similar general trend to that found by Samia et al. (2017) in Collazzone (Italy), who also found that the number of overlaps was greater than would be expected from a random model.
Yet, whilst the general trend in Collazzone may be similar to Nepal, there are some subtle differences.
Notably, the rate of decrease of landslide area with higher degrees of overlap in Collazzone is approximately half that observed in Nepal (-0.62 compared to -1.2 --2.0), suggesting that the Collazzone region has the potential to generate greater degrees of overlap given enough time, but that Nepal has a greater area of overlap than expected at lower degrees of overlap.It is unclear whether this reflects physical processes that differ between the two regions, or is simply an artefact of the different inventory characteristics used between the two studies.For example, one explanation could be due to a difference in landside typology and size between the two regions.In Italy, ancient deep-seated landslides are common (Cardinali et al., 2002) and included within the Samia et al. (2017) inventory.
Conversely, our inventories for Nepal are dominated by the shallow rockfalls and slides that are pervasive across the Himalaya (Dahal et al., 2013), with deep-seated relict events not included.Deepseated landslides are generally larger than shallow movements because there is more material available to be mobilised (Zêzere et al., 2005) and as such could have greater potential for higher degrees of overlap with new landslides.Therefore, it is possible that if we had included much older, large, deepseated landslides in our inventory we would have seen a similar rate of decrease to Collazzone.
However, if it is the inclusion of larger landslides causing the disparity, then we would expect the mean and maximum sizes of the Samia et al. (2017) inventory to be larger than ours.Yet this is not the case, with the Samia et al. ( 2017) inventory having larger average sizes than our sub-region inventories, but much smaller average and maximum landslide sizes than for our main inventory.An alternative explanation is that the rate of decrease of landslide area with degree of overlap is an artefact of inventory length.We note that the rate of decrease appears sensitive to inventory length, with the Samia et al.
(2017) inventory, our main inventory and our two sub-inventories having decreasing time-spans and decreasing rates.This is logical, as the longer the inventory, the greater the likelihood of observing increasing degrees of overlap.
Whilst the number of overlap results confirm that the Himalaya has more landslide overlap than expected if landslide occurrence was random, the overlap index and theoretical area analysis provide further quantification of the actual extent and characteristics of landslide overlap (and thus path dependency).The overlap index results tentatively support that path dependency occurs within the Himalaya, with a weak negative correlation between time between time-slices and overlap index across all three Nepal inventories (Fig. 6a -c).This trend suggests that the occurrence of spatially associated (i.e., overlapping) landslides relates to the time since earlier landsliding, with the likelihood of a new landslide overlapping with an earlier landslide decreasing with time.However, it should be noted that this relationship is not statistically significant, with overlap index showing large variations, particularly at 5 -15 years between time slices.
Again, these results follow a similar trend to those reported by Samia et al. (2017), who also observed a weak negative correlation between overlap index and time passed between time slices.However, the overlap index in the Collazzone region is at least a factor of 10 larger than in any of the Nepal Himalaya regions investigated here.As with the number of overlap results, we suggest that this is due to inventory differences in landslide typology and size.Indeed, whilst we did observe more overlaps than expected given a random distribution, as we only mapped recent shallow landslides, the overlaps between landslides were limited in their absolute potential area.Conversely, the first Samia et al., (2017) timeslice contained substantial relict deep-seated events, meaning that subsequent landslides had a high likelihood of occurring fully within these early failures.Whilst our main inventory and the Roback (2018) Gorkha coseismic inventory do include a couple of very large failures (e.g. the Langtang avalanche (Jones et al., 2020) and the Jure landslide (Regmi et al., 2017)) these occurred late in our time series, so there was less potential for subsequent large degrees of overlap.Interestingly, overlap index values between the Roback (2018) inventory and the time slices from the main region inventory are typically smaller than the index values between different slices of the main region inventory.This may be explained by the fact the Gorkha coseismic landslides had an overall smaller size distribution than expected for an earthquake of that magnitude (Collins and Jibson, 2015).
The unaffected area results also suggest that path dependency is occurring across Nepal, with the main region and sub-region 1 inventories showing that the AUA area diverges below the TUA over time (Fig. 7a -b).In other words, this shows that the occurrence of overlapping landslides is reducing the actually affected area of landsliding, as if landslides never overlapped then the actually and theoretically unaffected areas would be identical.Again, this result follows the same general trend as that observed in Collazzone.Interestingly, we find that sub-region 2 shows almost no divergence between the AUA and TUA.The cause of this lack of divergence remains unclear, and so should be considered in future work.
It is notable that the relative area values are several orders of magnitude higher in Collazzone than they are in the Nepal Himalaya, with a larger observed divergence between the actually and theoretically unaffected areas.This discrepancy is again likely explained by the difference in study region size and landslide inventories of our study and that of Samia et al. (2017).The region studied by Samia et al., (2017) was relatively small (~78 km 2 ) and their analysis included very large relict deep-seated events.
Conversely, our main study region was ~45,000 km 2 , whilst our sub-regions were ~513 -609 km 2 , and we only mapped recent small-scale landslides.Consequently, the landslide number densities per timeslice in Collazzone ranged from 0.21 -8.88 N/km 2 , whereas ours ranged from 0.003 -0.54 N/km 2 .As such, as this metric is proportional to study region size, the lower densities of landsliding in Nepal explain why our relative area values are so much lower.This also explains why the divergence between the actually and theoretically unaffected areas is lower for Nepal, because even though there are more overlaps than expected across both regions (e.g.Fig. 5), as a proportion of the study region size the actual area affected by overlaps is small.Furthermore, it is also observed that in all of the Nepal landslide inventories, major increases in both the AUA and TUA are coincident with cloud outburst storms in 1993 and 2002 (Fig. 7a) and the Gorkha earthquake in 2015 (Fig. 7a -c).This observation highlights that the occurrence of new landslides and the reactivation of pre-existing landslides is largely driven by extreme triggering events, rather than just gravitational processes.This observation corroborates other studies that have shown that shallow landsliding is commonly perturbed by infrequent, high magnitude landslide drivers such as earthquakes and storms (Chen et al., 2020;Dadson et al., 2004;Marc et al., 2015;Martino, 2015;Rathburn et al., 2017;Wolman and Miller, 1960).For example, the 2008 Wenchuan earthquake was observed to cause significant increases in post-earthquake landslide reactivation rates (Chen et al., 2020;Tang et al., 2016;Zhang et al., 2014).Therefore, our new data support previous observations showing that significant increases in rates of landslide reactivation, and thus path dependency, are most likely following extreme events.

Do spatially associated and unassociated landslides have different size, roundness and spatial distributions?
Overall, the number of overlaps, overlap index and theoretical area analysis show that landslide path dependency is occurring within the Nepal Himalaya.The second testable hypothesis is that path dependent (spatially associated) landslides will have different size and spatial distributions to spatially unassociated landslides, as found by Samia et al. (2017) in Collazzone, To investigate differences in spatial distributions we quantified how different landslide subsets (all landslides, spatially associated (partial) and spatially unassociated) in our main inventory and combined sub-region inventories were distributed according to four typical landslide control factors; elevation, slope, aspect and lithology.The first control factor to be considered is elevation; all landslides are proportionately more likely to occur at lower elevations (800 -2800 m) than at higher elevations (2800 -6000 m) in both study areas (Fig. 8a -b).This is unsurprising given that our inventories are of monsoon-triggered landslides, and that the lower elevation regions in the Himalaya typically experience more rainfall due orographic effects (Ichiyanagi et al., 2007).Furthermore, particularly across the main region, it is noted that landslide occurrence falls to zero by ~4500 m, despite there being significant areas of land at higher elevations.We suggest that this cut-off in landslide occurrence is owing to the onset of permanent snow, ice and permafrost, all of which become pervasive at elevations of ~5000 m (Chauhan and Thakuri, 2017;Fukui et al., 2007).Glacial, paraglacial and periglacial processes are unlikely to influence landslide occurrence across the sub-regions as there is little coverage of elevations greater than ~5000 m in these areas (Fig. 8b).
Secondly, with respect to hillslope angle, all landslide occurrences typically follow a normal distribution that peaks at 35 -40°.This result is expected given that higher slope gradients are a significant control on slope stability (Hasegawa et al., 2009;Regmi et al., 2013).Within the sub-regions, spatially associated (partial) landslides are more likely than spatially unassociated landslides to occur at slope gradients > 50 o .This is likely because hillslopes above their natural threshold angle are inherently more unstable (Blöthe et al., 2015), so past landslide material that remains on high angle slopes is more likely to fail again or be remobilised in the future.
The third control factor is aspect; all landslides across both regions are most likely to occur on south east, south, and south-west facing hillslopes (Fig. 8e -f).One explanation for this distribution could be due to thermal weathering, whereby south-facing hillslopes experience greater expansion and contraction as a result of high insolation, causing higher rates of rock and regolith damage, and therefore an increase in landsliding (Collins and Stock, 2016).Another potential explanation is the Nepal summer monsoon.The Asia Summer Monsoon typically travels from southeast to northwest (Nayava, 1974), causing more rainfall to be deposited on south and southeast facing slopes.Consequently, these aspects become saturated faster and will experience higher pore water pressures than north facing, more sheltered, hillslopes.This explanation is corroborated by Dahal et al. (2008) who also found increased landsliding on SE -SW facing slopes, and attributed this to the monsoon.In terms of the difference between spatially associated and spatially unassociated landslides, in the sub-regions, we see that spatially associated (partial) landslides are 10% more likely than spatially unassociated landslides to occur on S-facing hillslopes.This implies that spatially associated landslides are more strongly affected by the monsoon (and insolation).This interpretation is logical, as infiltration of rainwater and subsequent increases in pore pressure and decreases in hillslope stability will be exacerbated where failures have occurred in the past, thus leading to more failures in the same location going forward (Iverson, 2000).We would expect this effect to be particularly potent in the years immediately following a past failure, when bedrock and regolith may still be exposed to the elements (i.e. , not yet revegetated) decreasing slope stability (Löbmann et al., 2020).
Finally, when considering the influence of bedrock geology on spatial distributions, landslides of all types and across both study regions are more likely to occur in phyllites, dolomite, and Quaternary conglomerates and sandstones, but less likely to occur in granites.This is unsurprising, given that landslides would be expected more frequently where hillslope strength is lower, and that Quaternary sandstones and phyllites in Nepal are estimated to have Uniaxial Compressive Stresses (UCS) of ~ 27 MPa (Tamrakar et al., 2007) and 46 MPa (Samadhiya and Jain, 2003), respectively.Whereas granites are estimated to have higher UCS values of ~ 58 MPa (Tandon and Gupta, 2015).In addition, spatially associated (partial) landslides in both the main and sub-regions are more likely than spatially unassociated landslides to occur in phyllites.
In terms of size distributions, it is observed that the main region average landslide size is larger than it is for the sub-regions.This is likely a direct result of inventory resolution, with the higher resolution sub-region results being a more reliable reflection of true average landslide size.However, for both the main region inventory and combined sub-region inventories, we see that spatially associated (partial) landslides have larger mean areas and rollover values, but lower power-law scaling exponents.The higher rollover values for the spatially associated landslides indicates that the most frequent landslide area, and the area below which the frequency density of landslides begins to decrease, is greater.
Likewise, the lower power-law exponent for the spatially associated landslides indicates that larger area landslides contribute more to the overall distribution than they do for spatially unassociated landslides.
This result is similar to that found by Samia et al. (2017) in Collazzone, Italy, where it was found that there was a decrease in landslide size with decreasing degree of spatial association.
There are two potential explanations for this observation.One, that it reflects a physical tendency for spatially associated landslides to be larger.As posited by Samia et al. (2017) this could be due to boundary effects, whereby the edges of a past landslide have high instabilities due to over-steepening, decreased cohesion and shear strength, and increased infiltration, which make subsequent landslides in those locations larger.This explanation is supported by our control factor analysis, which shows that spatially associated (partial) landslides are more likely than spatially unassociated landslides at slope angles > 50 o .Alternatively, this could be an artefact of landslide mapping procedure.Landslides are identified via the visual comparison of pre-and post-time slice imagery.Typically, the easiest landslides to observe and map are those that replace vegetation with fresh bare earth, as these have high reflectivity differences that is easily picked up by optical imagery.However, spatially associated landslides are by definition occurring fully within or overlapping with a past landslide that may or not have been revegetated.If the past landslide has not been revegetated, then new landslides overlapping it may be difficult to identify, particularly if they are small, as the reflective difference between weathered bare earth and fresh bare earth is less evident than the difference between vegetation and fresh bare earth.
As such, smaller spatially associated landslides may be under sampled within the inventory compared to similarly small spatially unassociated landslides, an artefact that will be more pronounced in the lower resolution main region inventory.
Finally, differences in average roundness between the main and sub-regions is likely a resolution artefact.On average, the main region contains landslides with higher roundness values, which could be caused by lower inventory resolution making landslide shapes appears rounder.However, despite the differences between the two inventory types, in both cases, spatially unassociated landslides show a higher average roundness than spatially associated landslides.This supports Samia et al. (2017), who showed that spatially associated landslides were typically more elongated.Lower roundness values and thus more elongate landslides are typically associated with a debris flow like typology (Parise and Jibson, 2000), suggesting that spatially associated landslide events are more likely to be reactivated material failing as a debris flow than spatially unassociated events.

What spatial/temporal data are required to observe path dependency?
The second aim of this paper was to assess what spatial and temporal data are required to observe path dependency.This is a useful question to answer as it will allow future studies on this topic to target landslide data collection more efficiently for the purpose of investigating path dependency.
Our results and those of Samia et al. (2017) show that path dependency can be observed across a range of spatial and temporal scales.Here, we used two inventory types, one covering 30 years with annual time slices and landslides mapped at 15 -30 m resolution, and two others covering 10 years with biannual time slices and landslides mapped at 5 m resolution.The fact that path-dependency was observed across both inventory types suggests that only a small number of time-slices (at least 5) and landslides mapped at a low (15 -30 m) resolutions are required to determine whether path dependency is occurring.
However, we also observe that some path dependent results may be scale dependent.For example, we find that the rate of decrease of landslide area with degree of overlap is sensitive to inventory length.
Likewise, whilst the general trends in the overlap index and unaffected area results were similar across all spatial and temporal scales, the magnitudes of these metrics were strongly related to landslide density, which was itself a product of study region size and the type and scale of landslides included in an inventory; specifically, whether large relict deep-seated landslides are included.Thus, higher resolution mapping over smaller spatial scales combined with the inclusion of large deep-seated events will lead to inventories with a higher landslide density, and therefore a greater magnitude overlap index and TUA / AUA relative area.Furthermore, for the overlap index specifically, we conclude that longer inventories with more time-slices will be increasingly accurate, as they allow for greater inter-time slice comparisons, and thus will ensure that any short-term anomalous variations (such as that observed in Fig. 6) do not mask the overall trends in path dependency.
In addition, the results from Section 4.3 highlight that inventory resolution has a strong impact on the apparent influence of landslide control factors.Specifically, figure 7 shows that, when using the lower resolution main inventory data, the spatially associated and spatially unassociated landslides appear to have similar control factor distributions.In contrast, the higher resolution sub-region data suggest that spatially associated and unassociated landslides actually have slightly different control factor distributions.This highlights that whilst lower resolution landslide inventory data can detect overall rates of path-dependency, such data may not be sufficient to undertake a detailed quantification of pathdependency characteristics.Consequently, when investigating path-dependency for the purpose of undertaking landslide susceptibility, where an understanding of control factor is vital (e.g.Samia et al., 2018) the use of higher resolution, and thus more accurate, inventories is recommended.

Conclusions
The main aims of this paper were twofold.Firstly, to use the methods of Samia et al. (2017) to quantify if, and if so how, landslide path dependency occurs in the Nepal Himalayas, and thus provide evidence that the path dependency observed in Italy is also present in other geomorphic regions.As part of this aim, we also investigated whether path dependent (spatially associated) landslides have different size and spatial distributions to spatially unassociated landslides.Secondly, to provide insight into the spatial and temporal resolution of data required to observe and quantify path-dependency.
We find that there is strong evidence for landslide path dependency in Nepal, with all Nepal inventories having more overlaps between earlier and new landslides than would be expected from a random distribution.Furthermore, the unaffected area results show that the actual area of the study regions unaffected by landsliding is lower than the theoretical area that would be unaffected if there was no landslide overlapping, whilst the overlap index tentatively suggests that the amount of overlap between landslides in two time-slices decreases with time.This corroborates the results of Samia et al. (2017), suggesting that their recommendation for this process to be included in landslide susceptibility modelling should apply to geomorphic regions outside of Italy.
We also conclude that spatially associated landslides have larger size distributions than spatially unassociated landslides, and are more likely at slope angles > 50 o and south-facing aspects.This further supports the idea that path-dependency and path-dependent landslides should be treated more explicitly within landslide susceptibility models.
Finally, by investigating path dependency using inventories with two different spatial and temporal resolution landslide inventories, we confirm that path dependency can be observed with 2-year long time slices over 10-years and 30 m resolution landslide data.However, we also find that the rates and magnitudes of path dependent metrics are sensitive to inventory length, study region size and the size/type of landslides mapped, so these must be considered when interpreting path-dependent results.

Figure 1
Figure 1.a) Regional locations of Nepal and our main study region.b) the specific locations of our

Figure 2 .Figure 6 .Figure 7 .
Figure 2. Examples of the spatial associations between landslides in different time-slices using google Earth imagery (GoogleEarth, 2020).Note, the landslide polygons presented in this figure are used to