The Impact of Probability Density Functions Assessment on Model Performance for Slope Stability Analysis

Abstract

The development of forecasting models for the evaluation of potential slope instability after rainfall events represents an important issue for the scientific community. This topic has received considerable impetus due to the climate change effect on territories, as several studies demonstrate that an increase in global warming can significantly influence the landslide activity and stability conditions of natural and artificial slopes. A consolidated approach in evaluating rainfall-induced landslide hazard is based on the integration of rainfall forecasts and physically based (PB) predictive models through deterministic laws. However, considering the complex nature of the processes and the high variability of the random quantities involved, probabilistic approaches are recommended in order to obtain reliable predictions. A crucial aspect of the stochastic approach is represented by the definition of appropriate probability density functions (pdfs) to model the uncertainty of the input variables as this may have an important effect on the evaluation of the probability of failure (PoF). The role of the pdf definition on reliability analysis is discussed through a comparison of PoF maps generated using Monte Carlo (MC) simulations performed over a study area located in the Umbria region of central Italy. The study revealed that the use of uniform pdfs for the random input variables, often considered when a detailed geotechnical characterization for the soil is not available, could be inappropriate.

1. Introduction

In many areas of the world, rainfall-induced landslides represent a relevant threat to the population, infrastructure, buildings, and cultural heritage. In recent years, extreme rainfall events have induced an increasing frequency of slope movements [1,2,3,4,5]; therefore, the prediction of rainfall-induced landslides represents a major challenge for the scientific community.

The most damaging landslides are triggered by intense or prolonged rainfall [6,7,8,9], and the most common phenomena are shallow landslides [10].

Considering the negative impact of landslides on society [11,12], different approaches have been developed to protect and safeguard the territory; for example, geomorphological mapping, analysis of landslide inventories, heuristic terrain, and statistically based classification methods [13,14] have often been used to evaluate the landslide susceptibility of an area.

Landslide susceptibility (S) is the likelihood of a landslide occurring in an area on the basis of local terrain conditions [15]. It represents an estimate of “where” landslides are likely to occur.

Susceptibility (S), hazard (H) (representing the probability that a landslide of a given magnitude will occur in a given period and in a given area), and vulnerability (V) (considering people and infrastructures involved) define landslide risk, which is therefore based on what has happened in the past.

Current extreme regimes are not comparable with past rainfall events; therefore, the zoning of risk (function of S, H, and V) is not always sufficient to guarantee the protection of the territory and people.

A robust and sustainable response to planning policies is represented by landslide forecasting models that are able to simulate slope stability as a function of expected meteorological conditions and physical characteristics of the territory.

Spatial and temporal forecasting of shallow landslides triggered by rainfall can be performed according to different approaches: empirical methods, which analyze records of landslide events and attempts to determine spatial and temporal variations in the occurrence and frequency of landslides [16]; and physically based (PB) approaches [17,18,19,20] accounting for the local physical and mechanical properties that control the failure processes [21,22,23,24,25,26,27,28]. The latter approaches are preferred to forecast the spatial and temporal occurrence of shallow landslides triggered by individual rainfall events, and they are commonly used by the scientific community due to their capability to describe the natural physical processes through appropriate analytical equations (some examples are presented in [29,30,31,32]).

If a detailed description of the study area is available in terms of slope topography and physical, mechanical, and hydraulic soil properties, PB approaches can provide a high level of reliability (see [33]). Generally, the detailed reconstruction of slope topography does not represent a relevant problem; on the contrary, adequate characterization of the physical, mechanical, and hydraulic properties of soil cover is subjected to economic and practical limitations. Soils and rocks are described by parameters characterized by high variability in space both in horizontal and vertical dimensions [34]. For instance, mechanical properties show their uncertainty not only from site to site and within a given stratigraphy but also within homogeneous covers as a consequence of natural deposition processes [35]. In addition to soil property variability (random uncertainty), soil parameters are characterized by two different forms of uncertainty. The first, epistemic uncertainty, is linked to the impossibility of directly measuring a soil characteristic [36]; the second, model uncertainty, is related to the approaches used to describe a specific phenomenon [37]. Model uncertainty includes: (1) measurement uncertainty, related to systematic errors (bias) and random errors (precision); (2) statistical uncertainty, linked to limited information and influenced by the technique used; and (3) uncertainty due to the idealizations present in the physics formulation of the problem. While epistemic uncertainty can be reduced, random uncertainty cannot be eliminated [38].

