Do landslides follow landslides? Insights in path dependency from a multi-temporal landslide inventory

Landslides are a major category of natural disasters, causing loss of lives, livelihoods and property. The critical roles played by triggering (such as extreme rainfall and earthquakes), and intrinsic factors (such as slope steepness, soil properties and lithology) have previously successfully been recognized and quantified using a variety of qualitative, quantitative and hybrid methods in a wide range of study sites. However, available data typically do not allow to investigate the effect that earlier landslides have on intrinsic factors and hence on follow-up landslides. Therefore, existing methods cannot account for the potentially complex susceptibility changes caused by landslide events. In this study, we used a substantially different alternative approach to shed light on the potential effect of earlier landslides using a multi-temporal dataset of landslide occurrence containing 17 time slices. Spatial overlap and the time interval between landslides play key roles in our work. We quantified the degree to which landslides preferentially occur in locations where landslides occurred previously, how long such an effect is noticeable, and how landslides are spatially associated over time. We also investigated whether overlap with previous landslides causes differences in landslide geometric properties. We found that overlap among landslides demonstrates a clear legacy effect (path dependency) that has influence on the landslide affected area. Landslides appear to cause greater susceptibility for follow-up landslides over a period of about 10 years. Follow-up landslides are on average larger and rounder than landslides that do not follow earlier slides. The effect of earlier landslides on follow-up landslides has implications for understanding of the landslides evolution and the assessment of landslide susceptibility.


Introduction
Existing landslide research recognizes the critical role that is played by external triggers (e.g. extreme rainfall events, earthquakes and human interferences) along with intrinsic attributes of site (e.g. slope and lithology) that contribute to landslide occurrence Crozier 1986;Guzzetti et al. 2008). Much work has been done on the modelling of landslide susceptibility and hazard. By definition, landslide susceptibility is a nontemporal concept that refers to locations where landslides preferentially occur (Guzzetti et al. 2005;Varnes 1984), whereas landslide hazard describes the likelihood of landslide occurrence in time and space Varnes 1984) along with the magnitude of landslide occurrence . Landslide susceptibility and hazard have been studied using qualitative (Barredo et al. 2000;Ruff and Czurda 2008;Van Westen et al. 2003) and quantitative (Godt et al. 2008;Guzzetti et al. 2005;Lan et al. 2004;Remondo et al. 2003;Van Westen et al. 1997;Van Westen and Terlien 1996;Yeon et al. 2010) approaches.
Qualitative approaches emphasize the role of experience and expert knowledge in determining landslide susceptibility (Van Westen et al. 2003). Quantitative approaches assume that conditions that lead to landslide occurrence in the past and present are likely to cause landslides in the future, thus, the probability of occurrence of future landslides is determined using correlations among various conditioning factors and landslide inventories by statistical methods (Tien Bui et al. 2016;Van Westen and Terlien 1996). Deterministic quantitative approaches use detailed geotechnical and hydrological data in combination with statistical models to estimate the probability of slope failure (Aleotti and Chowdhury 1999;Van Westen and Terlien 1996).
Empirical landslide inventories that document the location and sometimes the date, shape and type of landslides play an important role in assessing landslide susceptibility and hazard Guzzetti et al. 2005Guzzetti et al. , 2012van Westen et al. 2006van Westen et al. , 2008. Most of the landslide inventories are prepared through the interpretation of stereoscopic aerial photographs (Guzzetti et al. 2012). In addition, field mapping (Brunsden 1993), analysis of surface morphology through highresolution DEMs and interpretation of satellite images (Guzzetti et al. 2012) are used for mapping landslide inventories. Overall statistics of landslide inventories have been studied, and it was documented that the frequency-area distribution of landslides follows a power-law distribution for medium and large landslides with an exponential rollover for small landslides (Guzzetti et al. , 2009Malamud et al. 2004;Stark and Hovius 2001;Turcotte et al. 2006;Wood et al. 2015).
Although a vast amount of research has been done on the prediction of landslide occurrence, there has, to our knowledge, been no empirical attention for the effect of earlier landslides on future landslides. This effect is nonetheless expected because landslides typically change the surface morphology (Schuster and Highland 2003), the sediment or regolith properties (Chen 2009), the vegetation (Singh et al. 2014) and the slope angle (van Westen et al. 2006), which are all factors that change landslide susceptibility. If true, such importance of landslide history for landslides susceptibility would be a form of path dependency (a concept from complexity theory (Phillips 2006;Temme et al. 2015))-indicating that the history of the landsliding process affects its future through one or more legacy effects. A likely reason for the lack of attention for quantifying the effect of earlier landslides on future landslides is that multi-temporal landslide inventories are very difficult to obtain (Atkinson and Massari 1998;Brenning 2005) and high-resolution multi-temporal datasets of intrinsic properties are virtually absent.
In this paper, we explore the possible effects of earlier landslides on future landslides. For this, we will use a rich multi-temporal landslide inventory from the Collazzone area in central Umbria, Italy (Ardizzone et al. 2013;Galli et al. 2008;Guzzetti et al. 2006).
Our main objective is to investigate whether earlier landslides determine the susceptibility for future landslides, i.e. whether there is path dependency in landslide occurrence. We test two hypotheses. The first hypothesis is that landslides follow landslides; i.e., that landslides increase the likelihood of another (in our terms: follow-up) landslide occurring in the same place. We operationalize this hypothesis without attention for the role of intrinsic attributes (e.g. slope and geology) to offer a clear contrast with existing approaches. The second hypothesis is that follow-up landslides differ from other landslides in terms of their shape, size and frequency-area distribution.

