A statistical study of precursor activity in rain-induced landslides

1 Introduction

Although landslide avalanches are well known, the formative mechanism of slip surfaces is usually not well understood. Catastrophic landslide velocities range from 1 m/day to 300 km/h (in case of earthquakes), and their intensity usually increases after prolonged rains or after the frost leaves the ground in the spring. For rain-induced landslides, the most important effect is the reduction in stability resulting from water pressure within the discontinuities and other defects [1–4]. The main failure mechanism during landslides is that shear takes place along either a discrete sliding surface, or within a zone, underneath the surface. If the shear force (driving force) is greater than the shear strength of the interface (resisting force), then the slope will be unstable. Instability could take the form of displacement that may or may not be tolerable, resulting in the collapse of the slope either suddenly or progressively [5]. This is the basic picture envisioned by Zaiser and Aifantis [6, 7], as well as by Zaiser and coworkers [8, 9], in considering the evolution of deformation/slip avalanches when material softening and stochasticity is stabilized by higher-order gradients of the constitutive variables. It is noted, in this connection, that the inclusion of higher-order gradients is necessary for the convergence of numerical solutions of either continuum or discrete models, as first indicated in the pioneer work of Aifantis [10, 11], whereas a more recent account can be found in Ref. [12].

Nowadays, it is well known that landslide phenomena show a fractal character and their occurrence frequencies could be analyzed with respect to the critical affected area through a power law. The relation between self-organized criticality (SOC) and natural hazards (i.e., earthquakes, forest fires, landslides, rockfalls, etc.) has been investigated over the last three decades. Indeed, the concept of SOC applies in avalanche dynamics modeling the tendency of natural systems to approach a critical state where a power-law behavior is followed having an exponent b with its value being robust for a wide range of system parameters. To explore the inherent mechanisms leading to the occurrence of various catastrophic phenomena, several statistical models using SOC features were suggested. One of them is the Olami-Feder-Christensen (OFC) model [13], which was introduced in 1992 as a simplification of the classical Burridge-Knopoff (BK) spring model [14] first introduced in 1967. The OFC model represents a classical non-conservative model revealing SOC.

For landslides, Piegari et al. [15] introduced a cellular automaton (CA) model based on the dynamical evolution of a space- and time-dependent safety factor field. They assumed an anisotropic version of the OFC model in terms of a finite driving rate. In particular, they studied a continuously driven anisotropic CA based on the OFC model, by characterizing landslide size distributions. Their results showed that different distributions for the landslide sizes were attained by different driving rates. The larger the driving rate, the more continuous is the changing of the frequency-size distribution from power law to Gaussian function, whereas the non-conservative redistribution of the load of failing cells is controlled by the anisotropic transfer coefficients.

Zhang et al. [16] introduced a modified version of the OFC model, to model avalanche size differences, by using a concept of weighted edge in order to improve the original redistribution rule and to obtain an inhomogeneous network with different local frictions and elastic behavior. Recently, Avlonitis and Papadopoulos [17] proposed, within the aforementioned gradient approach for softening materials, a generalized OFC model incorporating softening behavior within the interface of two moving plates. They have shown that as the interface enters into the softening regime, a significant change of the exponent b occurs, as by the governing gradient-dependent evolution equation.

In this work, we combine all the aforementioned findings to show that a gradual change in the power-law exponent b of an unstable slope may be used as an indication for possible rain-induced landslides.

2 Shear strength softening law

As essentially suggested by Zaiser and Aifantis [6, 7], as well as by Zaiser and coworkers [8, 9], and recently pursued in the present context by Avlonitis and Papadopoulos [17], a system of two distinct plates moving relatively to each other is assumed, with their interface being represented with a number of blocks (nodes) interconnected with springs of a certain elastic constant (Figure 1).

Figure 1

The generalized OFC model. Two distinct plates moving relatively to each other. Their interface is represented with a number of blocks (nodes) interconnected with springs of a certain spring constant.

A dynamic variable f_ij is assigned to each node of a two-dimensional grid (i, j denoting the number of rows and columns of the grid, respectively) that represents the interface of the plates. Assuming isotropy at the horizontal plane, f_ij=Kℓ²∇²u_ij+K_Lδu_ij, where u_ij is the node’s displacement at x_ij; K, K_L are the spring (elastic) constants; and δu_ij=υt-u_ij, where υ is the constant velocity of the upper moving plate and υt is the imposed driving displacement. The second-order gradient term in the right-hand side models (as usually in the continuum limit) the expression 2u_ij-u_i-1j-u_i+1j+2u_ij-u_ij-1-u_ij+1, and ℓ (is assumed unity) is a characteristic scale associated with space discretization. The dynamic variable f_ij is the total force (or stress divided by the unit area).

When the force on an arbitrary node is larger than a threshold value F_th, which is the maximal static friction, the node slips. This implies that for the node n_ij of the fault plane to become unstable, the following inequality must be fulfilled:

In Ref. [17], it was assumed that the threshold value F_th of the maximal static friction may be interpreted as a dynamic internal variable that varies as the process evolves, leading to a drastic change in the dynamic behavior of the spring-block model.

