Recurrence time distribution and temporal clustering properties of a cellular automaton modelling landslide events

Reasonable prediction of landslide occurrences in a given area requires the choice of an appropriate probability distribution of recurrence time intervals. Although landslides are widespread and frequent in many parts of the world, complete databases of landslide occurrences over large periods are missing and often such natural disasters are treated as processes uncorrelated in time and, therefore, Poisson distributed. In this paper, we examine the recurrence time statistics of landslide events simulated by a cellular automaton model that reproduces well the actual frequency-size statistics of landslide catalogues. The complex time series are analysed by varying both the threshold above which the time between events is recorded and the values of the key model parameters. The synthetic recurrence time probability distribution is shown to be strongly dependent on the rate at which instability is approached, providing a smooth crossover from a power-law regime to a Weibull regime. Moreover, a Fano factor analysis shows a clear indication of different degrees of correlation in landslide time series. Such a finding supports, at least in part, a recent analysis performed for the first time of an historical landslide time series over a time window of fifty years.


Introduction
Analysis of recurrence times in natural series is a challenging issue in many branches of science, as discovering temporal correlations between events is crucial to determining the probability of occurrence of future events.After the seminal paper of Bak, Tang and Wiesenfeld (BTW), who introduced the concept of self-organized criticality (SOC) to explain the scale invariance of many natural phenomena (Bak et al., 1987), several authors have searched for correlation functions of power-law type in time series.Although the original BTW model provides uncorrelated exponentially (Poisson) distributed recurrence times (Boffetta et al., 1999), several variants of the BTW model that supply power-law behaviour of recurrence time probability distributions have been proposed (Sanchez et al., 2002;Paczuski et al., 2005).Recurrence time statistics of solar flares is found to be characterized by power laws (Boffetta et al., 1999) as well as statistics of solar flare magnitudes (Baiesi et al., 2006).Powerlaw temporal correlations have also been found for recurrence time probability distribution of earthquakes above a given magnitude (Corral, 2004).On the other hand, several authors have demonstrated the validity of stretched exponential distributions to characterize return interval statistics between extreme events in the case of long-term correlated records (Altmann et al., 2005;Bunde et al., 2005;Eichner et al., 2007).Furthermore, recent works have shown that the Weibull distribution, which is the product of a power-law and a stretched exponential, best represents the recurrence time distribution of observed time series of tropical temperature and humidity (Blender et al., 2008) and historical landslides (Witt et al., 2010).The Weibull distribution is also found by Santhanam and Kantz (2008), who study long-range correlated time series and investigate the role played by the threshold in modifying the return interval distribution, especially in the power-law regime related to short return intervals.
In this paper, we analyse the recurrence time statistics of a cellular automaton (CA) model aimed at reproducing the main features of landslide event statistics by means of some characteristic parameters (Piegari et al., 2006a(Piegari et al., , 2009a)).
We note that long-term correlated time series are not assumed, but result from different CA simulations of landslide events.Interestingly, we find that the recurrence time probability distribution of the CA strongly depends on the rate ν at which the system is driven to instability.This rate controls how quickly the CA cells reach the instability threshold and, therefore, its values describe response time scales of slope systems rather than triggers of different origin.In this respect, we showed in a previous study (Piegari et al., 2009a) that similar values of ν can explain the similarity of frequency-size distributions of landslides induced by different triggering mechanisms (Malamud et al., 2004).By increasing the values of ν, a smooth transition from a powerlaw to a non power-law regime was found by analysing the behaviour of landslide frequency-size distributions (Piegari et al., 2009a).Here, we show that by increasing ν, the recurrence time probability distribution smoothly changes its shape as the effect of a smooth transition from a power-law regime to a Weibull regime.Finally, temporal correlations of synthetic time series of landslide events are examined by calculating the Fano factor, and critical variations of model parameters for recurrence time statistics are discussed in relation to previous studies on frequency-size probability distributions.

