A Heuristic Method to Evaluate the Effect of Soil Tillage on Slope Stability: A Pilot Case in Central Italy

: Among the various predisposing factors of rainfall-induced shallow landslides, land use is constantly evolving, being linked to human activities. Between different land uses, improper agricultural practices can have a negative impact on slope stability. Indeed, unsustainable soil tillage can modify the mechanical properties of the soils, leading to a possible increase of the instability phenomena. However, the effects of soil tillage on slope stability are poorly investigated. To address this topic, the PG_TRIGRS model (a probabilistic, geostatistic-based extension of TRIGRS) was applied to a cultivated, landslide-prone area in central Italy, thoroughly studied and periodically monitored through systematic image analysis and ﬁeld surveys. A heuristic approach was adopted to quantitatively evaluate the effect of soil tillage on the mechanical properties of the soil: after a ﬁrst run of the model with unbiased parameters, the slope stability analysis was carried out assuming several percentages of reduction of the effective soil cohesion to mimic an increasing impact of soil tillage on the strength conditions. Then, a comparison between observed landslides and the spatial distribution of the probability of failure derived from the application of PG_TRIGRS was carried out. A back analysis with contingency matrix and skill scores was adopted to search for the best compromise between correct and incorrect model outcomes. The results show that soil tillage caused a 20 to 30% reduction in soil cohesion in the analyzed area.