Study area and data
The Collazzone area extends for 78.9 km 2 in central Umbria, Italy (Fig. 1). Elevation in the area ranges from 145 to 634 m above sea level, and the slope, computed from a 10 × 10 m Digital Terrain Model (DTM), ranges from 0°to 63.7°with a mean value of 9.9°. The terrain is hilly with asymmetrical valleys, with lithology and attitude of bedding controlling the slopes. Only sedimentary rocks crop out in the area and include (1) alluvial deposits, Holocene in age; (2) travertine, Pleistocene in age; (3) continental deposits (gravel, sand and clay), Plio-Pleistocene in age ); (4) layered sandstone and marl, Miocene in age and (5) thinly layered limestone, Lias to Oligocene in age. The land use is characterized by cropland, forests, urban areas, pastures, vineyards, orchards and water. Farming in the area favours the development of slope failures and erosion. Soils range in thickness from a few decimetres to more than 1 m; they have a fine or medium texture . Soils have a xeric moisture regime characterized by cold and moist winters and dry summers. Precipitation is most abundant in October and November; with a mean annual rainfall of 841 mm in the period from 1951 to 2013. Snow falls in the area on average every 2-3 years. Landslides are abundant and range in age, type, morphology and volume from relict-partly eroded-large and deep-seated landslides, to young, mostly shallow landslides involving the soil mantle. Landslides are triggered predominantly by meteorological events, including intense and prolonged rainfall and rapid snowmelt.
Description of the multi-temporal landslide inventory A detailed multi-temporal landslide inventory is available for the Collazzone area ( Fig. 1 and Table 1). The inventory was originally prepared at 1:10,000 scale through the visual interpretation of five sets of stereoscopic aerial photographs taken unsystematically in the period 1941-1997 at scales ranging from 1:13,000 to 1:33,000. The landslide inventory was continued in the period from 1999 to December 2005 through field surveys carried out after periods of prolonged rainfall and in March and May 2010 using stereo satellite images (Ardizzone et al. 2013).
Landslides ages were estimated from the date of the aerial photographs and the morphological appearance of the landslide. In each of the five sets of stereoscopic aerial photographs used to prepare the multi-temporal inventory, landslides that appeared 'fresh' on the aerial photographs were separated from the other landslides that had occurred since the previous photograph. The date (i.e. year) of the aerial photographs used to identify the landslides was assigned to the 'fresh' slope failures. The other slope failures (i.e. the 'non-fresh' landslides) were attributed to the period between the date of the aerial photograph where they were identified and the date of the previous aerial photographs. Galli et al. (2008) demonstrated the high quality of the multitemporal landslide inventory map for the Collazzone study area. Nonetheless, we acknowledge that the mapping and age attribution are not perfect. Errors are probably largest in time slices that were obtained through field mapping due to the difficulty in translating lateral assessments onto vertical photographs, but even here we expect the uncertainty in determining a landslide boundary to be less than 3 m. For other time slices, the aerial photographs and satellite images that were used have resolutions of 1 m or better, and hence, uncertainty is of the same order. With an average size of landslides of over 5000 m 2 , it was expected that the uncertainty in mapping did not substantially affect our results.