Nevertheless, PB models are often used considering the quantities involved in landslide processes deprived of uncertainty, quantified by a single fixed value [39,40]; consequently, the derived predictions are expressed by a single value of the factor of safety, Fs, or the critical rainfall intensity, I_c.

On the contrary, when PB models are used with a probabilistic approach, the variability of the input quantities is modeled and the dependent variables are described as random as well.

For landslide forecasting through PB approaches, a deterministic approach, which assumes the input data without uncertainty [25,26,41,42,43,44] is less suitable than a probabilistic approach, which considers input data as random variables [45] defined through their probability density functions (pdfs) [46]. In the probabilistic approach, the safety level of the slope is given by the probability of failure (PoF), i.e., the probability associated with a value of factor of safety equal or less than 1 [47,48].

A key point in PB probabilistic modeling is represented by the definition of the theoretical pdf for the random quantities involved in the simulated physical processes. For a specific variable, the theoretical pdf must be able to reproduce its variability starting from the available measures. Typically, large variability in the physical and mechanical properties of the soil makes it quite difficult to identify a priori the theoretical pdf. In addition, when applying PB probabilistic models over large areas, the objective difficulties of having a significant number of measures to correctly assess the pdf for each variable leads to the adoption of simplified hypotheses. When the probability density functions of strength parameters cannot be determined from in situ or laboratory tests, the pdfs are estimated based on judgment, experience, or indications supported by other authors [49].

In this paper, the impact of the selection of the pdf on the PoF estimation is discussed through the application of a PB probabilistic model to a pilot study area characterized by a detailed geotechnical characterization. Available collected measures in this area are used to define the pdfs for soil shear strength properties: effective cohesion (c’) and friction angle (ϕ’). In this work, a comparison is carried out between the PoF evaluated, (i) starting from uniform pdfs for all random variables, called PoF_u, when the internal structure of the uncertainty is unknown and only the minimum and maximum values of the variable is known [50]; (ii) assuming the pdf able to consider the spatial variability of random quantities, called PoF_r. Analyses were carried out to develop a modification of the transient rainfall infiltration and grid-based regional slope stability analysis (TRIGRS) code [26] in its probabilistic version [51]. Reliability analysis was carried out using the Monte Carlo method [47].

This paper is organized as follows. In Section 2, an overview of the theoretical aspects of the approach and the equations that govern the physically based probabilistic model implemented is described; a detailed description of the study area is presented in Section 2.2, where geotechnical and hydrological assumptions considered for the reliability analysis are illustrated. Experimental settings and the results, obtained in terms of the PoF, are discussed in Section 3. Conclusions and future research developments represent the final section of the paper (Section 4).

2. Materials and Methods

2.1. Physically Based Landslide Forecasting Model

The probabilistic model, implemented in order to assess the impact of pdf definition on the PoF evaluation, represents a new extension of the probabilistic version [51] of the original TRIGRS code [26]. The PB approach, implemented in the MATLAB/Octave environment, is able to work through discretization of the study area on a regular grid, coupling a hydraulic model for the evolution of the pore water pressure during time and a mechanical model for the assessment of the temporal slope stability conditions.

In the deterministic analysis, the slope stability is expressed by the factor of safety, Fs, equal to:

$$F_s\left(Z,t\right)=\frac{tan\Phi \prime }{tan\alpha }+\frac{c^{\prime }-{\gamma }_w\psi \left(Z,t\right)tan\Phi \prime }{\gamma Zsin\alpha cos\alpha }$$

(1)

where a is the inclination grade of slope, c’ represents the effective cohesion of the soil, φ’ is the effective friction angle, ψ is the pressure head, and γ is the unit weight of the soil.