The model
The CA model is defined on a two-dimensional grid of L×L square cells.The state of each cell i is characterized by the value of a single dynamical variable, e i , representing the ratio of the load acting on the cell to the maximum allowable stress, i.e. the inverse of the so-called factor of safety (FS) e i = 1/FS i .If e i < 1, the acting load is within the allowable limit and the corresponding cell i is stable.If e i = 1, the acting load is at the allowable limit and, thus, values of e i ≥ 1 are associated with failure.We start from a random initial stable configuration, i.e. we attribute to each cell of the grid a uniformly distributed random value of e i with 0 < e i < 1∀i.The dynamics of the CA model is defined by the two typical rules of the SOC theory (Bak et al., 1987).The first rule is an overall driving that provides an increase in e i at the same rate ν approaching the system to the instability threshold, e th = 1.The model parameter ν controls the rate at which all sites are driven towards instability and, in the following, simulation results will be discussed by varying such a key parameter.The second rule is a relaxation rule defined by the equations: where nn(i) denotes the four neighbour sites of the cell i (i.e.nn(i) = up, down, left, right).To aid the reader, a list of variables used in the text is given in Table 1.Eq. ( 2) N(T i ) Number of events above threshold q and over time series with length W .

Sect. 4 q
Threshold value above which s ≥ q is considered an extreme event.
Sect According to Eq. ( 1), when a cell becomes unstable (i.e. e i ≥ e th ), it affects, via a chain process, the stability of the neighbour cells, as a fraction f nn(i) of e i toppling on nn(i).After a failure, we set e i = e min with e min = 10 −6 ; however, any other finite level would work (Jensen, 1998).In general, the transfer coefficients f nn(i) , which describe the distribution of over-threshold cell loads to the nearest neighbours, can be anisotropic.We set f left = f right and f down > f left , f right , f up to mimic the preferential direction for stress redistribution induced by gravity.We note that if anisotropic stress redistribution is related to specific slopes, different choices for the transfer coefficients, f nn , can be made to take into account the topography of individual investigated areas (Piegari et al., 2009b;Di Maio and Piegari, 2012).
During each iteration of Eq. ( 1), an amount of e i is lost from the system, i.e. the difference between e i and the sum of the amount f nn e i added to each of the four neighbour sites is not equal to zero.In other words, we take into account non-conservative (dissipative) cases where the quantity C = nn f nn , which fixes the degree of conservation of the system, is less than 1, contrary to other approaches (Hergarten and Neugebauer, 2000;Hergarten, 2003).The role of the model parameter C has been investigated in previous works, where we found that the shape of the frequencysize probability distribution of landslide events is strongly affected by the values of C (Piegari et al., 2006a), while it is not significantly affected by the values of the coefficients f nn in the range of values that supply power-law distributions (Piegari et al., 2006b).In particular, only if C < 1 do our synthetic frequency-size distributions reproduce those from catalogues (Piegari et al., 2009a), pointing out the relevance of dissipative phenomena to landslide triggering.Due to the peculiar role of C, in the next sections we will discuss simulation results by varying C, while we fix the ratios between the transfer coefficients equal to: f up /f down = 2/3 and f left,right /f down = 5/6 for all simulations.
Finally, we note that the inter-event occurrence time statistics of the CA is studied once the system has attained a stationary state in its dynamics, i.e. the mean value of the dynamical variable, e i , on the grid sites fluctuates about an average value.
In the following, simulation results are shown for the case of a square lattice of linear size L = 64.