Methods
In addition to our use of the term path dependency (indicating that previous process activity affects future process activity through legacy effects), it is useful to clarify some other terms that are important for our analysis. The mono-temporal inventories that together make up our multi-temporal inventory are separately called 'time-slices', whereas subdivisions of the multi-temporal inventory based on topology (see below) are called 'sub-inventories'. We use 'earlier' and 'follow-up' to describe landslides based on their relative order in the multi-temporal inventory. Note that follow-up landslides are not reactivated landslides. We consider a landslide a reactivated landslide when all or most of the landslide moved down again, under the same general condition as the first landslide. Instead, follow-up landslides are new landslides that have different size and shape than the pre-existing landslide. The particular topological relations of interest between landslides in earlier and follow-up time slices were called 'spatial association'.
To test our hypotheses, two sets of analyses were done on the multi-temporal landslide inventory. The first set of analyses focuses on the effect of overlap between landslides on the total area affected by landslides ( Fig. 2 (a)), on the degree of overlap between landslides from different time slices (Fig. 2 (b)) and on the number of overlaps between landslides over the entire inventory ( Fig. 2 (c)). The second set of analyses focuses on the properties of classes of landslides that vary in their spatial association with landslides from the earlier periods ( Fig. 2 (d)).

Degree of overlap
We propose three complementary methods for this purpose to express the importance of overlap between landslides: (1) unaffected area, (2) overlap index and (3) number of overlaps (Fig. 2).
The first method measures the cumulative effect that landslide overlapping has on the total area affected by landslides, by comparing a theoretically unaffected area with the actually unaffected area. This method allows assessment of the amount of reduction in affected area caused by (potentially followed-up) overlapping landslides. The theoretically unaffected relative area (−, where − denotes that this is a dimensionless parameter) if no overlap between landslides would occur is defined as:  where AS is the area of interest (m 2 ) and AL Ti is the total area (m 2 ) of all landslides in the time slices i of an inventory. The relative area that is actually unaffected by landslides over the entire period is calculated as: where ∪ denotes the spatial union among landslides in consecutive time slices, i.e. the total area covered by landslides after accounting for overlaps.
The second method provides a standardized measure of overlap between landslides in consecutive time slices. This overlap index (−), earlier termed relative area overlap (Maruca and Jacquez 2002), is calculated as: where t is the assumed date of occurrence of landslides from a first time slice (see below), and t + n is the average date of occurrence of landslides in a later time slice, and ∩ denotes the geometric intersection (i.e. the overlapping area) between two time slices. This is qualitatively similar to the error index that has been used to calculate the positional mismatch for pairs of corresponding landslides in the two inventories (Ardizzone et al. 2002;Carrara et al. 1992;Santangelo et al. 2015). The overlap index is non dimensional, ranges from zero (no overlap) to unity (perfect overlap) and is not a function of the size of the study area. The overlap index was calculated for consecutive time slices and for pairs of time slices that are two or three time slices apart. The value for the overlap index was then related to the time that passed between time slices. As explained above, in some time slices, there is more uncertainty about the date of occurrence of a landslide than in others. Uncertainty is larger in the time slices that describe longer periods, such as the time slice describing landslides that occurred between 1954 and 1977 (see Table 1). In some other cases, the date of occurrence is relatively well known, for instance in the time slice made of landslides that occurred after the 01-01-1997 rapid snowmelt event. These uncertainties, large or small, propagate into uncertainty about the time elapsed between landslides. We assumed uniform probability distributions of landslides in the time slices that describe periods-meaning that landslides may have happened at any moment during the considered time slice, with equal probability. This is despite the fact that there may have been significant rain, snowfall, and snowmelt events in these periods that may have caused landslides to be clustered in time. No information about this intra-period variability was available. Using the uniform probability distributions, we performed stochastic simulation, randomly placing 10,000 pairs of landslides in the periods of two time slices and then recording the time passed between each pair. The median time passed between time slices (t − t + n) was used as the time passed between time slices, and the first and third quartile were recorded to express uncertainty of the time between both time slices. This approach is unbiased and therefore no bias in consequent analyses is expected.
The third method quantifies the relation between area and number of overlaps. This method was used to quantify repeated overlapping over multiple time slices of the multi-temporal inventory. Every time slice of the inventory was rasterized at 1 m resolution, and the cells involved in landsliding in each inventory were given a value of 1. Then, all raster were summed, resulting in a raster where the values correspond to the number of overlaps over the entire multi temporal inventory. The analysis was performed using standard raster operations in ArcGIS. This measure does not change with the order of time slices, but is affected by the area involved in landsliding in the multi-temporal inventory and by the number of time slices in it. Therefore, it was compared with the number of overlaps from a null model in which the observed area of landsliding per time slice was assigned randomly to the study area, regardless of slope or lithology. Since spatial information was Fig. 2 Overview of the analyses performed on the multi-temporal landslide inventory to test hypothesis 1 (a, b, c) and hypothesis 2 (d)