To examine the origin of the slope’s weakening that emerged in the shear strength law, let us assume a constitutive stress-displacement law where both a hardening and a softening regime are involved (Figure 2). We further assume that as the density of the water within the slope’s interface increase (due to rain), the induced increase of the pore hydrostatic pressure at each site of the interface results to a corresponding modification of the flow stress (and of the stress-displacement law in general), the main effect being the decrease of the shear strength at this site.

Figure 2

Constitutive stress-displacement law, while a hardening and a softening regime are present.

It is assumed that as the evolution of the rainfall proceeds, the local shear strength evolves from τ₁, to τ₂, and then to τ₃. Note that there is a reduction of the slope, both in the hardening and softening regimes (Figure 2). In this case, as the shear strength of the interface evolves from τ₁, to τ₂, and then to τ₃, the hardening slopes’ E_i’s evolve from E₁, to E₂, and then to E₃.

In Figure 3, the evolution of the shear strength, for the general case, is depicted as a function of displacement. As in Refs. [6–9], two distinct constitutive regimes are assumed: the weakening and strengthening regimes where the shear strength decreases and increases its value, correspondingly. For the specific case of rain-induced landslides, the weakening regime can be understood in the following scheme: as the deformation increases, the pore deformation increases resulting to high water concentration in these sites. In turn, high water concentration induces a gradually increasing shear strength reduction due to the pore hydrostatic pressure. As result, the higher the deformation, the higher the shear strength reduction.

Figure 3

The shear strength as a function of the displacement.

It is noted that, in the more general case, a corresponding strengthening regime for the shear strength law may also be assumed [6]. Indeed, for interfaces in which after corruption, restrengthening between the sliding plates emerged, e.g., as is the case for earthquakes faults, it is possible to assume the inverse evolution as demonstrated in Figure 2. In this case, the shear strength evolves from τ₃, to τ₂, and then to τ₁.

Rather general and for the simplest case, a linear idealization of the strengthening-weakening shear strength curve is adopted,

where λ is a proportionality parameter modeling the strengthening or weakening process of the geomaterial (i.e., intact rock, rock mass, soil) and F^* a constant parameter. The sign of the proportionality parameter depends on which branch of the shear strength law (strengthening-weakening) of the fault plane the material deforms.

This is exactly the case examined in Ref. [17] where a non-local behavior of the displacement manifold is assumed, i.e., uij=uij0+c∇2uij (uij0 being a uniform over the entire space displacement field and c the gradient coefficient) and substituting in Eq. (2),

with Fth0=F*+λuij0 and K′=λc. As a result, relaxing the assumption of homogeneous (constant) shear strength and interpreting it as an internal spatially changing variable, the final inequality [Eq. (1)] for an arbitrary node n_ij of the interface to become unstable reads

where K′=λc, λ being a proportionality parameter modeling the strengthening or weakening process of the shear strength law (which does not coincide with the classical softening-hardening regime in the corresponding stress-strain law) and c is the gradient coefficient modeling the non-uniform character of the interface. Fth0 is an average local shear strength.

Equation (4) is of crucial importance as it redresses the generic weakness of the OFC model, while it models not only the elastic properties of the interface but also the structural heterogeneity through the gradient coefficient c. More important, it models the deformation of the interface under a softening or hardening mode expressed by the appropriate sign – minus or plus – of the proportionality parameter λ of the shear strength curve. The approach and theoretical framework just described above, as well as its implementation to a CA as discussed below, is directly motivated by similar works conducted for analogous avalanche-type deformation/slip processes [6–9, 12, 18–21].

3 CA simulations

As already mentioned, the aforementioned OFC modification resembles to a great extent (in fact, is directly motivated from) the stochasticity-enhanced gradient plasticity model first introduced by Zaiser and Aifantis [7] for modeling slip avalanches, later used for snow avalanches [8], thin film delamination [9], and compression of metallic foams [19] and micropillars [20]. The CA procedure used in the above-mentioned references is also used here, as briefly discussed below.

The system is loaded by considering the external stress (force) to increase with small steps. For each site, the inequality given in Eq. (1) is examined and, if it holds, then the site slips. When a quasi-static equilibrium is reached, the number of slipped or failing sites is measured and the external force is then increased again, with the same procedure repeated until an avalanche occurs and the whole system slips. The number s of slipped sites, i.e., the size of avalanche occurring, is then divided by the total number of avalanches measured to produce the probability of having an avalanche of a certain size (number of sites failing simultaneously).

In the following, we depart from Zaiser et al. [8] by relaxing the assumption that there is a constant rate between external stress and local displacement increments. As seen in Figure 2, this is attributed to the departure from linearity in the local stress-displacement constitutive relation of each site when water is present in the interface, as explained in Section 2.

This has the following consequence: the relative magnitude of stress and strain increment in the corresponding simulation code is not constant or, as it was assumed by Piegari et al. [15], different values of driving rates are present in relation to the different slopes E_i. The driving rates are analytically computed in the Appendix for all possible different scenarios of deformations, i.e., for hardening-softening stress-displacement law, as well as strengthening-weakening of the shear strength law. For the present study, we use the driving rates corresponding to a weakening behavior of the slope deforming under the softening regime, the rates being