Recurrence time statistics
The recurrence time, τ , is the time between two events whose size, s, is above a given threshold q (Fig. 1a).In the literature, τ is also named inter-extreme or inter-event occurrence time, as its definition requires the choice of a threshold.In the following, we analyse the recurrence time distribution of the CA for different thresholds q and by varying the key model parameters C and ν.In Fig. 1b, a simplified flow chart of the CA showing the algorithm for the calculation of the recurrence time distribution, p(τ ), is reported.Simulation results will be also discussed in light of the results obtained by Witt et al. (2010), who analyse for the first time temporal correlations and clustering in an extensive landslide data set.These authors discuss statistics of two landslide intensity time series obtained by introducing two kinds of threshold: the number of reported landslides in a day, D L , and the number of reported landslides in an event, S event , where an event includes mostly consecutive days with landsliding associated with the same triggering episode (rainfall in their study).In our CA model, the threshold value q is defined by s, i.e. the total number of cells that reach instability in a landslide event.We  note that all relaxation processes, caused by the reaching of instability for one or a group of cells (see Eq. 1) and induced by any trigger, last until all e i < e th and are considered to be instantaneous compared to the time scale of the overall drive ν.Consequently, the model is not able to distinguish the two measures of landslide intensity used by Witt et al. (2010), as they are the same in our case.Finally, it is worth noting that ν being not exclusively related to rainfall trigger, our numerical analysis may be helpful for interpretation of recurrence time statistics of landslide events related to different types of trigger mechanisms.

Analysis of varying threshold values
To investigate how the choice of q affects the recurrence time distribution p(τ ), we analyse the behaviour of p(τ ) considering different thresholds and keeping the values of ν and C fixed.In Fig. 2a and b, we show the behaviour of p(τ ) for q = 3, 5, 30 and 50, obtained by using two values of ν corresponding to different regimes of the model (see next section).As results from Fig. 2a, for a small value of ν, p(τ ) exhibits a linear regime in a log-log scale -that means a power-law regime in a linear scale -for recurrence time intervals that enlarge by increasing the threshold q.We note that in the limit of small values of ν, the model resembles the results of the OFC model (Olami et al., 1992), as extensively discussed in previous works (Piegari et al., 2006a(Piegari et al., , 2009a)).Thus, in this 2. Recurrence time distribution p() for different values of the threshold q, obtained by the system towards instability with rates (a)  = 10 -4 and (b)  = 0.003.In the simulations, ue of the conservation parameter C is kept constant, C = 0.4.regime p(τ ) is expected to be characterized by a power-law regime according to SOC-like models proposed for explaining recurrence time statistics of earthquakes and solar flares (Sanchez et al., 2002;Corral, 2004;Pacuzski et al., 2005).Interestingly, the power-law regime is better defined for the occurrence time statistics of the most extreme events (q = 30, 50).It is also worth noticing that variations of q in Fig. 2a affect the shape of p(τ ) either for short or for large return intervals.Conversely, when a larger value of the driving rate is considered (Fig. 2b), the behaviour of p(τ ) is affected by the threshold q mainly for short return intervals, as is found in Santhanam et al. (2008).Moreover, a power-law regime, which is still recognizable for short τ and large thresholds q, smoothly disappears with decreasing values of q.