Introduction
Shallow landslides induced by rainfall are very common phenomena that occur in hilly and mountainous areas, causing loss of human life and environmental and economic damage [1,2].The main triggering factors of shallow landslides are represented by intense or prolonged rainfall [3][4][5], while the main predisposing factors are represented by lithology, morphology [6,7], and soil conditions-such as land cover and land use [8,9].In particular, land use is constantly evolving, and its changes affect landslide occurrence [10][11][12][13][14]. Indeed, changes in vegetation cover have an impact on the landscape diversity [15] and relevant effects on the hydrological processes and mechanical structure of the soil, with either positive or negative consequences for slope stability [16][17][18][19].Moreover, agricultural practices, on one hand, contribute positively to landscape and on the other hand, can have a negative impact on the slope stability.As an example, unsustainable agricultural practices characterized by heavy mechanization, such as soil tillage, can generate an excessive pressure on the soil, making the soil more susceptible to instability and degradation phenomena [20].
Several scientific contributions, among others, show the positive effects of vegetation cover on slope stability [21][22][23][24][25]: from a hydrological point of view, the vegetation dissipates most of the kinetic energy of the raindrops, weakening the erosion action, with a degree of interception depending on the density of the leaves and the size of the plant [26]; from a Land 2022, 11, 912 2 of 15 geo-mechanical point of view, the most important effect is represented by the mechanical reinforcement exerted by the roots [22,[27][28][29][30], consisting of an increase in the shear strength, included in the Mohr-Coulomb failure criterion as additional soil cohesion [31,32].In particular, the main contribution of roots to soil strength is related to the presence of small roots in the most superficial soil layers, which increase the compound matrix strength.Such an effect has been largely known as root reinforcement [33].Another action involves large roots which intersect the shear surface and mobilize a soil-root friction force instead of the entire tensile strength [14].Moreover, root reinforcement depends on the density of roots in soils and root diameters.As an example, grasses provide significant reinforcement to the shallower layers of soil, while woody roots of trees and shrubs provide reinforcement over a greater depth of soil through a combination of both fine, fibrous roots and coarser, woody roots [34].
Conversely, the effects of agricultural practices, and in particular of soil tillage, on slope stability conditions are poorly investigated.In this case, some studies show that the modification of soil's mechanical properties can be related to the tools used to till the land.In particular, some studies investigated how the soil characteristics vary in response to practices with the aim of assessing the efficiency of agricultural machinery, e.g., [35].As an example, tillage with a rotary paraplow can be considered as a conservation technique leaving the soil state unchanged [36].Other works, e.g., [37], show that the effect of tillage tools on the soil depends on the nature of the soil (fine-grained or coarse-grained soil).In the case of paddy soils, different laboratory experiments were conducted to determine the effect of tillage operations on the soil's physical, rheological and mechanical properties [38,39].A negative impact of the tillage operation was found in soil bulk density, while a change in the rheological behavior of paddy soil according to the variation of the moisture content was observed [38].Some authors have performed field-scale analyses aimed at evaluating the effect of different tillage techniques (conventional or conservation) on soil erosion, also in comparison with non-tilled cases [40][41][42].Overall, considerable increases in soil erosion, runoff, and sediment loss are observed in the cases with conventional tillage as compared with the cases with soil saving technologies (conservation tillage or no-tillage) [41,42].On the other hand, no significant differences in hydraulic conductivity were observed [40].Recently, Straffelini and co-authors proposed a physical modeling approach to assess runoff and soil erosion in vineyards under different soil managements and observed that continuous tillage aggravated soil erosion as compared to reference tillage, single tillage, and nectariferous [43].
To our knowledge, no articles regarding particularly the impact of soil tillage on slope stability and shallow landslide occurrence in hilly environments are currently present in the literature.Indeed, most studies have focused on steep, often terraced landscapes [see, e.g., [20] for further references] or on erosion hazards, e.g., [41,42].Quantitative measurements of the link between tillage operations and soil mechanical characteristics are not straightforward.However, there is a general consensus that tillage, plowing, and leveling generate a modification of the mechanical properties of the soils involved, leading to a possible increase in the propensity to slope failures of an area.In particular, soil tillage could lead to a decrease of the soil cohesion (e.g., due to soil disaggregation) and an increase of the friction angle (e.g., due to soil compaction).Quantitative data on this topic are rarely available: a table reporting quantitative changes in the mechanical parameters was found only in Albiero et al. (2014) [44].
In this paper, we analyze the decrease of the soil cohesion due to soil tillage, through a back-analysis approach and with the application of a probabilistic, physically-based model for the triggering of rainfall-induced landslides [45] in an agricultural environment.The choice of a probabilistic model is motivated to take into account the natural changes of the physical and mechanical properties of soils and rocks, which are characterized by high variability in space both in horizontal and vertical dimensions [46].In probabilistic approaches, 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 ≤1.Probabilistic approaches can provide a high level of reliability when a detailed description of the study area is available in terms of slope topography and physical, mechanical, and hydraulic soil properties.For landslide prediction, probabilistic approaches, which assume input data as random variables defined through their probability density functions, are more suitable than deterministic approaches, which assume the input data without uncertainty [47].
The model is applied to the Collazzone area, a cultivated area located in central Italy, characterized by a high susceptibility to landslides.The area has been the subject of several studies, e.g., [48][49][50] and is periodically monitored through systematic image analysis and on-site surveys [51,52].This allowed a preliminary quantitative assessment of the effect of the crops on the stability conditions of the area.Through a back-analysis approach and with the support of sensitivity indices, a quantitative evaluation on the effect induced by soil tillage on the mechanical properties of the soil has been provided.The contribution is divided as follows: after the introduction, an overview of the theoretical aspects of the landslide model and the method used in this study are illustrated in Section 2. The description of the study area and the database available for the mechanical soil characterization are part of this paragraph.In Section 3, after defining the geotechnical and hydrological assumptions considered for the reliability analysis, the results of the model are shown and discussed.The conclusions and future research developments represent the final section of the paper (Section 4).

Study Area and Data
The Collazzone area extends for 80 km 2 in the Perugia province, Umbria region, central Italy (Figure 1A).The geology of the area has been investigated several times [53,54]: it consists of the alternation of recent fluvial deposits along the valley bottoms, continental gravel sand and clay, travertine deposits, sandstone and marl in various percentages and thinly layered limestone.The digital elevation model (with a 20 m resolution) reveals that the territory is mainly hilly with elevations ranging between 145 and 634 m a.s.l.; the slopes that stand on the area have a gradient varying between 5 • and 50 • , with the highest values in the northern part of the area (Figure 1B). the probability associated with a value of factor of safety ≤1.Probabilistic approaches can provide a high level of reliability when a detailed description of the study area is available in terms of slope topography and physical, mechanical, and hydraulic soil properties.For landslide prediction, probabilistic approaches, which assume input data as random variables defined through their probability density functions, are more suitable than deterministic approaches, which assume the input data without uncertainty [47].
The model is applied to the Collazzone area, a cultivated area located in central Italy, characterized by a high susceptibility to landslides.The area has been the subject of several studies, e.g., [48][49][50] and is periodically monitored through systematic image analysis and on-site surveys [51,52].This allowed a preliminary quantitative assessment of the effect of the crops on the stability conditions of the area.Through a back-analysis approach and with the support of sensitivity indices, a quantitative evaluation on the effect induced by soil tillage on the mechanical properties of the soil has been provided.The contribution is divided as follows: after the introduction, an overview of the theoretical aspects of the landslide model and the method used in this study are illustrated in Section 2. The description of the study area and the database available for the mechanical soil characterization are part of this paragraph.In Section 3, after defining the geotechnical and hydrological assumptions considered for the reliability analysis, the results of the model are shown and discussed.The conclusions and future research developments represent the final section of the paper (Section 4).

Study Area and Data
The Collazzone area extends for 80 km 2 in the Perugia province, Umbria region, central Italy (Figure 1A).The geology of the area has been investigated several times [53,54]: it consists of the alternation of recent fluvial deposits along the valley bottoms, continental gravel sand and clay, travertine deposits, sandstone and marl in various percentages and thinly layered limestone.The digital elevation model (with a 20 m resolution) reveals that the territory is mainly hilly with elevations ranging between 145 and 634 m a.s.l.; the slopes that stand on the area have a gradient varying between 5° and 50°, with the highest values in the northern part of the area (Figure 1B).For the purposes of this work, defining the mechanical parameters needed in input by the model: the geotechnical classification of the Perugia province proposed by Fanelli et al., (2016) [55] was adopted, which includes an estimation of the physical and mechanical properties of the soil types constituting the near-surface cover.According to this classification, which reports the geotechnical parameters and a general description of the stratigraphy, the central part of the study area is covered mainly by clays, while the rocks (e.g., marl and travertine) are distributed in the south-eastern part of the study area (Figure 1C).
The area studied represents a high percentage of cultivated land (Figure 2) and the soil is inventoried as arable land.According to the last Agriculture Report, published on the Umbria region website (https://www.regione.umbria.it/agricoltura/statistica,last accessed on 28 April 2022), about 75% of the total agricultural area is used.According to the reference soil groups of the international soil classification [56] and to the map of the soils of the Umbria region [57], most of the area can be classified as a calcaric cambisol.
For the purposes of this work, defining the mechanical parameters needed in inpu by the model: the geotechnical classification of the Perugia province proposed by Fanel et al., (2016) [55] was adopted, which includes an estimation of the physical and mechan ical properties of the soil types constituting the near-surface cover.According to this clas sification, which reports the geotechnical parameters and a general description of the stra tigraphy, the central part of the study area is covered mainly by clays, while the rock (e.g., marl and travertine) are distributed in the south-eastern part of the study area (Fig ure 1C).
The area studied represents a high percentage of cultivated land (Figure 2) and th soil is inventoried as arable land.According to the last Agriculture Report, published o the Umbria region website (https://www.regione.umbria.it/agricoltura/statistica,last ac cessed on 28 April 2022), about 75% of the total agricultural area is used.According to th reference soil groups of the international soil classification [56] and to the map of the soil of the Umbria region [57], most of the area can be classified as a calcaric cambisol.Despite the presence of vegetation, the area presents a high landslide susceptibility typically, the slope failures are triggered chiefly by meteorological events, including in tense and prolonged rainfall and rapid snow melting [53].Using aerial and satellite im ages, Fiorucci and co-authors [51] estimated the landslide mobilization rates in the area i the period 2004-2005 and speculated that the remarkably high yearly rate of landslid mobilization observed in the area in the analyzed period might be due to the agricultura and land use practices.Considering the interesting case study for the period 1941-2005, multi-temporal landslide inventory map analyzing different sets of aerial photograph and field surveys is available for the area [53]; even now, numerous annual surveys con tinue to be made in the area and an ongoing mapping has also been made using remot sensing product such as Lidar, monoscopic and stereoscopic satellite images [51,52,58 and field surveys carried out after intense or prolonged rainfall, when images were no available.The last survey was carried out by CNR-IRPI on 20 December 2020, followin the copious rain that affected the area in the first ten days of the month.In particular, o 8 December from 03:00 A.M. to 15:00 P.M. (local time), 50 mm of cumulative rainfall fe (Figure 3), corresponding to a mean intensity of 4 mm/h (0.07 mm/min) and a peak inten sity of 9.2 mm/h (0.16 mm/min).Despite the presence of vegetation, the area presents a high landslide susceptibility, typically, the slope failures are triggered chiefly by meteorological events, including intense and prolonged rainfall and rapid snow melting [53].Using aerial and satellite images, Fiorucci and co-authors [51] estimated the landslide mobilization rates in the area in the period 2004-2005 and speculated that the remarkably high yearly rate of landslide mobilization observed in the area in the analyzed period might be due to the agricultural and land use practices.Considering the interesting case study for the period 1941-2005, a multi-temporal landslide inventory map analyzing different sets of aerial photographs and field surveys is available for the area [53]; even now, numerous annual surveys continue to be made in the area and an ongoing mapping has also been made using remote sensing product such as Lidar, monoscopic and stereoscopic satellite images [51,52,58], and field surveys carried out after intense or prolonged rainfall, when images were not available.The last survey was carried out by CNR-IRPI on 20 December 2020, following the copious rain that affected the area in the first ten days of the month.In particular, on 8 December from 03:00 A.M. to 15:00 P.M. (local time), 50 mm of cumulative rainfall fell (Figure 3), corresponding to a mean intensity of 4 mm/h (0.07 mm/min) and a peak intensity of 9.2 mm/h (0.16 mm/min).The landslide information was obtained through a reconnaissance survey of the area (20 December 2020) driving and walking along main, secondary, and farm roads [59].The investigators stopped at viewing points to check slopes where single or multiple landslides were identified and took photos of each landslide or group of landslides using a camera provided with GPS, and prepared a rapid (raw) mapping of the landslide.In the laboratory, the geolocated photographs were used to improve the location of the individual landslides and to characterize the type and the size.Landslides identified in the field and in the photographs were mapped on Google Earth.The main weakness of the reconnaissance inventory is the completeness, since from viewing points some slopes are not entirely visible, and an undetermined number of landslides may not have been identified and mapped.
During this last field survey, an event inventory was defined, including 26 shallow landslides.For the aim of this work, the landslides located in areas classified as rocks according to Fanelli et al. (2016) [55] and in uncultivated areas were excluded from the analysis; thus only 19 landslides (shown in Figure 2) were considered for modelling and analysis.As can be seen from Figures 1 and 2, these 19 landslides are located in the central portion of the area; the soils involved (Figure 1C) are clays, characterized by a cohesive resistance.Figure 4 shows an example of a shallow landslide (soil slide) mapped during the field survey.The landslide information was obtained through a reconnaissance survey of the area (20 December 2020) driving and walking along main, secondary, and farm roads [59].The investigators stopped at viewing points to check slopes where single or multiple landslides were identified and took photos of each landslide or group of landslides using a camera provided with GPS, and prepared a rapid (raw) mapping of the landslide.In the laboratory, the geolocated photographs were used to improve the location of the individual landslides and to characterize the type and the size.Landslides identified in the field and in the photographs were mapped on Google Earth.The main weakness of the reconnaissance inventory is the completeness, since from viewing points some slopes are not entirely visible, and an undetermined number of landslides may not have been identified and mapped.
During this last field survey, an event inventory was defined, including 26 shallow landslides.For the aim of this work, the landslides located in areas classified as rocks according to Fanelli et al. (2016) [55] and in uncultivated areas were excluded from the analysis; thus only 19 landslides (shown in Figure 2) were considered for modelling and analysis.As can be seen from Figures 1 and 2, these 19 landslides are located in the central portion of the area; the soils involved (Figure 1C) are clays, characterized by a cohesive resistance.Figure 4 shows an example of a shallow landslide (soil slide) mapped during the field survey.The landslide information was obtained through a reconnaissance survey of the are (20 December 2020) driving and walking along main, secondary, and farm roads [59].Th investigators stopped at viewing points to check slopes where single or multiple land slides were identified and took photos of each landslide or group of landslides using camera provided with GPS, and prepared a rapid (raw) mapping of the landslide.In th laboratory, the geolocated photographs were used to improve the location of the individ ual landslides and to characterize the type and the size.Landslides identified in the fiel and in the photographs were mapped on Google Earth.The main weakness of the recon naissance inventory is the completeness, since from viewing points some slopes are no entirely visible, and an undetermined number of landslides may not have been identifie and mapped.
During this last field survey, an event inventory was defined, including 26 shallow landslides.For the aim of this work, the landslides located in areas classified as rocks ac cording to Fanelli et al. (2016) [55] and in uncultivated areas were excluded from the anal ysis; thus only 19 landslides (shown in Figure 2) were considered for modelling and ana ysis.As can be seen from Figures 1 and 2, these 19 landslides are located in the centra portion of the area; the soils involved (Figure 1C) are clays, characterized by a cohesiv resistance.Figure 4 shows an example of a shallow landslide (soil slide) mapped durin the field survey.

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 F S along the slip surface to be evaluated with respect to translational sliding as: 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, ψ = u/γ 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, F s depends on 12 parameters, briefly represented as: in which, in addition to the quantities identified above, θ s and θ r are the saturated and residual volumetric water content, E ed the soil stiffness, d w the initial pre-storm water table depth, a α and I LT the pre-storm infiltration rate parameters; h the thickness of the soil cover and k s 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 k s ; 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.
Land 2022, 11, x FOR PEER REVIEW the seasonal transformations undergone by the slopes can be seen; these changes likely associated with agricultural practices that are supposed to reduce soil stren Figure 5. Two aerial images of the central part of the study area (red dots represent shal slides).In October, the area is devoid of any vegetation (many brown areas); in April, seve are fully vegetated (green).Background images from Google Earth.
In the absence of quantitative information and specific studies on this topic into account the impact of soil tillage in the analysis of the stability conditions of a heuristic approach was here adopted.In particular, the slope analysis was car assuming a decrease of the effective cohesion c′ equal to 10, 20, 30, 40 and 50% in a areas (cf. Figure 2B) to mimic an increasing impact of soil tillage on the geome conditions.In all model runs, the same triggering rainfall event was considered In this way, only the effect of soil tillage on the stability conditions was evaluated

Statistical Analysis
As mentioned, the random variables considered in the PG_TRIGRS approac the effective cohesion c′; (ii) the effective friction angle φ′ and (iii) the saturated h conductivity ks.The random variability is described, for each variable, by theoretic ability density function (pdf); the pdf definition was evaluated on the basis of in si urements available for different soils of the Perugia Province [55].The correlatio cient (ρ) between c′ and φ′ is assumed to be equal to −0.5 [66][67][68], while the co between ks and the soil strength parameters is assumed equal to 0. Major detai stochastic characterization of the mechanical properties of the soil are presente work of Fanelli et al., (2016) [55] and Salciarini et al. [45,69]; for further informa 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.

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 k s .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][67][68], while the correlation between k s 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.In the absence of detailed information or recorded data, the initial pre-storm water table depth (d w ) was set equal to 50% of h, and the steady pre-storm infiltration rate (I LT ) was assumed to be negligible.The soil covers were considered completely saturated (S r = 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.

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 ).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.

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.

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.
Land 2022, 11, 912 9 of 15 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).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).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 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%.
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%.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.

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 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.

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 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 c .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).