Original Paper
Landslides not needed for this calculation, the null model was calculated in Microsoft Excel.
Effect of different spatial association Spatial associations between landslides from a given time slice and the immediately preceding time slice were used ( Fig. 2  (d)). Four classes of spatial association were defined (Fig. 3). The 'inside' class contains landslides that are completely inside landslides from the earlier time slice. The 'partial' class contains landslides that partially overlap landslides from the earlier time slice. The 'touching' class contains landslides that are outside landslides from the earlier time slice, but that touch landslides from this time slice, and the 'outside' class contains landslides that neither overlap nor touch landslides from the earlier time slice. This classification was chosen based on the assumption that different spatial associations between landslides relate to different mechanisms in which landslides affect the probability of follow-up landslides. Landslides that are 'inside' earlier landslides, for instance, may particularly point to saturation of materials on top of an earlier landslide (Igwe and Fukuoka 2015), and landslides that are 'touching' an earlier landslide may point to a slopeeffect because slopes are most strongly changed along the borders of earlier landslides (van Westen et al. 2006). The 'outside' class was considered likely unrelated to earlier landslides, by any mechanism.
Standard GIS tools were used to separate landslides from each time slice of the multi-temporal inventory, except the first time slice, into the four classes of association. The subclasses of association (sub-inventories) were then merged over the different time slices. Then, the following geometrical properties were calculated for landslides in each merged sub-inventory: minimum, mean and maximum area, roundness (a measure of shape) and the three parameters of the inverse gamma distribution that are often used in landslide research to describe frequency-area relationships (Malamud et al. 2004), see below. This allowed us to test our second hypothesis that landslides following up on earlier landslides have different geometric properties.
Minimum, mean, median and maximum sizes of landslides were calculated with standard GIS tools. Roundness (−) is introduced here as a simple measure of shape: where the theoretical circular perimeter (m) is the perimeter of the landslide, if it would have been perfectly round with the same area: where A L is the area of the landslide (m 2 ). Roundness values closer to unity indicate round shapes, values close to zero indicate more elongated shapes (Fig. 4). In earlier research, a landslide geometry generating algorithm was used to approximate the shape of landslides using geometric features such as upper and lower length and lateral boundary of the landslide scar (Chiang 2015). Taylor et al. (2015) approximated landslide shapes by ellipses and then used the length to width ratio of these ellipses to characterize shape. We argue that the proposed roundness is a simpler measure that Fig. 3 Examples of spatial association of landslides with landslides from the previous time slice Landslides makes no assumptions about landslide shape that nonetheless captures the most important difference between landslide shapes.
To estimate power-law parameters, we fitted the threeparameter inverse gamma distribution (Malamud et al. 2004) to the four sub-inventories, using LANDSTAT software (version 9) (Rossi and Malamud 2014). The three-parameter inverse gamma distribution is: Parameter α controls the steepness of the right tail of the probability density function. Parameter η controls the steepness of the left tail of the probability density function. Parameter λ controls the position of the rollover. From these parameters, the rollover, the most likely size of a landslide, can be calculated. Maximum likelihood estimation was used and uncertainties for the three parameters and the rollover were calculated using bootstrapping with 250 repetitions.
We used standard Analysis of Variance to test whether differences in landslide properties between sub-inventories are significant, complemented with T-tests for normally distributed properties (Student 1908) and U-tests for non-normally distributed properties (Mann and Whitney 1947).

