*2.2. Method*

In order to evaluate the effect of soil tillage on the stability conditions of the area, the method applied in this research includes the comparison between landslides observed in situ and the PoF distribution deriving from the application of the PG\_TRIGRS model. PG\_TRIGRS (probabilistic, geostatistic-based, transient rainfall infiltration and grid-based slope stability model) is a probabilistic extension [45] of the deterministic version implemented in the original TRIGRS code, developed by Baum et al., (2008) [60], which treats each cell of the study area, subdivided into a GIS grid, as an infinite slope. Referring to deterministic slope stability analysis, TRIGRS allows the coefficient of failure FS along the slip surface to be evaluated with respect to translational sliding as:

$$F\_s(Z,t) = \frac{\tau\_f}{\tau\_m} = \frac{\tan\phi\prime}{\tan\alpha} + \frac{c' - \gamma\_w \psi(Z,t)\tan\phi\prime}{\gamma Z \sin\alpha \cos\alpha},\tag{1}$$

where *τm* is the mobilized shear stress, *τf* is the available shear strength, *c* and *φ* are the effective cohesion and friction angle of the soil, respectively, *ψ* = *<sup>u</sup>*/*γw* is pressure head, *u* is the pore water pressure, *γw* is the water unit weight, *α* is the slope angle 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 [60].

In the deterministic analysis, *Fs* depends on 12 parameters, briefly represented as:

$$F\_s = f(a, h, d\_w, \gamma\_s, c', \text{qt}, \mathbb{E}\_{\text{ed}}, k\_s, \theta\_s, \theta\_r a\_{w, \text{t}} I\_{\text{LT}}), \tag{2}$$

in which, in addition to the quantities identified above, *θs* and *θr* are the saturated and residual volumetric water content, *Eed* the soil stiffness, *dw* the initial pre-storm water table depth, *aα* and *ILT* the pre-storm infiltration rate parameters; *h* the thickness of the soil cover and *ks* the hydraulic conductivity.

The parameters present in Equation (2) are deterministic quantities, considered exact values without uncertainty, but soils and rocks are described by parameters characterized by high variability in space both in horizontal and vertical dimension [61]. 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 [62]. When the randomness of the quantities is considered, the stability conditions are expressed by the PoF and the input quantities are random variables. This is the approach followed by PG\_TRIGRS which considers random variables *c*, *φ* and *ks*; the PoF is evaluated cell-by-cell, through the probabilistic point estimate method approach. The code, validated over different areas of the Umbria Region characterized by different landslide susceptibility levels [46,63,64], has been used to evaluate the PoF distribution linked to the rainfall event described in Section 2.1. In particular, different analyses were carried out considering different scenarios of reducing the mechanical properties of the soil linked to soil tillage.

Considering the spatial distributions of slopes and soils (Figure 1), it could be observed that the landslides were not localized in areas with higher slopes. Conversely, the landslides were mostly localized in human-modified clay soils in the central part of the study area, where the land is more easily cultivated than in the eastern part of the area, where finegrained soils alternate with rocks.

These findings highlight the impact of agricultural practices on the soil's mechanical characteristics. In fact, arable land often requires conventional tillage that is able to alter the intrinsic physical properties of the soil, such as soil structure, porosity, pore-size distribution, aggregation, particle size distribution, water retention capacity and permeability but above all the effective cohesion of soil [65].

Moreover, agricultural practices induce changes on the slopes. Indeed, in the rainiest periods there are freshly tilled soils without vegetation, while in the driest periods they are covered with vegetation. To this purpose, two images of the study area gathered in

April 2018 and October 2019 and are shown in Figure 5. Moving from autumn to spring, the seasonal transformations undergone by the slopes can be seen; these changes are very likely associated with agricultural practices that are supposed to reduce soil strength.

**Figure 5.** Two aerial images of the central part of the study area (red dots represent shallow landslides). In October, the area is devoid of any vegetation (many brown areas); in April, several zones are fully vegetated (green). Background images from Google Earth.

In the absence of quantitative information and specific studies on this topic, to take into account the impact of soil tillage in the analysis of the stability conditions of the area, a heuristic approach was here adopted. In particular, the slope analysis was carried out assuming a decrease of the effective cohesion *c* equal to 10, 20, 30, 40 and 50% in all arable areas (cf. Figure 2B) to mimic an increasing impact of soil tillage on the geomechanical conditions. In all model runs, the same triggering rainfall event was considered as input. In this way, only the effect of soil tillage on the stability conditions was evaluated.

#### **3. Results and Discussion**

#### *3.1. Statistical Analysis*

As mentioned, the random variables considered in the PG\_TRIGRS approach are: (i) the effective cohesion *c*; (ii) the effective friction angle *ϕ* and (iii) the saturated hydraulic conductivity *ks*. The random variability is described, for each variable, by theoretical probability density function (pdf); the pdf definition was evaluated on the basis of in situ measurements available for different soils of the Perugia Province [55]. The correlation coefficient (*ρ*) between *c* and *ϕ* is assumed to be equal to −0.5 [66–68], while the correlation between *ks* and the soil strength parameters is assumed equal to 0. Major details of the stochastic characterization of the mechanical properties of the soil are presented in the work of Fanelli et al., (2016) [55] and Salciarini et al. [45,69]; for further information, the stochastic soil characterization is summarized in Table 1.


**Table 1.** The table summarizes the stochastic parameters (mean value and coefficient of variation, CoV) used to characterize the soils in the study area (according to Fanelli et al., 2016 [55]). In the reliability analyses, rocks such as limestone, marl limestone, conglomerate, marl and travertine are not considered.

In the absence of detailed information or recorded data, the initial pre-storm water table depth (*dw*) was set equal to 50% of h, and the steady pre-storm infiltration rate (*ILT*) was assumed to be negligible. The soil covers were considered completely saturated (*Sr* = 1) and the soil thickness did not exceed 2 m. The stability analyses were conducted under saturated conditions, as conservative assumptions, therefore evaluations about the estimated soil moisture or the groundwater level change were not considered in the study.

#### *3.2. Modelling Results*

Starting from the soil characterization shown in Table 1, the model was applied to the study area considering the rainfall event described in Section 2.1 (Figure 3). The results show that the area is mostly stable (Figure 6). The stability conditions vary considering the critical rainfall event, but the PoF values are negligible in most part of the area. A small area with PoF between 0.3 and 0.5 is located in the northwestern part, which is characterized by high slopes and Turbiditic soils (with low values of c).

**Figure 6.** The image shows the PoF distribution considering natural conditions for the soil cover and the rainfall event described in Figure 3. For the soil characterization, the mechanical properties used are reported in Table 1.

Observing the landslide distribution, we noticed that the movements occurred on the portions of the study area characterized by modest slopes; the soil mainly involved is clay (according to Fanelli et al., 2016 [55]). This evidence suggested that agricultural techniques such as soil tillage could significantly reduce the effective cohesion of the soil covers. As previously mentioned, a reduction in the average value of cohesion was considered for the arable soils in the study area. Figure 7 shows the PoF spatial distributions, evaluated at the end of the considered rainfall event, assuming reductions of *c* in arable areas respectively equal to 20 (Figure 7A), 30 (Figure 7B), 40 (Figure 7C) and 50% (Figure 7D).

**Figure 7.** Spatial distribution of PoF considering a reduction of *c* in arable areas equal to 20% (**A**), 30% (**B**), 40% (**C**), and 50% (**D**).

The 20% reduction for *c* amplified the portions of the area originally identified by PG\_TRIGRS and shown in Figure 5; however, the central sector of the area, where shallow landslides were localized, was still affected by a negligible PoF. The central sector is characterized by gentle slopes, ranging from a minimum of 3◦ to a maximum of 18◦, therefore, the triggering causes of the observed landslides were likely due to the changes in the mechanical characteristics of the soil.

As expected, higher probabilities of failure are observed considering a 30% of reduction for the effective cohesion of the soils (Figure 7B). This decrease of c produces, in the central part of the area, a large sector characterized by PoF values between 20 and 30%; and most landslides are correctly predicted by the model. The additional decrease of 40% for c provides only an increase in extension of the critical areas already identified, without significantly improving the model prediction. Finally, considering a reduction of 50% for *c*, a large part of the area becomes unstable.

Figure 8 summarizes the quantitative results obtained from the analyses. Each row represents the distribution of the cells in the various PoF classes, while for each column a different reduction of *c* has been considered. As the reduction of *c* increases, as can be expected, more pixels move towards classes with higher PoF. The pixels decrease, with reference to first class (white), is minimum passing from a reduction of about 10% to a reduction of 20% of *c*, while it becomes important considering the *c* reduction of 30 to 50%. Pixels belonging to the orange class (0.2 < PoF ≤ 0.3) also increase passing from the reduction of *c* of 20 to 30%.

**Figure 8.** Pixel distribution in the five classes of PoF, according to the different reductions of *<sup>c</sup>*.

Relating to the spatial distribution of PoF with the landslides observed in situ, the *c'* reduction corresponding to the absence of landslides in the first class (0 < PoF ≤ 0.2) is 50%. In other words, considering a reduction of *c* equal to 50%, the pixels were the landslides that have been measured and all were located in the second class (PoF > 20%). The model identified accurately the pixels characterized by landslide occurrence but incorrectly classified the other stable pixels.

#### *3.3. Reliabaility Analysis*

In order to evaluate the reliability of the proposed method, we used the standard contingencies and metrics adopted for model evaluation. For each case of *c* reduction, we calculate the number of cells with: landslides correctly hindcasted by the model (true positives, TP); landslides mapped but not hindcasted by the model (false negatives, FN); positive outcome from the model and absence of mapped landslides (false positives, FP); no landslides mapped and negative outcome from the model (true negatives, TN). Given that the outcome of the model is probabilistic and not deterministic, we set three increasing

thresholds values for PoF to define a positive and a negative outcome. In particular, we identified the values of 0.2, 0.3, 0.5 as the threshold PoF value for calculating the contingencies. As an example, in the case of the PoF threshold fixed at 0.2, we considered a positive outcome of the model to be a cell with a PoF > 0.2 and a negative outcome a cell with a PoF ≤ 2; in the case of 0.5 threshold, a positive (negative) outcome was a cell with a PoF equal or higher (lower) than 0.5. Therefore, for each case of *c'* reduction and for each PoF threshold value, we calculated the four contingencies and some skill scores: the true positive rate, TPR = TP/(TP + FN); the false positive rate, FPR = FP/(FP + TN); the Hansen–Kuipers skill score, HK = TP − FN; the predictive power, Pp = TP/(TP + FP). Detailed results are reported in Table 2. As an example, with reference to the lowest Pof threshold (0.2), the best HK value is obtained considering a 40% reduction of *<sup>c</sup>*.


**Table 2.** Values of the calculated skill scores for the five cases of *c* reduction and the three PoF threshold values. For each case of *c* reduction, the average value for each skill score is also reported. The best average values for each skill score are in bold.

Finally, we calculated the mean values of the skill scores for each case of *c* reduction, averaging the values obtained for the three PoF thresholds values (Table 2). The average TPR closest to 1 (i.e., the optimal value) was observed for a *c* reduction equal to 50%; however, this condition also generated many FP, thus, other considerations were needed. Indeed, if we considered the average values of FPR and HK, the best model performance was obtained considering a 20 and 30% reduction of *c*, respectively. Overall, the best compromise between correct and incorrect model outcomes was obtained considering a 20% reduction of *c* between 20% (best values for HK) and 30% (best value for FPR and Pp).