Analysis of varying driving rates
In Fig. 3, we show the behaviour of p(τ ) by varying the driving rate ν and keeping the threshold q and the conservation parameter C fixed.As discussed in previous papers (Piegari et al., 2006a(Piegari et al., , 2009a)), the parameter ν, which describes the rate of approaching instability, is a key ingredient in the model, as its values are responsible for a change in the dynamics of landslide processes.For small ν, a few cells (a single cell in the limit of vanishing ν) initially reach the instability threshold and chain processes dominate landslide dynamics.By increasing ν values, an increasing number of cells initially reaches instability, so domino processes become less and less effective in causing avalanching processes, and instability is reached simply because a very large number of cells almost instantaneously reach the critical threshold.Such a change in the model dynamics controlled by ν corresponds to a smooth crossover in the frequency-size distribution from a power-law to a Gaussian-like regime (Piegari et al., 2006b).It is worth noting that only if ν > 0 does the automaton provide frequency-size distributions characterized by an exponential rollover for small landslides followed by a power-law regime for medium and large landslides, as observed for landslide frequency-size distributions from catalogues (Piegari et al., 2009a).In particular, since landslide inventories corresponding to different triggering mechanisms seem to obey the same frequency-size probability distribution (Malamud et al., 2004), we realized that triggers of different origin can affect the stability of a slope with the same rate of approaching instability, i.e. with the same ν and, therefore, the probability of occurrence of landslide events is controlled by the rate ν, rather than the nature of the triggering mechanism (Piegari et al., 2009a).
Interestingly enough, as shown in Fig. 3, the values of ν are also responsible for a change in the behaviour of the probability distribution of recurrence times.For small ν, p(τ ) is characterized by a power-law regime typical of SOC-like systems (Sanchez et al., 2002;Corral et al., 2004;Paczuski et al., 2005) in a range of return intervals that becomes shorter and shorter by increasing the values of ν. Surprisingly, if we use values of the model parameters for which a good matching between synthetic and real frequency-size probability distribution is obtained, i.e. ν = 0.003, C = 0.4 (Piegari et al., 2009a) and q = 5, we find that the behaviour of p(τ ) is characterized for small τ by a power law: with a negative exponent α = −2.67,while for medium and large τ it is described by a stretched exponential decay well approximated by a two-parameter Weibull distribution: with 1/β = 1.15 and the shape parameter γ = 0.82.Interestingly, the value of the power-law exponent α is in very good agreement with the experimental value found by Juanico et al. (2008), while the value of the Weibull shape parameter γ is in the range of values 0.67 < γ < 0.83 found by Witt et al. (2010) to characterize the landslide intensity of time series for S event (see Sect. 3).
In the case γ = 1, the Weibull distribution (3) reduces to an exponential (Poisson) distribution characteristic of uncorrelated time series.Thus, we deduce that in the investigated cases the simulated landslide events are correlated in time being γ < 1.Finally, it is worth noting that in the region of model parameters for which a good agreement with experimental data is obtained, simulated recurrence time distributions cannot be entirely approximated by Weibull distributions (see Fig. 3, square symbols), as the region of short τ requires the introduction of a different power-law exponent.

Analysis of varying conservation degrees
In previous works (Piegari et al., 2006a, b), we have shown that the model parameter C, which fixes the degree of conservation in the system, strongly affects the shape of the frequency-size probability distribution of landslide events.In particular, we recovered the behaviour of the frequency-size distributions from catalogues (Malamud et al., 2004) only for values of C ≈ 0.4-0.5, indicating that the inclusion of dissipative phenomena is obligatory for describing realistic scenarios.Thus, we investigate the role of C in determining the recurrence time distribution.In Fig. 4, we show the behaviour of p(τ ) for different values of C and keep the value of q and ν fixed.In particular, we set q = 5 and choose a   ν value equal to 0.003, for which the best agreement with data from catalogues was obtained for the frequency-size probability distribution (Piegari et al., 2009a).As shown in Fig. 4, variations of C affect the behaviour of p(τ ) mainly for large return intervals.Specifically, the probability of occurrence of large recurrence times increases by decreasing the degree of conservation in the system.As discussed in the previous section, the behaviour of p(τ ) for medium and large τ is well described by a Weibull distribution with values of the shape parameter γ less than 1.To analyse how the shape of the stretched exponential, which best approximates p(τ ), evolves with C, we also plot in Fig. 4 two Weibull distributions with 1/β = 1.15 and shape parameter values γ = 0.82 and γ = 0.94.We note that the larger the degree of conservation, the larger the value of γ becomes.Such a feature can be explained considering that in the SOC limit, realized for C = 1 and ν = 0, the model is expected to provide uncorrelated exponentially (γ = 1) Poisson distributed recurrence times (Boffetta et al., 1999).