Conclusions
The paper presents a study focusing on the effect of soil tillage due to agricultural activity on slope stability conditions in a hilly environment.To this aim, a back analysis has been performed using the PG_TRIGRS model on the study area of Collazzone, central Italy.The code has already been used for the prediction of shallow landslides in the Umbria Region, but it has not been tested on areas where landslides, triggered by a specific rainfall event, have been directly observed.The proposed probabilistic approach which assumes input data as random variables defined through their probability density functions, can provide a higher level of reliability than deterministic approaches, which assume the input data without uncertainty.
For this preliminary study, we decided to perform the stability analyses under saturated conditions and we have not considered the randomness of the initial water table depth d w in order to have only three random variables in input (c , ϕ , and k s ).This choice allowed us to reduce the computational times and to better identify the effect of soil tillage on the soil's mechanical properties.However, future upgrades of the method will include the piezometric variability, the unsaturated conditions, and the soil thickness analysis.The study area selected for this purpose is of particular scientific interest because it is characterized by intense agricultural activity that greatly increases its susceptibility to landslides.The latter appears to be related to agricultural management systems which modify the physical and mechanical properties of the soil.During the year the ground is tilled and denuded in the rainiest period, while in the summer and spring months the slopes appear covered by vegetation.The continuous change in soil conditions alters the mechanical characteristics of the soil in a non-negligible way.The use of unsustainable agricultural machinery should be limited to avoid loss of cohesion in the cover soil.
From the results obtained in this study, with reference to the area analyzed, soil tillage due to agricultural activity caused a reduction in cohesion of between 20 and 30%; this estimation agrees with the results obtained in studies aimed at evaluating the effects produced by specific tillage on soil's mechanical properties [44].With reference to the results, the slope stability model is able to predict surface landslides but it classifies, as a precaution, the areas not affected by landslides.In any case, in order to safeguard the Collazzone area from shallow sliding phenomena, it is necessary to practice planned and appropriately selected agricultural techniques.
The method presented in this work is quantitative and reproducible, thus can be applied in other areas with similar environmental contexts, also for comparison with in situ and laboratory tests to refine and optimize the evaluation of the effect of soil tillage on slope stability.