Degree of overlap
Over time, significant differences occur between the theoretically unaffected area (Eq. 1) and the actually unaffected area (Eq. 2) (Fig. 5). After the last time slice in the inventory, from May 2010, 7.8 % of the area has not experienced landsliding during the entire period of observation because of overlapping landslides. Clearly, there is a significant amount of overlap between landslides.
The overlap index (Eq. 3) appears to be negatively correlated with the time passed between time slices, although variation in overlap index is large, especially when observations are close together (Fig. 6). The overlap of landslides with those from an earlier time slice appears to decrease substantially over a period of about 10 years. This suggests that the cause of new landslides is not completely external. If the cause of overlap was completely external (for instance through repeated landsliding at particularly dangerous locations), then no relation between overlap index and the amount of time passed between slides should exist.
Instead, earlier landslides themselves appear to affect the probability of reoccurrence to a substantial degree-probably through legacy effects. Importantly, this effect decreases over time. This means that if a landslide happened longer ago, it is less likely to be overlapped by a future landslide. However, this is not always the case: some time slices that are close together in time, have low overlap index. Apparently, not all landslides affect their environment such that there is a larger probability for future landslides.
About 200 m 2 of the study area has been affected seven times by landslides in the 17 time slices of the multi-temporal inventory (Fig. 7). In the null model in which landslide cells are placed randomly in the study area, seven overlaps do not occur. Overall, much more of the study area has been affected by three or more overlapping landslides than in the random model, whereas the area where no overlapping occurs, or where only one or two overlaps occur, is less than in the random model. The relation between number of overlaps and the logarithm of area is about twice as steep in the random model than in the dataset, again a clear sign of spatial overlap between follow-up landslides.
The mean size of the landslides with spatial association 'inside', 'partial' and 'touching' is larger than those with spatial association 'outside' (Table 2, p < 0.001). The subtle decrease in mean area from 'inside' to 'partial' and 'touching' is not statistically significant (p > 0.05).

Shape factor
Landslides that are not spatially associated with landslides from the immediately preceding time slice appear to be more round than other landslides (Fig. 8). The type of spatial association is a significant determinant of shape (Table 2, p < 0.001). This is because of the higher roundness of the 'outside' type-all other types of association are not significantly different from each other (p > 0.05).

Original Paper
Landslides Frequency-area statistics The frequency-area distributions for the four sub-inventories representing different spatial association exhibit power-law scaling for large and medium landslides with a rollover for small landslides ( Fig. 9a-d). The three parameters of the inverse gamma distribution, and the resulting rollover, differ between subinventories (Table 2). When the strength of spatial association between landslides decreases (from 'inside' to 'partial' and 'touching'), the exponent of the inverse power-law (α) increases. The average value of (α) in the spatially associated sub-inventories is greater than in the 'outside' spatially un-associated sub-inventory. Interestingly, the value of the rollover is much smaller for the unassociated 'outside' sub-inventory than for the three other subinventories. The rollover also decreases with decreasing strength of spatial association. All differences between sub-inventories were significantly different (p < 0.001).