The Fano factor
To reveal fractal structures in temporal series, it is possible to analyse variations of the values of the Fano factor, which is defined as the ratio of the variance to the mean for windowed data (Fano, 1947).Let W be the time window and n the number of non-overlapping adjacent intervals of length T = W/n.The Fano factor is defined as:  4b and 4c) and an uncorrelated regime for  = 6 0.01 (Fig. 5d).The oscillatory behavior of F(T) in (d) is a purely numerical effect due to the finite 7 system size: it is shifted to larger time interval lengths T if larger grid size are used.where σ 2 [N (T )] and N (T ) are, respectively, the variance and the mean of the number of timings N (T i ) per ith interval, with i = 1, . . ., n.For an unclustered set of timings N (T i ), the inter-event occurrence times are Poisson distributed (Ross, 1983), and consequently the variance is equal to the mean for all interval lengths T : If such a condition is substituted in Eq. ( 4), one gets F (T ) = 1 for uncorrelated time series.In a different way, for fractal clusters of timings, i.e. timings unequally spaced with a fractal (scale invariant) structure, it has been shown that σ 2 [N (T )] = N (T ) ψ+1 (Teich et al., 1997).Then, if this last condition is substituted in Eq. ( 4), it turns out that the Fano factor depends on the time interval length T with a power-law exponent ψ: with constant A = N (W ) /W ψ .The case of Poisson distributed recurrence times, i.e.F (T ) = 1, is found if ψ = 0 for all T .
We have analysed the variation of F (T ) for simulated landslide event series with four thresholds q, and for different values of the driving rate ν.
Looking at the panels of Fig. 5, which show the evolution of F (T ) with increasing driving rate, it is apparent that F (T ) = 1 in all cases with ν < 0.01 (Fig. 5a, b and c).This result provides additional information on the system dynamics; specifically, by varying ν the system undergoes a transition between two regimes: a correlated regime (ν < 0.01), in which events are clustered in time with recurrence times governed by power-law and Weibull distributions (see previous section), and frequency-size distributions characterized by power-law decays (Piegari et al., 2009a), and an uncorrelated regime (ν > 0.01), in which events are Poisson distributed in time, and frequency-size distributions are characterized by bell-shaped curves (Piegari et al., 2009a).For the correlated regime, the Fano factor shows a power-law dependence on T , with values of the power-law exponent ψ that decrease with the threshold q and vary in the range 0.3 < ψ < 0.5.Interestingly, this range of ψ values is the same reported by Witt et al. (2010) for the Fano factor analysis carried out for two kinds of landslides series recorded in the Emilia-Romagna region (northern Italy).It is also worth noting that for values of ν ≈ 10 −3 , for which our synthetic distributions well match frequency-size distributions from catalogues (Piegari et al., 2009a), the behaviour of F (T ) shown in Fig. 5b reproduces well the behaviour of F (T ) derived from historical data (Witt et al., 2010), as well as the flattening of the F (T ) curve by increasing the threshold value.For larger values of ν (see Fig. 5c) and for interval lengths > 200, a clear change in the slope of the F (T ) function is found, with values of the power-law exponent ψ greater than 0.9.Finally, for ν > 0.01, the case of Poisson distributed recurrence times is found, being F (T ) = 1 for all T .In this regard, we note that the oscillatory behaviour of F (T ) is a purely numerical effect due to the finite size of the grid, and, therefore, it will be shifted to a larger T if the size of the system is increased.