Figure 1 .
Figure 1.(A) Localization of the study area (in red) within the Umbria region (blue borders); (B) slope distribution in the area; (C) geotechnical classification of soil types according to Fanelli et al., (2016) [55]; the green circles represent the landslides.

Figure 1 .
Figure 1.(A) Localization of the study area (in red) within the Umbria region (blue borders); (B) slope distribution in the area; (C) geotechnical classification of soil types according to Fanelli et al., (2016) [55]; the green circles represent the landslides.

Figure 2 .
Figure 2. Maps showing (A) the distribution of land use and (B) a picture of the cultivated an uncultivated areas in the study area.

Figure 2 .
Figure 2. Maps showing (A) the distribution of land use and (B) a picture of the cultivated and uncultivated areas in the study area.

15 Figure 3 .
Figure 3. Bar chart showing the hourly rainfall measured by the rain gauge located in Collazzone on 8 December 2020.Data provided by the administration of the Umbria region.

Figure 4 .
Figure 4. Example of shallow landslide (soil slide) in the study area.Photo taken on 20 December 2020, by CNR-IRPI.

Figure 3 .
Figure 3. Bar chart showing the hourly rainfall measured by the rain gauge located in Collazzone on 8 December 2020.Data provided by the administration of the Umbria region.

Figure 3 .
Figure 3. Bar chart showing the hourly rainfall measured by the rain gauge located in Collazzon on 8 December 2020.Data provided by the administration of the Umbria region.

Figure 4 .
Figure 4. Example of shallow landslide (soil slide) in the study area.Photo taken on 20 Decembe 2020, by CNR-IRPI.

Figure 4 .
Figure 4. Example of shallow landslide (soil slide) in the study area.Photo taken on 20 December 2020, by CNR-IRPI.

Figure 5 .
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.

Figure 8 .
Figure 8. Pixel distribution in the five classes of PoF, according to the different reductions of c′.

Figure 8 .
Figure 8. Pixel distribution in the five classes of PoF, according to the different reductions of c .

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 toFanelli et al., 2016 [55]).In the reliability analyses, rocks such as limestone, marl limestone, conglomerate, marl and travertine are not considered.

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 toFanelli et al., 2016 [55]).In the reliability analyses, rocks such as limestone, marl limestone, conglomerate, marl and travertine are not considered.

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.