Improving Spatial Landslide Prediction with 3D Slope Stability Analysis and Genetic Algorithm Optimization: Application to the Oltrep ò Pavese

: In this study, we compare inﬁnite slope and the three-dimensional stability analysis performed by SCOOPS 3D (software to analyze three-dimensional slope stability throughout a digital landscape). SCOOPS 3D is a model proposed by the U. S. Geological Survey (USGS), the potentialities of which have still not been investigated sufﬁciently. The comparison between inﬁnite slope and 3D slope stability analysis is carried out using the same hydrological analysis, which is performed with TRIGRS (transient rainfall inﬁltration and grid-based regional slope-stability model)—another model proposed by USGS. The SCOOPS 3D model requires deﬁnition of a series of numerical parameters that can have a signiﬁcant impact on its own performance, for a given set of physical properties. In the study, we calibrate these numerical parameters through a multi-objective optimization based on genetic algorithms to maximize the model predictability performance in terms of statistics of the receiver operating characteristics (ROC) confusion matrix. This comparison is carried out through an application on a real case study, a catchment in the Oltrep ò Pavese (Italy), in which the areas of triggered landslides were accurately monitored during an extreme rainfall on 27–28 April 2009. Results show that the SCOOPS 3D model performs better than the 1D inﬁnite slope stability analysis, as the ROC True Skill Statistic increases from 0.09 to 0.37. In comparison to other studies, we ﬁnd the 1D model performs worse, likely for the availability of less detailed geological data. On the other side, for the 3D model we ﬁnd even better results than the two other studies present to date in the scientiﬁc literature. This is to be attributed to the optimization process we proposed, which allows to have a greater gain of performance passing from the 1D to the 3D simulation, in comparison to the above-mentioned studies, where no optimization has been applied. Thus, our study contributes to improving the performances of landslide models, which still remain subject to many uncertainty factors.