Discussion and conclusions
We have analysed the recurrence time statistics of a cellular automaton modelling landslide events by performing a numerical analysis in the parameter space and estimating Fano factor behaviours.The model is an extended version of the OFC model, which itself is a kind of paradigm for SOC in non-conserved systems (Hergarten and Krenn, 2011), but it works differently from the original OFC model as a finite value of the driving rate ν is applied.By increasing the values of ν, which quantifies the rate of approaching instability, the system undergoes a smooth transition from a correlated regime characterized by power laws to an uncorrelated regime in frequency-size statistics (Piegari et al., 2009a).We found that such a transition is an indicator of the change in predominant mechanisms to propagate instability: for small values of ν, chain processes dominate the landslide dynamics, while for large values of ν such processes are no longer active, and large events simply occur because a large number of cells simultaneously reach instability.By comparing synthetic frequency-size probability distributions with those from landslide catalogues, quite good agreement was found for ranges of driving rate values that characterize the crossover region between the two different regimes (Piegari et al., 2009a).In this paper, we have shown that synthetic recurrence time distributions and Fano factor behaviours resemble those from historical time series (Witt et al., 2010) in a range of ν values that again characterize a crossover region between a correlated (power-law) regime and an uncorrelated (Weibull) regime.The variation of the driving rate in a uniform way might raise doubts about the applicability of our analysis to real-world cases.According to a previous study (Piegari et al., 2009a), the model grid can simulate subsectors of the region where the Witt et al. (2010) data were collected, over which the triggering mechanism can be considered, in a reasonable first approximation, homogenously active.A constant driving rate is representative of an averaged rapidity with which instability is reached somewhere in the system.Actually, approximate constant driving rates for any couple of the following events can be calculated and, in principle, they are different from each other.However, in the real world, luckily, the number of occurrences of strong perturbations, described by large values of the driving rate, is definitely much smaller than that of small perturbations, and these latter ones will dominate the inventory.Accordingly, we reasonably assume that the small approximate constant driving rates dominate the statistics obtained from the inventory.Thus, if small approximate constant driving rates are (Gaussian) statistically distributed, we use, in a first approximation, its mean value to characterize the statistics of the whole inventory.The soundest justification for our assumptions is provided by the good agreement found between synthetic and real data for a given range of approaching instability rates.This suggests the existence of average characteristic response times for loading critical conditions in slope systems, which statistically characterize a crossover region between a correlated and an uncorrelated regime.
Summarizing, if statistical analyses are performed on landslide inventories, a qualitative comparison between real and synthetic data hints at averaged values of model parameters that statistically characterize the data set.From the comparison with the extensive landslide data set used by Witt et al. (2010), our numerical analysis suggests that statistics of such landslide data can be described by a crossover region between a correlated regime and an uncorrelated regime, where recurrence time distributions are characterized by power-law and Weibull behaviours for short and long return times, respectively.Finally, in such a region of the parameter space, clear indications of temporal correlations and clustering by the Fano factor behaviours support, at least in part, the analysis performed by Witt et al. (2010).

1 2
Figure 1.(a) Graph showing the definition of the recurrence time ; (b) simplified flow cha Fig. 1.(a) Graph showing the definition of the recurrence time τ ; (b) simplified flow chart of our CA model describing the algorithm for calculation of the recurrence time distribution p(τ ).

Fig. 2 .
Fig. 2. Recurrence time distribution p(τ ) for different values of the threshold q, obtained by driving the system towards instability with rates (a) ν = 10 −4 and (b) ν = 0.003.In the simulations, the value of the conservation parameter C is kept constant, C = 0.4.

1Figure 3 .
Figure 3. Recurrence time distribution p() for different values of the driving rate .In the 2

Fig. 3 .
Fig. 3. Recurrence time distribution p(τ ) for different values of the driving rate ν.In the simulations, the values of the conservation parameter C and the threshold q are kept constant (C = 0.4, q = 5).Also shown are the power-law distribution with negative exponent α = −2.67 (short dash line) and the Weibull distribution with 1/β = 1.15 and shape parameter γ = 0.82 (solid line). 1

Figure 4 .
Figure 4. Recurrence time distribution p() for different values of the conservation value C 2

Fig. 4 .
Fig. 4. Recurrence time distribution p(τ ) for different values of the conservation value C.In the simulations, the values of the threshold q and driving rate ν are kept constant (q = 5 and ν = 0.003).Also shown are two Weibull distributions with 1/β = 1.15 and shape parameters γ = 0.82 (black line) and γ = 0.94 (magenta line).

Figure 5 .
Figure5.Fano Factor analysis for simulated landslide event series with four thresholds q, and for

Fig. 5 .
Fig. 5. Fano factor analysis for simulated landslide event series with four thresholds q, and for different values of ν.The behaviour of F (T ) shows fractal organization of landslide events, with power-law exponents ψ > 0 for ν < 0.01 (a, b and c) and an uncorrelated regime for ν = 0.01 (d).The oscillatory behaviour of F (T ) in (d) is a purely numerical effect due to the finite system size: it is shifted to larger time interval lengths T if larger grid sizes are used.

Table 1 .
Variables used in the text.