Origin of the power-law exponent in the landslide frequency-size distribution

Landslide statistics is characterized by a power-law frequency-size distribution (FSD) with power exponent α centered on 2.2-2.4, independently of the landslide trigger. So far, the origin of the α-value, critical to probabilistic hazard assessment, remains hypothetical. We present a generic landslide cellular automaton (LSgCA) based on the rules of Self 10 Organized Criticality and on the Factor of Safety (FS) concept. We show that it reproduces the power-law FSD for realistic parameter ranges (i.e. cohesion, soil unit weight, soil thickness, angle of friction, slope angle, pore water pressure) with LSgCA simulations yielding α = 2.17±0.49, which is in agreement with α = 2.21±0.53 obtained from an updated metaanalysis of the landslide literature. The parameter α remains stable despite changes in the landslide triggering process, with the trigger only influencing the spatial extent of the landslide initiation phase defined from an FS contour. Furthermore, 15 different FS formulations do not significantly alter the results. We find that α is constrained during the initiation phase of the landslide by the fractal properties of the topography, as we observed a positive correlation between fractal dimension and α while α did not change during the propagation phase of the LSgCA. Our results thus suggest that α can be estimated directly from the FS map for probabilistic landslide hazard assessment. However full modeling (including the propagation phase) would be required to combine the spatial distributions of landslide and exposure in probabilistic risk analysis. 20