In transient flow conditions, the factor of safety varies with Z and t due to the evolution with time and space of the pressure head ψ generated by the rainfall infiltration process.

The evaluation of ψ evolution is based on the solution of the mass conservation equation for solid and liquid phases. In the particular case of fully saturated soil covers, and under 1d conditions, the equation reduces to the simple diffusion equation:

$$\frac{\partial \psi }{\partial t}=D_{\alpha }\frac{{\partial }^2\psi }{\partial Z^2}$$

(2)

where D_α is the soil coefficient of 1D consolidation, corrected for the slope inclination a.

The resolution of Equation (2) requires the definition of the initial and boundary conditions. The initial condition is defined by the position, at time t = 0, of the pressure head profile $\psi \left(Z,0\right)$, while the flow boundary conditions are defined at the layer base $Z=d_b$ and at the ground surface (Z = 0). The latter are determined starting from the time evolution of the rainfall intensity i(t). The closed-form solution for Equation (2) is described in the works by Baum et al. (2002) for the saturated case [25] and in Baum et al. (2008) for the unsaturated case [26].

In conclusion, the stability analysis is controlled by the balance of mass for the pore water, which requires the definition of: (i) saturated (θ_s) and residual (θ_r) volumetric water content, (ii) the soil stiffness (E_ed), (iii) the initial prestorm water table depth d_w and the prestorm infiltration rate parameters a_α and I_LT, (iv) the thickness of the soil cover h, and (v) the hydraulic conductivity k_s.

The Fs results depend on 12 parameters, represented in a synthetic form in Equation (3):

$$F_s=f(\alpha ,h,d_w,{\gamma }_s,c^{\prime },\phi \prime ,E_{ed},k_s,{\theta }_s,{\theta }_{r,}{a}_{\alpha ,},I_{LT})$$

(3)

In the stochastic analysis, the stability conditions are expressed by the PoF. In principle, all the 12 parameters in Equation (3) can be considered a random variable; however, the number of random variables can be reduced without affecting the forecast reliability under appropriate hypotheses. In the considered case study, the following assumptions can be made:

(1)

On the safe side of the prediction, and without detailed characterization for the unsaturated parameters, the soil cover can be considered in fully saturated conditions; thus, the balance of mass for the pore water (dψ/dt = D_α d²ψ/dZ², in which ψ is the pressure head linked to infiltration process, and D_α is the soil 1D consolidation coefficient) reduces to the simple diffusion;

(2)

Because γ_s is typically affected by low uncertainty, it can be considered constant and evaluated from the literature data;

(3)

If a high-resolution DTM is used, the slope steepness can be assumed accurate enough to be characterized by no uncertainty;

(4)

The water table depth should be monitored at different points of the study area.

Therefore, the PB probabilistic model used in this paper considers the randomness of the following quantities:

• Soil mechanical properties (c’, φ’);
• Soil saturated hydraulic conductivity k_s and soil stiffness E_ed, both considered in the infiltration problem through the coefficient of consolidation, $D_{\alpha }=\frac{\left(k_s{E}_{ed}\right)}{{\gamma }_w}$;
• Thickness of the soil cover layer, h.

The model is organized into two blocks: the first assesses the pdfs, and the second performs the Monte Carlo simulations for the reliability analysis.

The comparison is carried out in terms of PoF estimations, by considering two different pdfs: (1) the assessed pdf_r, evaluated by considering the actual variability of the random quantities; and (2) the uniform pdf_u, which ignores the internal structure of the uncertainty considering only the possible or probable range of parameters variation. In the first case, the PoF_r is obtained; in the second case, the PoF_u is estimated.

To evaluate the influence of the pdf’s theoretical distribution, the Monte Carlo simulation was used for the reliability analysis. This approach requires knowledge of the pdf distributions of the random input variables in order to generate pseudorandom numbers for each random quantity.

The exact method consists of N deterministic analyses able to define the value assumed by a random dependent variable (the PoF in this case) connected to independent random input quantities. The PoF estimation, in addition to being connected to the accuracy of the estimation of the independent variables, is strictly linked to the number of simulations performed (N) and to the number of variables considered (m). N can be expressed as:

$$N={\left(\frac{h_{\alpha /2}}{2\epsilon }\right)}^{2\cdot m}$$

(4)

where ε is given by the confidence level 1 − α, and h_α/2 is the number of standard deviation units between the mean value and that defining the assigned confidence level. High confidence levels are required to ensure a good accuracy of the results.

For unrelated variables, and for the problem analyzed, the Monte Carlo method involves the generation of 5 sequences of random numbers. The N realizations provide a sample of possible values for the safety factor, Fs. In particular, the PoF_u and PoF_r were evaluated considering the number of Monte Carlo realizations for all cases equal to N = 100,00, corresponding to a possible estimate error of about 25%.

2.2. Study Area

The study area selected for the investigation is known as Nuvole di Morra (Figure 1a), a district located in the Città di Castello municipality (northern sector of the Umbria region, central Italy). The slopes are computed directly from a rather accurate digital elevation model, TINITALY/01, which is obtained starting from separate DEMs of single administrative regions of Italy available with a 10 m cell size grid. The slopes have a steepness varying between a maximum of 30° and a minimum of 5°, with the highest values in the west and southeast parts of the area.

The area, characterized by a high susceptibility to landslides, was affected by a damaging landslide in 2005 triggered by prolonged and heavy rainfall on 10 December 2005. The movement involved an area of about 45 Ha with a soil volume of 500,000 m³. The landslide appeared to be chiefly a reactivation of pre-existing landslides; however, areas not previously affected by old landslides were also affected by the new slope movement.

The landslide (Figure 1a) can be divided into three zones: the source area (red in Figure 1a), located in the northwest sector of the area; the entrainment area (blue in Figure 1a), located in the southwest sector; and the indirect movement area (green in Figure 1a), located in the east part.

Relevant damages were caused to structures and infrastructures, and several families were transferred to other safe places (Figure 1b,c). After the event, geotechnical surveys were planned to support the structural remedial work design and implementation.

In particular, in correspondence to 8 perforation tests (S1 to S8 in Figure 1a), 4 tests of Class Q5 (undisturbed sample) and 10 tests of Class Q3 (disturbed samples of Class 3) were analyzed. In addition to laboratory tests, 2 standard penetrometric tests (SPTs), 4 geometric multichannel analysis of surface waves (MASWs) tests, 4 seismic refraction profiles, and 3 permeability tests were performed.

In this study, the available in situ measurements were integrated with the information derived by the data set of Perugia Province [52,53] to obtain a detailed geotechnical characterization that makes this area particularly suitable for the evaluation of the impact of the pdf on the model performance.

2.3. Physical and Mechanical Properties

On the basis of the documentation related to the study area, the soil analyzed can be considered homogeneous [54]. According to the geological map of the Umbria region, the outcropping lithotype is “turbidites pelitic arenaceous” (also reported in the classification described in [52,53]). The definition of the stochastic parameters (mean value and coefficient of variation (COV)) for c’ and φ’ (

Table 1. Physical and mechanical soil characterization used in the reliability analysis.

Random Variable	Symbol	Unit	pdf	Mean	COV
Effective friction angle	φ′	(deg)	Normal	30	0.20
Effective cohesion	c	(kPa)	Log-normal	5	0.25
Saturated hydraulic conductivity	ks	(m/s)	Log-normal	5 × 10−7	1
Oedometric modulus	Eed	(kPa)	Beta	1 × 10−4	0.18
Thickness of the soil cover	h	(m)	Normal	14exp(−0.07∗α)	0.20

) is derived from the integration between the regional data set information [51] and in situ test measurements, as mentioned previously in Section 3.

For the soil stiffness (E_ed) and the hydraulic conductivity (k_s), considering the few available measurements, literature data were used. In particular, for k_s, characterized by a high variability, a log-normal distribution was considered, while for E_ed, often delimited by two fixed extreme values, the beta distribution, which combines low values of physical quantity and a low probability of occurrence, was considered [34].

The thickness of the soil cover (h), for each cell of the grid, was evaluated using the following empirical formula:

H = 14exp(−0.07∗α)

(5)

where α is expressed in degrees [42], and h is a further random variable of the probabilistic model, considered normally distributed (Table 1).

The unit weight of the soil, γ_s, was assumed equal to 19 kN/m³, and the soil cover was considered homogeneous. Considering the significant rainfall events that affected the autumn of 2005, the stability conditions were precautionarily evaluated assuming the saturated conditions of the soil (Sr = 1).

In the absence of detailed information or recorded data, the initial prestorm water table depth (d_w) was set equal to 50% of h, and the steady prestorm infiltration rate (I_LT) was assumed to be negligible.

2.4. Rainfall Data

The role of the pdfs on the PoF evaluation is discussed with reference to the rainfall event that affected the Umbria region between September and December 2005. The large amount of rainfall triggered many shallow landslides; among these, Nuvole di Morra is the most relevant event.

Meteorological rainfall data for the study area were obtained through the rain gauges of the Umbria region’s monitoring network, which provides semihourly data, recorded continuously throughout the year 2005. The rainfall data recorded by the rain gauges of Città di Castello and Trestina (located near Città di Castello) were considered (Figure 2). In particular, the Morra landslide was triggered by rainfall observed between 25 and 28 November. The rainfall event considered an input in the probabilistic physically based landslide forecasting model was characterized by a total duration (t) of 32 h and cumulative rainfall of about 90 mm.

3. Results and Discussion

Figure 3 shows the spatial distribution of the safety factor FS evaluated by means of the TRIGRS code (deterministic analyses), considering the average deterministic values of the physical and mechanical soil properties reported in Table 1.

As shown in Figure 3, in the absence of rainfalls, the area is stable, with values of Fs always greater than unity. As can be expected, the areas characterized by high susceptibility to landslides correspond to the sectors characterized by the highest slope (west and southeast areas of the polygon in Figure 3).

The modified PG_TRIGRS code was applied to the study area to obtain predictions on its stability with different rainfall conditions. Using the rainfall event described in Section 2.4, the code provides the predictions shown in Figure 4. As expected, the stability conditions of the study area change during time, with the evolution of the rainfall event. To observe the effect of the rainfall on the model of spatially distributed predictions, the PoF is evaluated at four different time instants, td (td = 8 h, td = 16 h, td = 24 h, and (iv) td = t = 32 h), considering the uniform pdf for the random quantities (Figure 4).

Considering a subdivision of the PoF variability range into four classes (0 ≤ PoF < 0.15; 0.15 ≤ PoF < 0.35; 0.35 ≤ PoF < 0.5; PoF ≥ 0.5), as expected, at an early time (Figure 4a), there are no areas associated with a very high PoFu (PoF > 50%). As time increases (Figure 4b), the PoF_u starts to increase, and there is a transition from low (white cells) to high PoFu levels (orange and red pixels). Starting from time t = 24 h (Figure 4c), cells with a PoF_u variable between 15% and 50% start to concentrate in the southeast and west sector of the landscape.

At the end of the storm (Figure 4d), the stability conditions are similar to those shown for t = 24 h. In the last 8 h, the rain intensity of the event is low, and it corresponds to a cumulative rainfall of 8 mm.

In the second simulation, the pdf based on the actual measures of the study area was used. Similar results are obtained for PoF_r; in particular, at the beginning of the event (Figure 5a), areas in the western sector depicted by PoF_r > 50% are observed. The first 8 h of rain are characterized by a cumulative rain intensity of about 20 mm.

As time increases, the pixels that fall in the very high PoF_r class do not increase (Figure 5b), but a transition from a low PoF_r (white cells) to medium PoF_r values (beige and orange cells) in the southeast sector is visible (Figure 5c,d).

In this case, the stability conditions reached at t = 24 h are similar to those obtained at the end of the event.

A similar spatial PoF distribution can be observed in Figure 4 and Figure 5. The portions characterized by a higher probability of failure are distributed in the same sectors of the area; however, the pixels that as time progresses pass into higher PoF classes are different (Figure 6 and Figure 7).