From Eq. (5), the driving rate d(s)=dτ/du=-A is estimated, where in a non-dimensional space and for the range of b values encountered in real fields, its value ranges from 10^-3 to 10^-2.

Combining results from 100,000 such simulations provides the power-laws shown in Figure 4A–C for different values of the slope λ. It can be seen that for different slopes λ, i.e., for different relative magnitudes of stress and displacement, a different exponent b of the corresponding power-law distribution emerged, and more specifically as the shear strength decreases, i.e., the slope λ decreases, the absolute value of the exponent b also decreases. The larger the reduction of the shear strength, the lower the exponent b of the emerging power law (Figure 4D).

Figure 4

(A–C) The distribution of the slip events vs. their size s. The dependency of the b value from the weakening rate λ, for driving rate d(s)=0.007 and structural heterogeneity w=5. (D) The larger the reduction of the shear strength, i.e., the smaller the weakening rate λ, the lower the exponent b of the emerging power law.

In Figure 5D, the dependence of the b value on the structural heterogeneity of slope’s interface is examined. Here, different structural heterogeneities are modeled by introducing different distributions of the deformation values of the initial sites. Indeed, the different curves correspond to different Weibull distributions of the initial site deformations, u_ij. As in Refs. [22, 23], different simulations were performed for different shape parameters w of the Weibull distribution {p(x; m, w)=(w/m)(x/m)w−1exp[-(x/m)w], x≥0} of u_ij, thus imposing gradually stronger heterogeneity: a higher w corresponds to a higher variability in the corresponding distribution and, as a result, to a higher heterogeneity. It can be seen that the stronger the heterogeneity is, the higher the b value is, in accordance with our derivation in Eq. (6) and because, as we have noted before, the dissipation parameter α and the b value always evolve reversely.

Figure 5

(A–C) The distribution of the slip events vs. their size s. The dependency of the b value from the structural heterogeneity of slope’s interface, for driving rate d(s)=0.007 and weakening rate λ=0.019. (D) The larger the w (structural heterogeneity), the lower the exponent b of the emerging power law.

4 Discussion

The simulation results depicted in Section 3 are consistent with those of Ref. [17], as both approaches are based on the same theoretical foundation. To see this more clearly, we revisit the redistribution parameter

as it was defined in the original OFC simulation code. Here, K″ is the coefficient of the Laplacian of displacement. While in the initial OFC model K″ coincides with K, the elastic constant of the interface at the horizontal direction, here we have K″=K·ℓ²-K′, where K′=λc as it was defined in Eq. (3) with ℓ=1. For the problem at hand, as weakening of the slope evolves owing to rain-induced load, the corresponding weakening rate λ^s decreases (see Appendix). As a result, K′ decreases, K″ increases, and, finally, from Eq. (6), α also increases. Noting that α always evolves inversely proportionally with the exponent b, the expected outcome from the formalism of Ref. [17] is that as the slope evolves, the value of exponent b of the corresponding distribution of microslip events within the interface will gradually diminish. Indeed, this is in accordance with the simulation results depicted in Figure 4D.

Within the proposed framework presented in the previous sections, a robust study of catastrophic rain-induced landslides can be achieved both theoretically and through simulations. It was shown that not only the effect of water concentration in the slopes’ interface can be modeled, but also the role of the mechanical properties of the slope, such as heterogeneity of the shear strength within the interface, can be accounted for. The basic phenomenology of the slope’s behavior is sufficiently considered by the approach adopted in Ref. [17], as well as the simulation algorithm presented here through the notion of altering driving rates.

Indeed, and on another example, the effect of the slopes’ tilt can also be studied by means of the notion of evolving driving rates. As the tilt of the slope increases, the weight of the overburden on intact rock mass or soil (at each instant of time) results in a higher incremental load per unit time to the interface. It is expected that slopes with a higher tilt will be characterized by a higher probability for large slip events. In Figure 6, the simulation outcomes for three different tilts are shown, where the three different tilts in the simulation code were identified with three different rates of stress increase (per simulation run) altered by one order of magnitude. It can be seen that for higher rates of incremental stress evolution, i.e., for higher slope tilts, the absolute b value is lower, which, in turn, means that large slip events are more probable than small slip events at the interface. As a result, the higher the tilt is, the higher the probability for large slip events is, in accordance with real field observations.

Figure 6

(A) The distribution of the slip events vs. their size s. (B) The dependency of the b value from the driving rate d(s), for the heterogeneity parameter w=15 and weakening rate λ=0.005. The larger the driving rate d(s), the lower the exponent b of the emerging power law.

As a conclusion, it is argued that the above analytically predicted and, through simulation, evaluated behavior of the slopes’ interface for rain-induced landslides can be used to build a novel tool for predicting (the onset of) catastrophic phenomena. More specifically, it is proposed that the decrease of the b exponent of the slopes’ power-law distribution can act as a quantifiable precursor indicator for catastrophic rain-induced landslides.