Discussion
In this discussion, we will first focus on our findings relating to the two hypotheses, before considering wider implications.

Do landslides follow landslides?
We found that the multi-temporal landslide inventory of our study site recorded substantial overlap between landslides and an associated reduction in the fraction of the study area affected by landslides (Fig. 5). The overlap is more for landslides happening sooner after an earlier landslide and decreases over a timescale of about 10 years (Fig. 6). This indicates that landslides in our study site do follow landslides. Importantly, it also suggests an internal control on the landslide system next to the range of intrinsic controls such as lithology, slope steepness and vegetation, which are commonly the focus of landslide susceptibility studies (Guzzetti et al. 2005;Van Westen et al. 1997). We propose that the mechanisms that can explain the nature of this internal control (i.e. a positive landslide-landslide effect that decreases over time) fit into two categories. First, the probability of landslides overlapping earlier landslides may decrease over time because the deposits of an earlier landslide stabilize due to, e.g. a more stable slope geometry, the regrowth of vegetation and the repair of soil structure and cohesion. Second, the probability can decrease over time because follow-up landslides have already occurred, erasing topographic instabilities that remained after an earlier landslide. Targeted field observations with high temporal resolution could quantify the relative importance of these two reasons by focussing on vegetation regrowth and detailed topographic evolution.
Regardless of the mechanism, our findings suggest that there is path dependency among landslides: older landslides act in some way or other as initiators for follow-up landslides for a certain period. This is also consistent with the number of landslide Fig. 6 Overlap index as a function of time between time slices of the multi-temporal landslide inventory; the data shows a higher overlap index for time intervals up to approximately 10 years Landslides overlaps over the entire period captured by the multi-temporal inventory (Fig. 7). This number of overlaps is larger than in the null model where landslides occur randomly in the study area-although when considered separately from Fig. 6, this observation can also be explained by spatial differences in landslide susceptibility. Such landslide self-organization  into emergent patterns is not captured in traditional cause-effect studies of landslide susceptibility. It also adds an important consideration to the discussion about landslide selforganized critical behaviour Hergarten 2003;Turcotte et al. 2002). Apparently, not all landsliding potential ('metastable regions', Guzzetti et al. 2002, p171) is removed by a landslide-instead first landslides appear in some cases to increase the potential for follow-up landslides. The interplay between the traditional cause-effect approach and path dependency in a self-organization approach should be explored and quantified in follow-up studies.
The period over which the positive landslide-landslide effect is observed-about 10 years in our study area-is very likely specific for the local settings and other areas with similar climate, geology, topography, soils and vegetation. We maintain that this period is a measure of landslide path dependency that reflects the rate of processes that reduce landslide susceptibility after a first landslide, or the rate at which follow-up landslides occur in our study site (see above). In areas with more frequent landsliding, or in sites with more rapid vegetation growth that restores stability, it may be shorter (all other factors being equal).