It is possible to notice, as expected, that with increasing rainfall time, the areas characterized by a lower PoF decrease (in both cases of PoF_r and PoF_u), as long as the areas with bigger PoF values increase. It can also be noticed that there is an increase in the rainfall time of areas characterized by PoF_r > 50 (red classes), while the red classes of PoF_u retain small and almost constant values.

In particular, from the histograms drawn in Figure 6 and Figure 7, it can be noted that:

-

A decreasing trend with time of rainfall can be observed for gray classes: after 8 h of rainfall, the areas involved by PoF_r < 15% are about 85%, and those by PoF_u < 15% are about 70%, while after 32 h, gray PoF_r classes reach 70%, and gray PoF_u classes about 50%;

-

An increasing trend with time of rainfall can be observed for both beige (15% < PoF < 35%) and orange (35% < PoF < 50%) classes;

-

In the case of the more dangerous classes (red classes), it can be noticed that the red PoF_r portions pass from about 5% (td = 8 h) to 8% (td = 32 h), while the red PoF_u areas are negligible during all the duration of the rainfall (they reach about 1% at the end of the rainfall).

To promote a quantitative comparison between two approaches, the quantity:

$${\Delta }_E=Po{F}_u-Po{F}_r$$

(6)

was evaluated. For the sake of brevity, only the comparison in relation to the end of the event is shown (Figure 8).

The areas where the quantity ${\Delta }_E$ is negative, thus PoF_r > PoF_u, correspond to the areas characterized by the highest PoF (red class in Figure 4 and Figure 5); the portions characterized by a positive or null value of ${\Delta }_E$ correspond to areas with a lower probability of failure (white and beige pixels in Figure 4 and Figure 5).

Considering the topographic features measure after the Morra landslide, the sector in which ${\Delta }_E$ is negative corresponds to the triggering area of the landslide, and the portion in which ${\Delta }_E$ is positive corresponds to the secondary area of the sliding phenomenon.

In conclusion, the PB probabilistic model with a realistic pdf provides a higher PoF_r in the trigger zones; the PB model with a uniform pdf provides a higher PoF_u in the less susceptible sectors. The use of a uniform pdf then produces: (i) a nonprecautionary estimate of the failure conditions in the most vulnerable areas and (ii) false alarms in the stable slopes (through a PoF overestimation in areas characterized by small instability conditions).

Therefore, the link between the reliability of the results and the pdfs definition for the input random variables represents a crucial aspect in stochastic PB models, even if the attention of the scientific community has been mainly paid to the comparison between the results of the statistical and deterministic analysis [11,12,55,56].

4. Conclusions

The physical process that leads to slope failure is characterized by multiple uncertainties, which should never be considered null. The uncertainty of the quantities involved is considered in probabilistic approaches where the safety level of the slope is expressed as a random variable, quantified by the probability of the failure PoF.

The PoF evaluation is strictly related to the theoretical pdfs considered for the input quantities of the forecast models. In this work, it was shown that the definition of pdfs for the physical and mechanical properties of the soil influences the reliability analysis; this effect has been quantified with the exact method of Monte Carlo. The probability of failure was computed in relation to a particular case study where a detailed geotechnical characterization is available. This feature supports the possibility to define a pdf that is other than uniform. In addition, the small extension of the study area ensures a good compromise between the reliability of the results and the calculation times, which can prohibit results for too large a scale. In relation to the Morra area, the results show that differences between PoF_u and PoF_r are not negligible. In the portion of the area characterized by high susceptibility to landslides, the use of the pdf_u seems to lead to nonprecautionary PoF estimates (PoF_u < PoF_r). On the contrary, in areas characterized by lower susceptibility, the PoF_u appears to overestimate the conditions of potential instability (PoF_u > PoF_r). This trend occurs in relation to all the time instants considered to be independent of the rainfall intensities observed.

In order to generalize the results, the comparison in terms of the PoF on other study areas is necessary and large areas, characterized by different morphologies, will be considered in future research.

Properly defining pdfs for random parameters not only improves the reliability of the results but promotes a more informed use of physics-based probabilistic approaches.