Of the sheer number of strategies to model landslides (e.g., Densmore et al., 1998;Hergarten and Neugebauer, 1998), cellular automata (CA) that follow some of the rules of Self-Organized Criticality (SOC) provide a simple and natural approach to generate power-laws.Indeed, SOC, epitomised by the Bak-Tang-Wiesenfeld CA (Bak et al., 1987) and generalized by the Olami-Feder-Christensen CA (Olami et al., 1992), is characterized by the emergence of a power-law FSD due to a simple bottom-up triggering process analogue to an avalanche.Although SOC CAs and other CA variants are broadly used for landslide modelling (e.g., D'Ambrosio et al., 2006), the sandpile model (Bak et al., 1987), with slope α = 1.0-1.2does not explain the steeper slope observed for landslides (Hergarten and Neugebauer, 2000;Turcotte et al., 2002).Pelletier et al., (1997) retrieved α = 2.6±0.1 from a percolation model controlled by a threshold shear stress dependent on a slope based on physical parameters, but only considered the area of the landslide initiation phase, not of the landslide itself.
Interestingly, Hergarten and Neugebauer (2000) obtained α = 2.1 by applying a two-variable product in a SOC system with one relaxation variable and one time-dependent weakening variable, which was related to the factor of safety as an ad-hoc assumption in Hergarten (2013).Piegari et al., (2006;2009) only obtained reasonable α estimates by allowing arbitrary parameter variations to match the data (see discussion in Hergarten, 2013).It goes the same with Guthrie et al., (2008), who calibrated their input parameters to fit an observed landslide FSD.Those SOC models (Hergarten and Neugebauer, 2000;Piegari et al., 2006;2009), but also other models (e.g., Densmore et al., 1998), neglected changes in topography by considering only one individual slope, questioning their validity for landslide inventories whose events take place on complex topography.Moreover, none of those models lead to landslides that seem very realistic (Hergarten, 2013).Hergarten (2012) proposed a simple model inspired from the basic models of SOC where avalanches propagate on the surface when a certain threshold is exceeded (Bak et al., 1987;Olami et al., 1992) and applied it on real topographies.It was however limited to rock falls and did not define any physical trigger ("random impacts" were used instead unclear why the landslide FSD seems to be independent of the triggering mechanism and what process or specific physical parameter explains α = 2.3±0.6.Finally, the available models remain mostly abstract and difficult to implement in probabilistic multi-hazard assessment (as in, e.g., Mignan et al., 2014;Liu et al., 2015).
The aim of this article is to present a generic landslide cellular automaton (in short, LSgCA) that is consistent with the broad range of α-values observed in Nature and is realistic enough for inclusion in future generic multi-risk (GenMR) analyses (e.g.Mignan et al., 2014;2017;2018;Komendantova et al., 2014;Liu et al., 2015;Matos et al., 2015).The purpose of the LSgCA is not to replace more sophisticated existing landslide modelling tools but to provide a transparent and simple, yet robust, approach to understand what parameter has the greatest influence on α.This is of importance in the probabilistic hazard assessment of landslides to be able to extrapolate the size of larger potentially damaging events (e.g., Guzzetti et al., 2005).We first present the LSgCA approach, where the initiation phase of landslides, the propagation phase and topography simulation are described (section 2).We then update the meta-analysis originally made by Van Den Eeckhaut et al., (2007) by adding about 40 new cases, hence more than doubling the number of data points available to estimate α (section 3.1).
Finally we investigate at what stage of the LSgCA process and for what parameterizations the α distribution observed in Nature is retrieved (section 3.2).
An FS map can then be produced that depends mainly on the spatial distribution of θ, i.e., on the topography.The impact of earthquake shaking can be modelled via the concept of Newmark displacement (Newmark, 1965): function of the critical acceleration ratio a c /a max and earthquake magnitude M, with critical acceleration  !=  − 1 sin and g = 9.81 m/s 2 (Jibson, 2007).Landslides are initiated at locations that exceed D N = 15 cm (Jibson et al., 2000).The earthquake peak ground acceleration a max is derived from empirical relationships, here using As investigated in section 3, we do not expect the use of different FS models, nor different triggers, to have an impact on the landslide FSD α-value since they only impact the spatial extent of the landslide initiation zones defined by the thresholds FS ≤ 1 or D N ≥ 15 cm.Indeed, a power law behaviour can emerge whether at the initiation phase by a percolation process on a fractal lattice (e.g., Pelletier et al., 1997) or during the propagation phase in a SOC-style process (see review by Mignan (2011) for the case of earthquake power laws).In the first case, α depends on the fractal dimension of the lattice, in the second, on the propagation rules that depend on soil conditions, and not on the triggering mechanism, at least not at the lowest-order.

Propagation phase (post-failure motion)
The proposed CA is based on the concept of sandpile (Bak et al., 1987) and defined for a square grid composed of cells (x,y) binned in Δs increments.Variables are the altitude z(x,y) and soil depth h(x,y).Input parameters are the initial topography (z, h) and soil conditions (φ, C, γ, etc.) from which FS is computed (Tables 1-2).We use the Moore neighbourhood nomenclature (Gray and New, 2003), and compute the maximum slope θ max in the direction dir Moore (θ max ).In contrast to the SOC random field where triggering occurs at random cells, triggering is here initiated in cells for which conditions FS ≤ 1 or D N ≥ 15 cm are satisfied (section 2.1).The landslide footprint is defined as LS RS (x,y | FS ≤ 1) = 1 for rainstorm triggers and as LS RS (x,y | D N ≤ 15 cm) = 1 for earthquake triggers, and with LS(x,y) = 0 elsewhere.
For LS(x,y) = 1 and h(x,y) > 0, the propagation rule is where the mass movement is defined by The mass is progressively transferred from cells to cell by Δh such that a stable slope θ stable , if h is high enough, is reached.If h is too low, i.e., h < Δh, a scarp forms equivalent to h = 0 once the propagation rule is applied.θ stable is the angle for which FS tends to 1 from above for rainstorms or when D N tends to 15 cm from below for earthquakes.For each modified cell dir Moore (θ max ), LS(dir Moore (θ max )) = 1.The process continues until the above condition fails at all cells or when the number of unstable cells reaches a plateau (i.e., stability criterion to avoid infinite loops).It should be noted that the process is non-SOC since it naturally dies off after any event.However, the constant load defined from successive earthquakes and/or rainstorms might lead to a SOC-like behaviour on the long-term (Hergarten, 2003).
Figure 2 shows a landslide generated by the LSgCA in a smooth virtual region, as defined in Komendantova et al., (2014), Liu et al., (2015) and Mignan et al., (2017) for generic multi-hazard testing, here with an earthquake as trigger (black segment).Note that the propagation phase, in itself, does not lead to a power-law distribution (nor does the initiation phase for the matter, as explained in the previous section).In contrast to SOC where the spatial field is random and the load increased by random increments, the proposed CA can apply on smooth surfaces and produce relatively realistic landslide footprints.Loading, by an increase or addition of an elliptic instability patch (for the earthquake trigger), may yield only one landslide, as shown in Fig. 2 (see video in the supplementary materials).
For the rest of this study, we will apply the LSgCA on more realistic topographies and investigate whether the α distribution observed in Nature is explained from the topography itself (i.e., percolation model) and whether it is altered by the landslide propagation phase (i.e., SOC-like model).We model the fractal topography using the diamond-square algorithm (Fournier et al., 1982), based on fractional Brownian motion with σ n = 2 -nH δ, δ being the z-deviation at first iteration, n the iteration number, and H the Hurst exponent.The fractal dimension is simply D f = 3-H.We will test 2.1 ≤ D f ≤ 2.9, representing increased terrain roughness with 7 iterations leading to a 2 7 +1 = 129 × 129 cells grid with the constant thickness h.The resulting topography is considered "pristine" and is unrealistic due to its many steep peaks (Miller, 1986) and to the constant h.However, we can then apply the LSgCA to "erode" the topography, which yields a new topography that is more realistic (i.e., scarps become devoid of soil while basins contain more soil).We define the topography on a 10 by 10-km grid with a 10-meter resolution.Individual landslides are defined as continuous patches verifying LS(x,y) = 1.The α parameter is then calculated using Maximum Likelihood Estimation (MLE), as described in Clauset et al., (2009).

Updated meta-analysis
To validate our model based on FSD statistics (section 3.2), we first update the review made by Van Den Eeckhaut et al., (2007).We do so by adding observations made during the last decade.Inputs are listed in Table 3 and the resulting distribution plotted in Fig. 3.This new meta-analysis includes a total of 70 cases, including 40 new compared to the 2007 study.Note that α = α cum +1 when the FSD is cumulative but also that α = 1.5αV when the landslide size is volume V instead of area A (with V ~ A 3/2 ).When compared to other models, α = α 2 +1 = ρ+1 with α 2 from the double Pareto distribution (Stark and Hovious, 2001) and ρ from the inverse Gamma distribution (Malamud et al., 2004). 1 (Harp and Jibson, 1995;Xu et al., 2015;Bucknam et al., 2001) are also plotted.

LSgCA application on fractal topography
We ran simulations in which the ground parameters were drawn randomly from the range of values presented in Table 2.
The four FS formulations of Table 1 were systematically tested for sensitivity analysis.We first evaluated α at different steps of the landslide process (with constant D f = 2.4): (1) at the initiation phase, directly based on the FS map, assuming that the landslide area limits itself to FS < 1 patches.This is equivalent to a percolation model where a FSD is defined from a certain metric threshold (e.g., Pelletier et al., 1997); (2) at the end of the propagation phase on an "eroded" topography, i.e., after the LSgCA has already stabilized all slopes from the pristine modelled topography.Landslides observed afterwards can be explained from weathering, occurring in places that remain unstable (such as scarps); (3) as in (2) but with the water saturation increased to mimic a rainstorm.The propagation phase, being based on a SOC-like process, has in theory the potential to create or alter the landslide FSD.All those steps are illustrated in Fig. 4.
We observe that a power law distribution with a median α = 2.2, as observed in Nature (section 3.1), already emerges at the initiation phase, just based on the FS map, which is in agreement with e.g.Pelletier et al., (1997).Since the FS map directly depends on the topography, we also investigated the role of the fractal dimension on α (see below).Comparison of the FS maps in the left column of Fig. 4 shows that the number of instabilities indeed decreases from pristine topography to "eroded" topography (after application of the LSgCA) and that new instable areas appear if the LSgCA uses the factor of safety where the slope is completely submerged or the water table is present at the surface, instead of the factor of safety in the dry condition where pore water pressures are zero.We also find that whatever the triggering process, such as weathering in the second row or rainfall in the third one, a similar power law distribution is obtained with again a median α = 2.2 (here for 100 simulations).
We did a more systematic analysis based on 1,000 simulations per FS formulation and per step in the landslide process (representing a total of 12,000 simulations).The resulting α distributions are plotted in Fig. 5, and compared with the real distribution.The boxplot indicates the range, the lower quartile, the median, and the upper quartile of the produced α.Once again, we find a good match between the LSgCA and the observations.We also verify that α remains stable over the different steps of the landslide process, demonstrating that the SOC-like behaviour of landslides is not at the origin of α ≈ 2.2.Taking into account all the simulation results shown in Fig. 5, we obtain α = 2.17±0.49.Regarding the inherent variations of α, they are due to the stochasticity of the process, meaning the random variations in the ranges of tested input parameters (Table 2) and in the simulated topography.
Despite the fact that the distribution of α isn't solely dependent on a single parameter, rather to a combination of sets of input parameters and complex geometries, it must be the roughness that plays an important role in the FSD statistics of landslides (as identifiable from the FS maps in Fig. 4).In a second test, we evaluated the role of D f changes (in the range of 2.1 ≤ D f ≤ 2.9) on the estimation of α in the percolation phase (i.e., initiation phase, no propagation).The higher values of D f in the percolation phase imply an increase in the surface roughness.This suggests that we limit the range of our analysis to 2.2 ≤ D f ≤ 2.5 to produce more realistic topographies in order to interpret the behaviour of α (lower and higher values represent either too smooth or too rough topographies -see below).Similar to the analysis in the previous sections, the ground parameters were chosen randomly from the typical range of parameters in Table 2. Once again, the four FS formulations were tested (Table 1).Results are shown in Fig. 6.The boxplot indicates the range, the lower quartile, the median, and the upper quartile of the produced α obtained from a total of 36,000 simulation (1,000 per FS, per D f ).The surface topography of the case D f = 2.1 is too smooth to trigger landslide and produce a meaningful α distribution.The median of the α values observed in the range 2.2 ≤ D f ≤ 2.5 tend to increase for all FS equations.For the case of D f = 2.6, the median value of the α drops down to an approximate value of α = 1.5, due to the high instability of fractal topography and therefore the failure of an individual mega large landslide.A similar result is obtained for D f = 2.7.Above, the fractal topography fails to produce a realistic surface where landslides analysis is able to produce valid results.In the acceptable range of the fractal dimension (2.2 ≤ D f ≤ 2.5), a higher D f means an increase of the ratio of small landslides to the larger ones until they coalesce onto the full topography grid, explaining our observations.

Conclusions
In this study, we presented a landslide generic cellular automaton (LSgCA) that models the propagation phase of landslides.
After revisiting the literature, we updated the observed power-law exponent distribution to α = 2.21±0.53 and validated our LSgCA by obtaining a similar distribution in simulations with α = 2.17±0.49.Although realistic values had already been obtained in previously published models (e.g., Pelletier et al., 1997;Hergarten and Neugebauer, 2000), we here demonstrated that α remains at first-order constant at the different stages of the landslide process.Moreover, despite the α variability being due to different input parameter values and different topography iterations, an increasing trend in the median α was observed as a function of increasing fractal dimension D f , suggesting this parameter as the main parameter influencing α.Indeed, since α ≈ 2.2 is already observed in the initiation phase where the only information is the FS map, α can only find its origin in the underlying fractal topography.As a corollary, it also explains why the landslide trigger type (e.g., weathering, rainfall, earthquake) does not seem to influence α.
We believe that the LSgCA solves the existing problem of abstract physical-based scenario analyses, unrealistic topographies and/or the lack of landslide propagation phase in the landslide CA literature (Piegari et al., 2006;Hergarten and Neugebauer, 1998;Densmore et al., 1998;Guthrie et al., 2008).LSgCA provides a simple tool for further studies in probabilistic multi-hazard assessments by defining landslide footprints that could be included in stochastic landslide sets conditional on rainstorm or earthquake triggers (e.g., Mignan et al., 2014;2018;Liu et al., 2015;Matos et al., 2015).But by proving that the origin of α lies in the topography, landslide hazard can be approximated directly from the FS maps, before any dynamic modelling.

Acknowledgment
This work was supported by the Swiss Competence Center for Supply of Electricity (SCCER SoE) T4.1 "Risk, safety and societal acceptance" and by the Swiss National Science Foundation Program "NRP70" as a part of the project "Risk Governance of Deep Geothermal and Hydro Energy" with the grant number 40 7040_153931.
Table 1.List of Factor-of-safety equations used in the present study (Cruikshank, 2002).
Infinite slope in the dry soil condition Infinite slope with seepage parallel to slope Infinite slope with seepage and tree roots       3).

Figure 2 :
Figure 2: LSgCA application in a virtual region, defined in a 100 km by 100 km grid with smooth topography (as defined in Komendantova et al., (2014), Liu et al., (2015) and Mignan et al., (2017) for generic multi-hazard illustrations).See section 3.2 for the testing of more realistic, fractal, topographies and the supplementary material for the animated version of this landslide.

Figure 4 :
Figure 4: Effects of different steps of the landslide process on the α α distribution.Red cells represent unstable patches in the FS maps.Brown patches represent the thickness of landslides in the landslide maps (scarps are represented in dark grey).FSDs are shown for 100 simulations with the median highlighted in black.

Figure 5 :
Figure 5: Distribution of α α obtained at different steps of the landslide process, for different FS formulations.Each distribution was estimated from 1,000 simulations.Note the stability of α α and the match with observations.