Are follow-up landslides different?
We also hypothesised that follow-up landslides have different properties in terms of size, shape and frequency-area statistics. This is to some extent the case. Mainly, the mean size of landslides is lower for landslides that are not spatially associated with earlier landslides, than for those that are. Apparently, on average, followup landslides are larger (although there are also less of them, Table 2). In addition, our empirical data suggest that there may be a slight decrease of size with decreasing strength of spatial association between landslides (from 'inside' to 'partial' to 'touching' landslides).
The larger average size of follow-up landslides may be explained by conditions that have changed after occurrence of an earlier landslide, leading to increased landslide susceptibility. These changed conditions may relate to changes in surface morphology (Schuster and Highland 2003), sediment properties (Chen 2009), the vegetation (Singh et al. 2014), slope angle and land use (van Westen et al. 2006).
The shape of follow-up landslides differs from the shape of other landslides ( Table 2). Landslides that are spatially associated with earlier landslides are less round (on average) than those that are not (p < 0.01). The difference in roundness can be explained by two mechanisms. First, weaker materials in an earlier landslide may move faster, and therefore further than other materials. This would then be reflected in less round shape. Second, especially for landslides of the 'touching' class, an earlier landslide may be higher in the landscape, deflecting landslides downslope, which leads to them being more elongated.
Both spatially associated landslides ('inside', 'partial' and 'touching' sub-inventories) and spatially un-associated landslides (outside sub-inventory) follow a power-law scaling in the right tail of inverse gamma distribution for large and medium landslides with an exponential rollover in the left tail of distribution. However, the exponent of the inverse power-law (α) and rollover (λ) differs between spatially associated landslides (even between sub-inventories) and spatially Fig. 7 Area versus number of overlapping landslides. The area experiencing more than four overlaps is larger for the observed dataset than for the random dataset λ ± S D 5 2 ± 6 6 9 ± 5 9 0 ± 1 1 4 9 ± 1 Rollover ± SD (m2) 1073 ± 157 1317 ± 158 1858 ± 241 670 ± 25 un-associated landslides (Fig. 9a-d). The exponent of inverse power-law (α) on average in spatially associated sub-inventories (i.e. 'inside', 'partial' and 'touching') is greater than un-associated sub-inventory ('outside'). For landslides that are spatially associated with earlier landslides, the exponent of power-law (α) increases when the strength of the spatial association (from 'inside' to 'touching') decreases ( Table 2). The exponent of the inverse power-law (α) for 'touching' landslides (α = 1.78) seems consistent with the most reported range of values in the literature (1.5 < α < 2.5) (Borgomeo et al. 2014;Guzzetti et al. 2002;Malamud et al. 2004;Van Den Eeckhaut et al. 2007). The other landslides (i.e. 'inside', 'partial' and 'outside') exhibited lower values for exponent of the inverse power-law (α) (i.e. 1 < α < 1.5) ( Table 2). The small exponent of the power-law in the sub-inventories suggests that larger landslides are contributing to each sub-inventory (Borgomeo et al. 2014;Van Den Eeckhaut et al. 2007). The rollover (λ), which represents the size of the most frequent landslide, is larger for spatially associated landslides than for non-spatially associated landslides ('outside', Table 2). Interestingly, the rollover for the 'touching' landslides is three times larger than for the 'outside' landslides and larger than for the 'inside' and the 'partial' landslides. Analogously to our explanation for lower roundness, this may be explained by a boundary effect that can be caused by steeper slopes on the sides of an earlier landslide. The lateral boundary of a landslide is an intrinsically more disturbed part of a landslide, where the mechanical properties are weaker and infiltration is commonly larger. These factors contribute to instability, and hence to larger landslides. Fig. 8 Roundness as a measure of shape for different classes of spatial association. Error bars indicate the standard error Fig. 9 Comparison of the inverse gamma probability density function fit to frequency densities from four landslide sub-inventories representing different spatial association types Landslides Implications for landslide susceptibility assessment The traditional division between landslide susceptibility and hazard is that the former describes spatial differences in landslide probability, and the latter describes the likelihood of landslide occurrence in time and space Varnes 1984). Our results suggest that a time-related internal control-the time that passed since an earlier landslide occurred-also plays a role in determining how susceptible a location is to landsliding. Apparently, susceptibility is not time-invariant and contains a temporal element that traditionally was seen as part of the definition of hazard. Therefore, spatial and temporal probabilities of landsliding are not independent (Guzzetti et al. 2005). This means that the susceptibility of a location changes over time, as the effects of an earlier landslide slowly disappear. This can be expressed as: where the susceptibility for landsliding of a certain location s at time t is not only a function of conditioning attributes such as slope or lithology, but also of the time passed since an earlier landslide in the same location or close by. In our dataset, the function describing the effect of earlier landslides on susceptibility is positive and decreases over time (e.g. Fig. 6). An earlier landslide appears to make a location more likely to experience landsliding again over a period of about 10 years, after which susceptibility appears to return to its previous value (Fig. 10). We explored possible mechanisms for this effect above.
The effect of the type of spatial association with (or the distance to) an earlier landslide is complex (Table 2). In places, a closer spatial association of a landslide with an earlier landslide leads to properties closer to spatially unassociated landslides (e.g. for the rollover), whereas in other places it leads to properties that are more different from unassociated landslides (e.g. the mean area) ( Table 2). This complexity should be further researched since it may lead to a deeper understanding of the processes causing follow-up landsliding.
Other geomorphic settings could lead to different functions describing the effect of earlier landslides. For instance, exhaustion of sediment after landsliding, which occurs in locations where only little soil overlies bedrock (e.g. the steep basalt slopes of the Drakensberg in South-Africa (Singh et al. 2008)), would lead to a decreasing of landslide susceptibility instead of an increase. Given that weathering is required for new soil to form, it would also probably be a longer term effect than the 10 years that we observed (Fig. 6). The exhaustion of soft material apparently plays no important role in our study site-fitting with its geology of predominantly soft, easily weatherable rocks. A complicating factor is that when the duration of the effect of landslides on landslide susceptibility increases much beyond a decade, other slow changes will also start affecting susceptibility, such as land use changes (e.g. deforestation, afforestation and changes in agricultural practises) and climate changes. This makes it more difficult to observe the effect that landslides themselves have on susceptibility.
The decreasing temporal effect of earlier landslides on susceptibility is to some extent comparable with the legacy effect of earthquakes on the occurrence of earthquake-induced landslides (Lin et al. 2007;Marc et al. 2015;Parker et al. 2015). For instance, Marc et al. (2015) demonstrated that landslide rates were significantly elevated within 0.7-4.5 years after four earthquakes with magnitude higher than 6 Richter. After this, rates dropped gradually back to pre-earthquake levels.
Implications for mapping and monitoring of landslide populations Most landslide inventories are not multi-temporal, although they do often describe landslides that occurred over a longer period or in response to several events. From such inventories, it is not possible to estimate the time-dependent effect of landslides on landslide susceptibility. In settings where that effect is strong, this may lead to susceptibility assessments that are biased. In statistical terms, this would be because the landslide observations used to estimate a susceptibility model or map are not independent. In geomorphic terms, if the landslide-landslide effect is positive and decreases over time (Fig. 6), resulting susceptibility maps would probably be biased to a smaller range of conditioning attributes than would be otherwise the case. Some places that are susceptible to landsliding would then possibly be considered safe. In case of a negative landslide-landslide effect, they would probably be biased to a larger range of conditioning attributes. In this case, some places that are safe would possibly be considered susceptible. How substantial these biases should be the subject of further study. Clearly, this calls for more multi-temporal landslide inventories. Fig. 10 A landslide susceptibility map bases on slope units (coloured polygons) with the occurrence of landslides over time (ellipses with black outline). In the susceptibility map, green means low risk, yellow means moderate risk, light red means high risk and dark read means very high risk. The implication of our work is seen in the clustering of landslides after a first landslide happens within a time scale of about 10 years. This is reflected in a temporarily higher susceptibility in the affected slopes. The two temporarily highly susceptible slopes (dark red polygons) experience different landslide dynamics as a result of different landslide histories

Conclusions
In our study area, landslides follow landslides-and they do so more often than expected based on a random control. This landslide path dependency is strongest over a timescale of about 10 years, after which the effect of earlier landslides on follow-up landslides appears to decrease. This legacy effect has a considerable consequence on the relative area that is affected by landslides. Follow-up landslides differ from non-follow-up landsides in terms of size and shape: they are typically larger and less round. Also, the exponent of the power-law (α) and the value of rollover in spatially associated landslides on average are greater than in spatially un-associated landslides. From these findings, we conclude that landslide susceptibility in our study area should be considered as a dynamic measure that reflects path dependency and selforganization in landslides. Susceptibility changes because the effect of previous landslides on susceptibility, which in our case is positive, decreases over a period of about 10 years. Although our research concentrated on the Collazone region in Italy, we expect these results to be relevant for other landslide prone areas as well.