Introduction
Shallow rainfall-induced landslides frequently cause human losses and substantial damages to infrastructures in many mountain and hilly regions worldwide [1]. Landslides can have devastating effects on the downstream area, especially when they evolve into debris flow [2][3][4]. As a result, many efforts have been devoted to the development of techniques and methodologies useful for the space-time prediction of rainfall-induced landslides.
Models for determining the rainfall conditions that trigger landslides can be broadly divided into two categories, namely rainfall triggering thresholds and numerical physically based models. The former are defined as the rainfall conditions the exceedance of which is likely to trigger landslides [5][6][7]. Due to the empirical nature of this approach [8][9][10][11][12][13][14][15][16], the quality and reliability of input data can affect the reliability of the prediction [17][18][19][20]. The latter simulate the hydrological and geotechnical processes responsible for the trigger [21][22][23][24][25], and can be used also for hazard mapping and thus for land planning purposes. However, their application may be hampered by the limited availability of data on soil properties [26,27]. These models are composed of two parts: a hydrological model to determine the soil response to rainfall in terms of pore pressure changes and a slope stability model to estimate the induced change in the ratio of resisting to driving forces acting on potential sliding masses [28,29].
The great majority of landslide models used infinite slope stability analysis [21,22,[30][31][32][33][34][35][36], according to which the failure of each cell is assumed to be independent of the other ones in the catchment, resulting in unstable areas that have low connectivity, which are thus quite unrealistic. Hence, there is an increasing scientific interest towards the development of software that implements three-dimensional slope stability analysis [37][38][39][40][41][42][43]. SCOOPS 3D (software to analyze 3D slope stability throughout a digital landscape) [44] belongs to the 3D models category and uses a three-dimensional (3D) approach to assess the stability of many potential landslides within a user-defined size range and considering landslide triggering as a cascade of failures of interconnected soil columns. When applied to the catchment scale, based on a digital terrain model (DTM), 3D slope stability analysis requires rather complex algorithms aimed at searching iteratively the unstable surfaces (landslides), which involve several numerical parameters. Setting these parameters in the optimal manner may be a tedious task and discourage the application of these models, or even induce performances which do not reflect their potentialities.
In this paper, we firstly aim to contribute to the investigation of the real advantages and disadvantages of 1D vs. 3D slope stability analysis at the catchment scale. Given the limitation mentioned above, we propose to use optimization techniques to define the numerical parameters of 3D slope stability models. In this context, we propose the use of genetic algorithms, a technique known for its efficiency and stability in finding global optimum solutions.
Specifically, we investigate the performance of two different models from the published literature, namely: TRIGRS (transient rainfall infiltration and grid-based regional slope-stability) [24,45], and the above mentioned SCOOPS 3D [44]. This methodology is applied on a river basin in the Oltrepò Pavese area (northern Italy), where several landslides were triggered after an extreme rainfall event (160 mm in 48 h) on 27-28 April 2009.
The last version of TRIGRS, i.e., v.2.1, provides as output a dedicated 3D format of pore water pressures suitable for the slope stability analysis in SCOOPS 3D [24]. Besides, TRIGRS also performs a 1D stability analysis, based on infinite slope assumption, which can be compared to the results of the 3D analysis through SCOOPS 3D. In our work, TRIGRS is used for the hydrological analysis. Then, the resulting pore pressure field is used first as input to the infinite slope stability model embedded into TRIGRS model itself, and then as input to SCOOPS 3D for the geomechanical analysis. Though the present is not the first application of SCOOPS 3D methodology to a real case study [46,47], here it is applied to a case study where a detailed knowledge of the occurred landslide geometries is available, due to a post-event processing of aerial images. The use of the multi-objective optimization and the availability of detailed spatial information on observed landslides enable a more in-depth assessment of the concerned models as compared to the two above-cited studies.
The paper is organized as follows. First the proposed methodology is presented in the "Materials and Methods" section, focusing on the models adopted in the models (TRIGRS, version 2.1 [24]; SCOOPS 3D, version 1.1 [44]) and on the genetic algorithm NSGAII (Nondominated Sorting Genetic Algorithm II) used in the optimization of SCOOPS 3D. Then the "study area" section describes the relevant features of the study area. Next, a comparison of the performance of the 1D and 3D approaches is presented in "Results and Discussion" section. Finally, conclusions are drawn in the last section.

Materials and Methods
The methodology adopted in this work aims to compare the performance of a 3D model vs. a 1D model in the prediction of landslide areas as a result of an intense rainfall event. This comparison requires a hydrological analysis to be performed as a preliminary step. In this work, this is done through the hydrological module of TRIGRS program, which uses a linearization of the 1D vertical infiltration Richards equation (see Section 2.1). Taking as input the pressure head fields of the hydrological analysis, the geomechanical analysis is performed using the infinite slope model implemented within TRIGRS (see Section 2.2) and the SCOOPS 3D model (see Section 2.3). While using the software SCOOPS 3D, a novel parameterization method based on the use of multi-objective optimization is proposed (see Section 2.3.1) in order to improve its prediction capability.
To test the performance of TRIGRS and SCOOPS 3D in terms of landslide prediction against observations in a real case study using the inventory of the occurred landslides geometries, the ROC analysis (receiver operating characteristic) [48] is used. ROC curves show the full picture of the relationship between true-positive rate and false-positive rate across all possible threshold values [49], namely, in this case, the link between the probability to obtain a true-positive result in the class of the observed landslides and the probability to obtain a false-positive result in the class of not occurred landslides [50].
In more detail the performance of the models can be assessed through indices based on the confusion matrix or the receiver-operating characteristics (ROCs), that is, in terms of the count of true positives (TP), true negatives (TN), false positives (FP) and false negatives (FN) (e.g., [51] As a function of the variables reported in Table 1, the three reference standard ROC indices-namely, true positive rate, false positive rate and true skill statistic-are listed below (Equations (1)-(3)): The highest performances correspond to FPR = 0 and TPR = 1, when, relatively to a given rainfall event, all observed landslide cells are equal to all unstable cells predicted by the model, i.e., model produces no false or missing predictions. A description of the methodology is presented in Figure 1. Additional remarks concerning the simulation framework are reported in Section 2.4.

Pressure Head Computation by TRIGRS Model (The Transient Rainfall Infiltration and
Grid-Based Regional Slope-Stability Model-V.2.1) TRIGRS (transient rainfall infiltration and grid-based regional slope-stability model) is a research computational model developed by U. S. Geological Survey (USGS) to describe the timing and distribution of rainfall-induced shallow landslides [52]. This is an event-based spatially-distributed model to assess the pore water pressures on a cell-by-cell basis and to outline variations in factor of safety (FS) resulting from extreme rainfall events that have durations ranging from hours to few days. Different outputs of the model can be displayed in a geographical information system (GIS) [52] and can be saved at multiple times during the simulation. The software is founded on analytical solutions for Richards' one-dimensional (1D) partial differential equation representing the vertical subsurface flow in vertically isotropic and homogeneous material. Version 1.0 of TRIGRS [45] was based on the hydrological model of Iverson [53] for a finite basal boundary bedrock depth by Savage et al. [54] whereas the version 2.0 and the version 2.1, which is used here, are based on less restrictive hypotheses than the previous version, i.e., an analytical solution to the Richards' vertical infiltration equation with a Gardner [55] retention curve. A scheme for a simulated hillslope cell is shown in Figure 2, representing the soil profile as a system consisting of two layers with a saturated zone beneath the water table that is overlain by a capillary fringe in addition to an unsaturated zone extending to the ground surface [52]. One-cell sketch for the transient rainfall infiltration and grid-based regional slope-stability model (TRIGRS) model [52].
For the simulation of the vertically rainfall infiltration at the ground surface through the unsaturated zone, the one-dimensional Richards equation is applied as follow: where Z is the vertical downward coordinate, t is the time, θ(Z, t) is the soil water content, ψ(Z, t) is the pore pressure, K(ψ) is the hydraulic conductivity and ∂ is the ground slope surface. According to Srivastava and Yeh [56], Equation (1) can be linearized using the exponential soil water retention curve proposed by Gardner [55] described by the following equations: where K s is the saturated hydraulic conductivity, α is the inverse height of the capillary rise, −ψ 0 = 1 α is the vertical height of the capillary fringe above the water table, θ r and θ s are the residual and saturated water contents, respectively.
Using Equation (5) and the computation of K(Z, t) based on a generalized form of the solution of Srivastava ad Yeh [56], the pressure head in the unsaturated zone ψ(Z, t) is obtained: where δ is the ground surface slope and α 1 = αcos 2 δ.

Infiltration, Runoff and Flow Routing
The infiltration at each cell of the domain, I, is evaluated as the sum of the precipitation rate P plus the runoff rate from the upslope cells, R u , as long as the infiltration does not exceed the hydraulic conductivity K s : When the precipitation and the supplied runoff from adjacent cells exceed infiltrability, runoff R d appears and is diverted flowing down to adjacent downslope cells. The distribution among the downslope cells is computed as follows: Further details on the TRIGRS governing equations have been reported by [24,52].

Slope Stability Analysis by TRIGRS Model (The Transient Rainfall Infiltration and
Grid-Based Regional Slope-Stability Model-V.2.1) Following Iverson (2000), the slope stability analysis in TRIGRS is assessed under the infinite slope assumption, assuming failure plans parallel to the ground surface [57]. The factor of safety FS is computed on cell-by-cell basis for an arbitrary depth Z with the following formula: where ϕ is the soil friction angle for effective stress; c is the soil cohesion for effective stress; γs and γw are the soil unit weight and the unit weight of groundwater, respectively; ψ(Z, t) is the pressure head as a function of time t and of depth Z; δ is the slope angle. Therefore, pressure head affects slope stability and, as shown by Equation (12), an increase in pressures provides a decrease in the safety factor FS. In accordance with most of the studies, FS = 1 is assumed as limiting equilibrium stadium for landslides triggering [19,21,28,53,58], therefore the failure is predicted when FS < 1 and stability holds where FS ≥ 1.
It should be pointed out that TRIGRS, as it works on cell-by-cell basis, needs GIS software for preparing the necessary input gridded data in ASCII format. This allows investigating several scenarios by varying the geotechnical and hydraulic parameters which could be assumed homogeneous in the whole study area or not.

Slope Stability Analysis by SCOOPS 3D
Model (Software to Analyze Three-Dimensional Slope Stability throughout a Digital Landscape) SCOOPS 3D is a research computational model developed by U. S. Geological Survey (USGS) [44] which uses a three-dimensional method of columns approach for analyzing slope stability throughout a digital landscape (digital elevation model, DEM). This model utilizes a considerable number of spheres to cut the terrain in order to generate several intersections, corresponding to various potential slip surfaces ( Figure 3). For each potential slip surface, SCOOPS 3D applies the 3D limit equilibrium method (LEMs) to assess the slope stability. The use of spatially distributed spheres ( Figure 3) as predefined potential slip surfaces allows considering the three-dimensional characteristics of topography. These spatially distributed spheres can be user-defined by a series of model parameters. SCOOPS 3D inspects the DEM for potential failure based on well-defined size criteria and a user-defined search grid assuming several soil layers with different properties and several groundwater inputs (dry, piezometric surface, pore-pressure ratio etc.). For a given trial surface, the three-dimensional factor of safety, FS 3D , is computed using two limit equilibrium methods: the Fellenius ordinary method [59] and Bishop's simplified method [60]. Being cell-by-cell computation method, each cell within the domain is analyzed by different trial surfaces until the lowest value of FS 3D is found. This means that, for each intersection, the FS 3D is uniformly distributed along the trial surface and the final slip surface for an individual landslide consists of many portions of spheres which represent the minimum stability for every cell. At the end of the research procedure, SCOOPS 3D generates a new terrain map with all the materials above the potential slip surface removed. Although the slope stability is related to the hillside deformations, only the balance of the involved forces is analyzed considering a rigid mass potentially sliding regardless of deformations and displacements. In general, all limit-equilibrium methods define the factor of safety as the ratio of the average shear resistance (strength), s, to the shear stress τ, required to maintain limit equilibrium along a predefined trial surface: As usual, in Equation (13), the slope is considered unstable when F s3D < 1.
Bishop's simplified method [60] is the most conventional and diffused one since it provides reliable FS 3D results that are close to more recent and rigorous LEMs [44,[61][62][63]. It is based on the assumption according to which all forces acting on vertical faces of each column can be neglected in the equilibrium equations. In SCOOPS 3D, the same assumption is made for columns according to the method proposed by Hungr [62] and the shear resistance, τ, acting on a pontential slip surface in defined by Coulomb-Terzaghi law [64] τ = c + (σ n − u)tanϕ (14) where c is the soil cohesion, ϕ is the angle of internal friction, σ n is the total normal stress acting on the failure surface and u is the pore water pressure acting on the shear surface. Therefore, the three-dimensional factor of safety is computed summing, for all columns within the potential failure surface, the FS 3D with the following formula: where, for the i,j column in a potential failure mass, A h i,j is the horizontal area of the trial slip surface, R i,j is the distance from the axis of rotation to the center of trial slip area, W i,j is the weight, u i,j is the pore-water pressure, c i,j is the cohesion of the trial slip surface and finally ϕ i,j is the angle of internal friction on trial slip surface; m α i,j is a term for part of the computation of normal force acting on the trial slip surface of the i,j column in a potential failure mass, used in Bishop's simplified method of analysis; k eq and e i,j are terms related to the option to simulate earthquake or seismic loading effects (equal to 0 in this study). More details and a complete description of the SCOOPS 3D model are reported by Reid et al. [44].

SCOOPS 3D Search Grid Configuration and Optimization by NSGAII Genetic Algorithm
As explained in the previous section, SCOOPS 3D model utilizes a considerable number of spatially distributed spheres (defined by model parameters) to generate a number of trial surfaces. The set-up of the search grid is a crucial step in order to optimize the computational speed and to ensure high performance of the model, thus it needs to be systematically and efficiently configured. The following parameters require to be set for the SCOOPS 3D simulations: The radius increments ∆ r used to explore trial surfaces. The value of ∆ r is used to increase the radius of the spheres that create a series of trial surfaces ranging in size between the user specified minimum and maximum area of the potential slip surfaces. When the potential failure mass exceeds the maximum area or intersect a DEM boundary, the generation of the spheres with radius increment ∆ r is stopped.
The minimum and the maximum elevation, z s,min and z s,max , of the search-lattice nodes with respect to the DEM elevation. The elevation of the first search-lattice node is calculated as the sum of z s,min and a multiple of z srchres (a user defined parameter usually equal to the DEM resolution). This is so that, the lattice node being used is always above the DEM.
The minimum and the maximum horizontal surface area limits, a min and a max , for potential failure masses. The area of a potential failure mass must fall in this range for a trial surface to be considered valid.
At the moment, no well-established procedures exist in the scientific literature for the parameterization of SCOOPS 3D. As an example, Tran et al. [46] analyzed several ranges of the concerned parameters until no change was found on the predicted Fs3D map. Here, the search for the optimal parameters was thus carried out by using a multi-objective approach, namely the NSGAII multi-objective genetic algorithm [65].
Genetic algorithms are a heuristic search and optimization method inspired by natural evolution. The basic concept is that evolution will find an optimal solution for the analyzed problem after a number of successive generations, similar to natural evolution. Any genetic algorithm includes the following elements: the chromosome encoding (string representations of the decisional variables in each solution), the fitness function (computa-tion evaluating the quality of the chromosome as a solution), selection (designed to use fitness to guide the evolution of chromosomes), recombination (recombination of selected chromosomes through crossover and mutation processes to form members of the offspring population) and the evolution scheme (iteration of the scheme until the stopping criteria are reached) [66]. Due to its random nature, the genetic algorithm improves the chances of finding a global solution and proves to be very efficient and stable in searching for global optimum solutions [67].
The NSGAII [65] is the second generation evolutionary multi-objective optimization genetic algorithm that improves the previous NSGA thanks to its novel features, namely: a more efficient non-domination sorting algorithm, no sharing parameter (i.e., the niche radius), and an implicitly elitist selection method that greatly aids in solving high order problems (i.e., problems with more than two objectives) [68]. NSGAII is one of the most popular multi-objective optimization algorithms and, unlike the single objective optimization approach, it simultaneously optimizes various objective functions [69]. In NSGAII, a population of individuals featuring a number of genes equal to the number of decisional variables is considered. Like for the evolution of living organisms, the population of individuals evolves through processes of mutation and crossover, until they converge, in a certain number of generations, to the set of final solutions. This set will take the shape of a Pareto front or Pareto band, in the case of two or more than two objective functions. The diversity inside the Pareto front is ensured by means of the crowding distance, a parameter that encourages the spacing between the solutions in the front. NSGAII is an elitist algorithm, in that it prevents an undominated solution to be deleted from the front. Inside the Pareto front, indeed, no solution is better than any others, since all solutions are in the optimal trade-off between the two objective functions.
Thanks to the abovementioned features, the NSGAII genetic algorithm turns to be a suitable tool in order to thoroughly find the optimal search input parameters responsible for the failure surfaces search for SCOOPS 3D.
In more detail, the two objective functions considered for the present optimization problem are the true positive rate (TPR, Equation (1)), and the false positive rate (FPR, Equation (2)), to be maximized and minimized, respectively. The aforementioned SCOOPS 3D parameters ∆ r , z s,min , z s,max , a min and a max represent, on the other hand, the decision variables ranging between the minimum and the maximum constraints defined by the user.
Inside the Pareto front, the choice of the ultimate solution can be made using several criteria. As an example, the set of parameters that yields the expected performance in terms of one of the two objectives (either FPR or TPR) can be chosen. Otherwise, the set of parameters that yields the Pareto front solution closest to the theoretical performance of the best possible model (FPR = 0 and TPR = 1) can be chosen. Finally, another criterion lies in selecting the solution farthest up from the bisector line which corresponds to a totally random prediction in a ROC plane.
In view of these consideration, a Matlab code was designed in order to implement the NSGAII optimization and to call and run SCOOPS 3D in a single routine. This avoided launching manually the model multiple times, with evident savings in terms of simulation time.

Simulation Framework
The methodology described in the previous sections ( Figure 1) can be applied as follows for a generic case study, i.e., a certain catchment in which landslides occur during an intense rain event. In a preliminary step aimed at constructing the input of the slope stability analysis, on the one hand, the hydrological properties of the soil are evaluated together with the recorded rainfall data. On the other hand, the geomechanical soil properties are evaluated and the terrain morphological analysis is performed. Then, once the pore pressure field is computed through the TRIGRS model, the soil stability analysis is performed based on the following simulations: Simulation I: 1D slope stability analysis using the TRIGRS model.
Simulation II: 3D slope stability analysis using the SCOOPS 3D model. As for the first simulation, TRIGRS requires several hydrological soil properties (total unit weight of water γ w , saturated hydraulic conductivity K s , saturated hydraulic diffusivity D 0 ) and the soil water characteristic curve for wetting of the unsaturated soil (Gardner's parameter [49] α, saturated volumetric water content θ s , and residual volumetric water content θ r ). All the parameters of TRIGRS can be derived from in-situ measurements of the soil in the analyzed case study. In more detail, field surveys together with laboratory test results have enabled to derive the hydrogeological properties needed for the abovementioned simulations. The geotechnical characterization was based on standard soil analyses carried out according to the ASTM (American Society for Testing and Materials) standards. The performed tests included an assessment of the physical parameters of the materials and triaxial tests which allowed the determination of the shear strength parameters. The hydrological properties in the study area, instead, were determined through a laboratory reconstruction of the soil water characteristic curve (SWCC) and the hydraulic conductivity function (HCF). These functions were reconstructed using undisturbed soil samples.
As for the second simulation, instead, in addition to the geotechnical and hydraulic properties, the hydrological data obtained from the TRIGRS model and the configuration of the search grid as a result of the optimization procedure through NSGAII genetic algorithm are required.
A catchment was selected with a surface area of about 2 km 2 where many shallowseated landslides occurred due to rainfall on 27-28 April 2009 (Figure 5a). Based on aerial photographs taken immediately after the event and subsequent field surveys, the shapes of the observed landslides were individuated and inventoried. Observed landslides featured very variable size, ranging from about 15 m 2 to 6300 m 2 with an average value of 500 m 2 . The distribution of terrain slope of the catchment is also shown in Figure 5b.
One crucial step in the application of the TRIGRS model is the definition of the initial conditions, mainly the finite depth of the basal boundary (that is the covering soil above of the impermeable bedrock) and the depth of the water table, both referred to the ground surface. To this end, generally steady state initial conditions are assumed. Here, it was assumed that the water table had un upper limit located approximately 0.1 m above the bedrock, accordingly to the in situ measurements collected by a monitoring station installed in the study area after the event taken as reference [71]. In more detail, some probes were installed in the soil at different depth to measure some of the main hydrological features such as the soil water content and the soil water pressure. The thickness of the basal boundary, instead, can vary as a function of many different and interplaying factors, such as the underlying lithology, the slope gradient, the hillslope curvature, the upslope contributing area, and other factors [72]. Since many applications require the availability of soil thickness measures on a dense scale, several methods for estimating the spatial patterns of soil thickness have been proposed in the scientific literature [72][73][74][75]. In this application, the empirical method proposed by Saulnier [76] has been adopted according to which the effective soil depth is defined as a function of the slope angle within the catchment [46,77]: where y min and y max are, respectively, the minimum and the maximum values of effective soil depth, and x i is the slope angle at point i. The minimum and the maximum values of soil thickness were derived from field data collected after the event investigated in this study [78] (Table 2). The location of the measuring points concerning the soil thickness used for the computation of the distribution of the soil depth within the study area is also reported in Figure 5a and the correlation between measured soil thickness points and the estimated ones is shown in Figure 6. During the mentioned landslide event dates, the Cigognola rain-gauge station recorded 160 mm of rain in 48 h (20% of the annual average amount), with a peak intensity of 22 mm/h at 09:00 p.m. on 27 April. After this maximum intensity was reached, many shallow landslides were triggered, causing one fatality, several damages to agriculture, and road blockages [78].
The geotechnical and hydraulic properties in the study area are summarized in Table 2. For a more detailed description of the geological features of the study area, please refer to Bordoni et al. [71].
Lastly, Table 3 summarizes the minimum and maximum constraints concerning the decision variables used in the optimization process of the slip-surfaces search parameters for SCOOPS 3D. These values were set based on some practical considerations available in SCOOPS 3D reference book [44]. The results of the slope stability analysis in Simulations I and II were considered at the time step of the maximum rainfall intensity, namely, 9 p.m. on 27 April 2009, when the most devastating landslides were triggered [70]. Figure 7 shows the spatially-averaged max pressure head (Panel (b)) and the amount of potential unstable cells (Panel (c)) within the study catchment during the rainfall event (Panel (a)) as a result of the 1D slope stability analysis (simulations by TRIGRS software v. 2.1). It can be noted that the number of cells corresponding to a factor of safety (FS1D) < 1 increases with total rainfall. Figure 8 shows, in terms of FPR and TPR, the results of the comparison between 1D and 3D slope stability analysis as output of the two main simulations carried out. As expected, the output of TRIGRS is made up of a single pair of values FPR and TPR. The output of SCOOPS 3D is, instead, made up of a Pareto front of optimal solutions in the trade-off between FPR and TPR.  The figure shows that, in general, the 3D approach gives better results than the 1D method in terms of model performance, demonstrating that the 3D approach is able to better describe the landslide triggering mechanisms through the assumption of more realistic slip surface geometries. Although the point in the ROC plane corresponding to 1D analysis is higher than the line of no discrimination (random prediction), it lies significantly below the points corresponding to all 3D simulations, i.e., independently from the choice of the slip-surface searching parameters. The Pareto front representative of the Simulation II features TPR values ranging between 0.54 and 0.87 and FPR values ranging between 0.28 and 0.53. In particular, if two extreme points are considered within the ROC plane, namely (a) TPR = FPR = 0 (no unstable areas predicted by the model-most extreme underestimation of landslide area) and (b) TPR = FP = 1 (no stable areas predicted by the model-most extreme overestimation of landslide areas), as shown in Figure 8, the point representing the performance of the TRIGRS model is closer to the extreme condition (a) than SCOOPS 3D model to the condition (b). This means that, by using SCOOPS 3D, it is more likely to exactly predict unstable areas than to obtain a false positive result in the class of not occurred landslides.

Results and Discussion
The wide range of performance variation associated with the different values of the slip-surface search parameters highlights the importance of their correct choice, and thus of the optimization performed. Among the solutions in the Pareto front, the best can be assumed as the one closest to the ideal performance point (FPR = 0, TPR = 1), which is FPR = 0.39 and TPR = 0.76 (true skill statistic = TPR − FPR = 0.37).
The results described above in terms of FPR and TPR are reported in the following Table 4 as well as the optimized value of the slip surfaces search parameters for SCOOPS 3D in correspondence with the ultimate solution of Simulation II. The results of the analysis conducted in terms of factor of safety FS are shown in Figure 9. In more detail, Figure 9a,b present the spatial distribution of the factor of safety as a result of both simulations carried out. According to Equations (12) and (13), the potentially unstable cells are defined by a FS < 1, namely the dark red cells of the maps. As expected, the 1D slope stability analysis, i.e., Simulation I, yield a more disconnected field of the FS, if compared to the FS 3D distribution, i.e., Simulation II. This reflects the fact that, for the infinite slope assumption, the possible failure of each cell within the study area is independent of the other ones, neglecting connectivity of the real failures.
For a more objective comparison between the two maps, the conditional distributions of the factor safety, respect to observed landslide and stable areas have been computed, and represented as cumulative empirical frequency plots in Figure 9c. As can be seen, the 3D model is able to better distinguish stable from unstable conditions, with a general higher cumulative frequency and factor of safety values ranking from 0.5 to 1.2 in the case of the observed landslides and until 2.5 in the case of no landslide occurrence. Relative to the 1D model, instead, lower cumulative frequencies are accompanied by a missing clear discrimination in terms of factor of safety between the areas interested by the observed landslides and the remaining part of the study area, as depicted by the dashed lines in Figure 9c.
Furthermore, in order to evaluate the performance of the landslide models adopted in this paper, the landslide ratio of each predicted FS class was used, namely the LRclass index [79]. LRclass index is defined as the ratio of the percentage of contained slope failure locations in each FS class to the predicted percentage of area in each class of FS category. A larger value of %LRclass corresponds to a lower over-prediction by the model [79]. Table 5 shows the results obtained concerning the LRclass analysis to allow an additional comparison with Tran et al. [46]. The unstable areas predicted are 37.5% by the 3D method and 12.8% by the 1D method. The %LRclass = 82 of the 3D approach is pretty higher than %LRclass = 64.5 of the 1D method. Figure 9. (a) Factor of safety spatial distribution as a result of 1D slope stability analysis (Simulation I); (b) factor of safety spatial distribution as a result of 3D slope stability analysis (Simulation II); (c) cumulative frequency of the factor of safety resulting by both 1D and 3D models. Overall, these results show that there is a large uncertainty in landslide prediction. In this regard, it is worthwhile to compare the performance we obtained with similar studies in the literature. In particular, relatively to the 1D model, Baum et al. [80] showed performance of the TRIGRS model for a study area north of Seattle, Washington in terms of ROC statistics, reporting a FPR = 0.08 and a TPR = 0.35 which are better than in our case (FPR = 0.12 and TPR = 0.21), likely for the availability of more detailed geological data than in our case. If we compare our results with those of Tran et al. [46], in terms of the %LRclass index for FS < 1, again we find better performances for the 1D model (85.6% vs. 64.5% in our case), while for the 3D model we obtain only slightly worse performances (87.4% vs. 82%). From another standpoint, this comparison, reveals a greater gain in performance when passing from the 1D to the 3D model. When comparing our results to those of He et al. [47], again in terms of the %LRclass index for FS < 1, the results for the 1D model are closer (71.28% vs. our 64.5%), while we obtain slightly better performances for the 3D model (82% in our case, vs. 80.16%). The greater gain in performances from 1D to the 3D that we obtain demonstrates the success of the optimization process we propose.

Conclusions
This work presented the comparison of 1D (TRIGRS) and 3D (SCOOPS 3D) slope stability models combined with the hydrological analysis performed with the former software. The applications concerned the simulation of the landslides triggered by an intense rain event, with 160 mm in 48 h, in the Oltrepò Pavese in the year 2009. Compared to other case studies in the scientific literature, this paper is the first, to the best of our knowledge, where SCOOPS 3D is tested against a detailed inventory of observed landslides with accurately measured shapes. This resulted in a more accurate evaluation of the modelling performance and in the opportunity to set up a more reliable model parameterization, obtained through the application of multi-objective optimization in the trade-off between true positive rate (TPR) and false positive rate (FPR), to be maximized and minimized respectively.
Overall, the analysis of the results has pointed out that SCOOPS 3D can enable better prediction of landslide prone areas, leading to TSS values up to 0.37, in comparison with a value of 0.09 for the 1D model TRIGRS. In comparison to other studies, we find the 1D model performs worse than in other studies, likely for the availability of less detailed geological data. On the other hand, for the 3D model we find even better results than the two other studies present to date in the scientific literature. This is to be attributed to the optimization process we proposed, which allows to have a greater gain of performance passing from the 1D to the 3D simulation, in comparison to the above-mentioned studies, where no optimization has been applied. Thus, our study contributes to improving the performances of landslide models, which still remain subject to many uncertainty factors.