Deterministic Physically-Based Distributed Models for Rainfall-Induced Shallow Landslides

1. Introduction

Among landslides, those involving the shallow soil (up to 2-3 m of depth) are most frequently induced by rainfall and are highly dangerous as no premonitory signs are present over territories. The landslide occurrence is the final step of a chain of processes that start from the rainfall infiltration and lead to the slope collapse. Therefore, a tool including process quantification is necessary to get prepared against potential landslides hazards.

Over time, different techniques have been employed to predict shallow landslides occurrences in response to rainfall. The use of physical-mathematical models, especially after the spread of computers, has gained large consensus among engineers and the scientific community, not only for applications at local or slope scale, but also at larger extensions.

To detect landslide events, two main approaches are recognized: physically-based modeling and statistical modeling. PBDMs belong to the first category, and are also referred to as process-based models. Statistical approaches are further named “data-driven models” as they use data on past landslide events. Although they represent interesting techniques for large scale hazard assessments, they highly depend on the resolution of past landslides inventories, which in most cases are not very consistent, thus leading to high levels of uncertainty propagation. Moreover, under climate change conditions, data-driven approaches have to be undertaken carefully since new, extreme climatic conditions may not be represented by past data.

Physically-based distributed models (PBDMs) have been widely used to accurately quantify landslides-related processes at different geographic scales in recent years [1]. Whether the aim is to design emergency plans to warn people or informing policy makers about consequences of extreme climatic events, PBDMs represent valid tools as they can detect landslide occurrences in advance.

On the other hand, physically-based models start from the quantification of individual processes, helping to design mitigation measures thanks to the possibility of assessing single aspects’ contribution to the overall stability. When dealing with future landslide prediction and climate change impacts on slope stability, this modeling approach is considerable noteworthy, especially when vegetation presence is considered [2].

PBDMs can be used to back-analyze real landslides to derive the required soil parameters as well. These models can be highly sensitive to the input parameters variability. However, their use for both small scale and large scale applications is approved as they can maintain a low computational effort. To avoid errors in assigning unknown soil parameters, some physical models comprise probabilistic treatments of input uncertainty or are coupled with other approaches to get probability maps of failure as outputs.

When developing new models for shallow landslide detection, decisions upon which processes to consider must be taken. Normally, PBDMs for shallow landslide prediction are composed by a hydrological module interfaced with a geotechnical module, typically based on the estimation of the slope safety factor through the limit equilibrium method. The hydrological processes involved in landslides initiation are many and including all of them in a single model can make it impossible to be practically used over large areas in reasonable computational times. Because of this complexity, in the majority of the existing models, when considering hydrological aspects in detail, the mechanical aspects are considered in simplified ways, and vice-versa. The simplified approach based on the infinite slope scheme (ISS) is adopted in the majority of PBDMs for shallow landslides . In this approach, the assumption of a planar failure geometry is consistent with shape, size, and failure mechanisms of rainfall-induced shallow landslides.

In recent years, models including antecedent hydrological information have raised importance in landslide hazard assessment [3], and including vegetation-related processes has become a central issue in shallow landslides prediction. This is because canopies and roots modify, both hydrologically and mechanically, the stability condition of slopes [4,5,6]. Moreover, the rooted portion of soil has peculiar characteristics and its behavior, together with a proper consideration of unsaturated soil mechanics, should be taken into account in shallow landslides prediction. In fact, the interest in nature based solutions (NBSs) and bioengineering techniques has raised importance in preventing landslides over territories, not only because of environmental purposes but also because of economic aspects [7,8,9]. Thus, assessing quantitatively the role of vegetation on slope stability is of importance, although a comprehensive methodology is still lacking.

In this work, a narrative review about the main physically-based distributed models developed over time is presented, highlighting which aspects or techniques were used and introduced over time by the different models. The review is restricted to the models that adopt the infinite slope scheme and are based on the limit equilibrium method. The aim is drawing a chronological memory of rainfall-induced shallow landslides modeling, also to address new research efforts in directions that were not yet explored. The paper is organized as follows. In the second section, a brief discussion about shallow landslides modeling methods is presented, stretching some conceptual aspects about hydrology and plants’ contribution. The third section focuses on models considered relevant to the aim of the review, including only models that don’t use any kind of probabilistic approach. This choice is meant to be helpful for researchers, in order to classify rainfall-induced physically-based models for shallow landslides in two different groups: “purely” physically-based distributed deterministic models (PBDDMs), i.e. the specific topic of this review, and hybrid forms. At the end, a conclusive section considers aspects that still have to be properly explored in PBDDMs. Some suggestions are also provided.

2. Shallow Landslides Modeling Methods

As already mentioned, the majority of the existing PBDDMs for rainfall-induced shallow landslides prediction are composed by two interfaced modules: one computes the hydrology, the other computes the slope stability.

The most observed mechanism for shallow landslides activation is the formation of nil or positive pore water pressure because of soil saturation. Saturated lenses develop in response to rainfall events because of infiltration and seepage processes and they can trigger landslides when enlarged enough. Water can accumulate in the soil because of the presence of an impermeable or a low-permeability layer, either during the capillary uprising process or the vertical downward seepage [10].

Because of these observations, the first slope stability models were developed based on saturated soil conditions [11], mainly considering the pore water pressure as a function of depth [12]. Under saturation conditions, soil properties and hydrological parameters can be assumed as constant and the role of matric suction (i.e negative pore water pressure) is absent. Real meteorological conditions can even be ignored in models based on saturated soil mechanics, otherwise real rainfall amounts are considered as contributing to the uprising of an existing ground water table.

Some limitations of these simplified approaches have been recognized and practical applications related to slope stability problems also involve unsaturated soil mechanics [13]. In particular, as shallow landslides have been observed to occur even under negative pore water pressure, i.e. when the soil is in partially saturated condition [14,15], it is important to consider the evolution in time of water pressure and unsaturated soil mechanics in PBDDMs.

The unsaturated soil mechanics is normally based on non-linear soil properties functions, including the constitutive relationship between soil water content and soil suction (also defined as water potential, or hydraulic head), namely the Soil Water Retention Curve (SWRC, or Soil Water Characteristic Curve, SWCC), and the one which relates soil suction or soil water content to soil hydraulic conductivity, namely the Hydraulic Conductivity Function (HCF) [16,17]. The relationships are unique for a certain pore size distribution of soil and are empirically defined. However, because of hysteretic behavior, these properties follow different patterns during the year, namely during wetting or drying phases [18]. Moreover, the presence of vegetation may change the pore size distribution because of root growth, thus leading to an evolution of the SWRC and of the soil permeability during time with respect to SWRC or HCF of a bare soil [19].

When considering unsaturated conditions, transient or stationary input rainfall can be considered for the hydrological balance computation, and water movements can be assessed through a transient analysis. The most used approach to get transient values of soil water content is the use of partial differential equations, in particular, numerical solutions of Richards’ equation for unsaturated seepage process [20,21], at different spatial domains (i.e. 1-D, 2-D or 3-D). This approach allows to consider different hydrological processes into a single balance equation, although in 3-D domain its computation can be consuming, both in terms of energy and time.

When dealing with rainfall-induced landslide hydrology, several processes should be quantified. In any case, the major water input is represented by rainfall infiltration, which enters the soil matrix and undergoes gradient and gravity-driven movements. Then, subsurface flows and redistribution processes take place.

When a numerical problem is involved, it is essential to define boundary conditions, that is defining how the considered physical system acts at the borders of its spatial and temporal domains. This allows the problem to be consistent with the physical reality.

As already mentioned, defining the antecedent condition of soil moisture prior to rainfall is fundamental to correctly simulate water movements during a precipitation event. The global soil water balance equation should be solved for a certain period to assess a consistent soil moisture condition prior to the analysis in the period of interest.

In this sense, as a key hydrological process, plant growth and evapotranspiration activity should not be neglected [22]. In fact, vegetation alters the soil water content through different mechanisms, including the absorption of water over the rooted portion of the soil, the modification of SWRCs and parameters such as the saturated hydraulic conductivity, canopies rainfall interception, and the creation of preferential flow patterns through roots and stems. This latter aspect is very difficult to be considered simultaneously with matrix flow equations, out if the model allows to consider two different domains at the same time [23]. Although preferential flow paths have different nature and originate in different ways (i.e., desiccation cracks, voids left by decaying roots etc), any of them is replicable by models only if approaches such as dual-permeability (or double porosity) are adopted [24,25].

Notwithstanding the effectiveness of Richards’ equation in soil water balance computation, this approach may not be the most efficient in specific cases. For example, when a landslide is triggered because of a wetting front depth propagation, the hydrological mechanisms can be approximated through simplified approaches for infiltration (e.g. [26,27]).

For what concerns geotechnical modules of PBDDMs, the most simplified and used approach for slope stability is the computation of a Factor of Safety (FoS). This approach is based on the limit equilibrium method and considers the relation between stabilizing and destabilizing actions on a slope. When the FoS is equal to 1, the slope is in a critical equilibrium state; as the FoS drops below 1, the entire slope is estimated as unstable.

Several methods to compute the FoS exist [28]. In a general sense, the FoS is computed over one or more potential failure surfaces, to detect under which conditions and/or at which depth the whole slope becomes unstable. PBDDMs should adopt an automatic procedure for analyzing different slip surfaces, especially when the soil strength and pressure profiles differ along with depth.

Since shallow landslides are normally translational and their length/depth ratio is generally low, the infinite slope method (ISM), which assumes a ground-parallel planar failure surface, represents the most used approach, although it has some limitations [29]. Among others, Lu and Godt developed an equation for FoS that allows the classic saturated soil mechanics theory based on the effective stress concept to be easily extended also to the unsaturated regime [30]. However, models that aim to simulate rotational movement also exist as well as models that can approximate the landslide runout [31].

It is known that vegetation contributes to soil reinforcement in different ways [9]. For what concerns plant roots, it is known that the rooting system can extend the soil shear strength either through water absorption and mechanical reinforcement, normally quantified as an additional cohesive term extending the soil cohesion. More in detail, the mechanical improvement of soil shear strength through the root network tensile strength and its interaction with soil and bedrock. Both large and fine roots contribute to the global reinforcement exerted. The overall root cohesion can be derived through different methods, namely the Wu and Waldron model (WWM), the Fiber Bundle Model (FBM) and the Root Bundle Model (RBM) [4]. The root reinforcement can be either basal or lateral when exerted, respectively, at the base of a failure plane or on the lateral sides of the landslide body. Nevertheless, in soil conditions near saturation it is not totally clear if these effects are present or not [32].

PBDDMs can be applied at slope scale (small scale) or territorial scale (large scale). In the first case, it may be possible to carry on field campaigns to get specific soil hydraulic and geotechnical parameters. In the second case, especially when the analyzed area is larger than a single slope or a small catchment, field works can become very expensive and assigning reliable soil parameters can be challenging. It is possible to refer to literature data for similar cases or, if models are capable to be provided with, some parameters uncertainty treatment should be also involved [33]. For models that do not involve probabilistic approaches, parameters values based on texture or pedotransfer functions are also normally involved if field campaign data are not available.

3. Models

SHALSTAB [34] is one of the first published PBDMs, in 1994. It computes slope stability for cohesionless soils over an infinite slope under a steady-state flow condition and it is intended to be used with ESRI-ArcGIS. It was one of the first GIS-based PBDDMs. By assuming a steady saturated flow, the model is able to calculate the rainfall amount that is necessary to trigger a landslide over a specified area, using a contour-based Digital Elevation Model (DEM) methodology. Its output is an estimation of the critical saturated soil height, computed through an hydrological model called TOPOG [35]. The model neglects the effects of the degree of saturation in the vadose zone. Because of its structure and mathematical formulations, SHALSTAB is not suitable for forecasting the timing of landslide triggering.

In 1995 the dSLAM model was published [36]. The model aim was quantifying the slope instability in steep and forested areas. With dSLAM, mechanical root reinforcement models were introduced in shallow landslides PBDDMs, while hydraulic vegetation effects were ignored. The model was designed for translational slide overlying a lithic contact. dSLAM is a contour-based model and does not account for rainfall spatial distribution.

In 1998 the SHETRAN model was extended with a shallow landslide erosion and sediment yield component [37]. It is a basin scale model that considers spatial variability of rainfall input and hydrological responses. It considers snowmelt as a triggering factor as well. The model computes the impact on sediment yields at the basin outlet. Two levels of resolution are comprised: first, a spatial assessment of shallow landslide susceptibility is carried out and a critical soil saturation condition is obtained; subsequently, a time-varying simulation based on the hydrological grid-based physical SHE model [38] is conducted at a coarser resolution; thus the model can be computationally sustainable for basin-scale simulations. The SHE model comprises evapotranspiration activity and canopies interception.

In 2000, Iverson published its model intending to assess the effects of transient rainfall on timing, rates and locations of landslides [39]. The hydrology is computed considering the soil infiltration capacity equal to the saturated hydraulic conductivity, in order to derive analytical values of pressure head, assuming a pre-existing steady-state pore pressure. The aim of the model is deriving rainfall thresholds for shallow landslides triggering, corresponding to the peak-value of pore pressure, by deriving transient pore pressure distributions that sum to the pre-existing one. Although a simplified transient analysis is included, Iverson’s model is better suitable for shallow landslides related to short duration rainfall (approximately, from 1 h to 70 h), as the model neglects lateral water transmission [40].

TRIGRS [41] is a grid-based PBDDM that aims at locating timing and size of landslides. At the time, most of the landslide prediction models were producing only susceptibility maps, not involving a complete transient and distributed analysis. In its first version, TRIGRS solves on a pixel-by-pixel basis the one-dimensional vertical version of Richards equation, assuming differentiated bare soil areas characterized by a unique homogeneous isotropic layer. The safety factor is computed in the infinite slope general framework. In general, TRIGRS is the most used PBDDM for shallow landslides at regional scale, even if no graphical interface is provided. Its use is linked to a Geographical Information System (GIS), where input data can be prepared and outputs can be visualized. Many versions of TRIGRS have been developed over time, including the vegetation effect, the condition of unsaturated regime, a parallelized version, and so on. TRIGRS provides the minimum Factor of Safety (FoS) calculated at selected time steps of a rainstorm, the pore water pressures, and the depth of minimum FoS at a certain time step. The soil is described by a SWRC. The unique consideration of 1-D water movements, although accurate under certain conditions, is not appropriate when complex topography is present over a large area. Moreover, the input DEM spatial resolution seems to strongly influence the results [42].

In the same year 2008, the SLIP model was developed [43]. It is a PBDDM that adopts a simplified hydrological approach to simulate soil water balance, assuming that the saturated portion of soil becomes larger with increasing rainfall. The soil consequently desaturates through a percolation process. To be quickly applied to large areas with low computational effort, the model avoids the need of being provided with complex hydrological formulation and approximates a transient seepage computation. In its recent versions, the SLIP model also included a simplified approach for accounting plants interception and roots cohesion [44]. SLIP represents a valid example of a simplified although effective solution for Early Warning Systems over large areas, that is an urgent need due to climate change.

In 2011, the SUSHI model by Capparelli e Versace was published [45]. It describes water movements in a bi-dimensional domain, allowing the consideration of an irregular soil stratigraphy and different soil parameters. A 2-D Richards’ equation is involved in the model, following the assumption of isotropic soil, through the adoption of a specific capillarity coefficient called C(Ψ). This coefficient represents the rate at which water is absorbed or released because of pressure head changes in the soil. A fully implicit method and the finite difference method are employed in the SUSHI model. Evapotranspiration is accounted for, by adopting a uniform root distribution and, deriving a sink transpiration term. The vegetation is modeled considering a fixed Leaf Area Index value (LAI) throughout the year.

Arnone et al. [46] published in 2011 a model called tRIBS-Landslide, then modified to include vegetation as well in 2013 [47] and in 2016 [48], turning the model into tRIBS-VEGGIE-Landslide. Their methodology is based on a Triangular Irregular Network (TIN) mesh and accounts for post-failure movement by considering certain angles as thresholds for determining if the landslide body will move up to a run-out distance. Landslide movement is assumed to follow the same flow directions evaluated by tRIBS hydrological component (the deepest descend), which is based on a transient computation of infiltration and redistribution processes. The basin morphology can be modified by landslides material, with consequent impact on most of the simulated processes. For what concerns vegetation, the model version of 2016 considers a root density distribution with depth through a Root Bundle Model, the evapotranspiration estimated through the resistivity formulation based on daily LAI that must be provided as mean annual variation. The model has not been validated in a real case study.

In 2014, Milledge et al. pointed out that the existing models were highly computationally demanding such that not practically applicable across landscapes [49]. The proposed model, called MD-STAB, simulates lateral resistances acting on landslide margins using earth pressure theory and lateral root distribution, which is modeled as an exponential function of soil depth, in a three-dimensional limit equilibrium force balance. This assumption allows the model to consider roots that cross laterally the shearing surfaces. The possibility of considering forces that act on the lateral sides of rigid blocks contradicts the infinite slope assumption, for which the inter-slice interactions are ignored. This model was extended, including derivation of root reinforcement from field measured forest stand characteristics [50]. The model ignores infiltration, soil suction and capillary rises and the groundwater level is steady and parallel to the slope surface.

In 2017, SOSlope was published by Cohen and Schwarz [51]. The model focuses on the effect of roots and soil strength on slope stability in forests. The hydrological aspects are considered through a simplified and empirical dual-porosity model. Through this approach, SOSlope can approximate the water dynamics in both the soil matrix and the preferential flow domains. The model is suitable for assessing fundamental aspects such as the role of forest structure (e.g., tree size, tree spacing), root distribution, and root mechanical properties on the triggering mechanisms of shallow landslides. SOSlope considers both lateral and basal root reinforcement. The model is particularly suitable for highly-detailed forest management purposes and outputs can be used in GIS environments.

Lizàrraga and Buscarnera in 2018 developed a model that uses suction-dependent plasticity and limit equilibrium theories to derive slope factor of safety in unsaturated soils [52]. Their model considers mechanical aspects usually not included in other models, such as suction-hardening and liquefaction potential. The aim is to simultaneously quantify the susceptibility to frictional failure and wetting-induced liquefaction in shallow slopes, in order to incorporate these considerations in regional-scale landslide hazard mapping [53], coupled with a transient hydrology computation based on Richards’ equation. Laboratory data are used to determine input parameters. The model application points out the strong interplay between infiltration mechanisms (i.e., slow or fast) and the mode and depth of slope instability. Vegetation effects are neglected.

In Figure 1 the normalized number of citations of the different models considered are represented. It can be seen that the most cited models are the Iverson model, SHALSTAB and TRIGRS. Out of them, only the dSLAM model has a normalized citation number higher than 25. These four models can thus be considered as milestones of PBDDMs knowledge and development over time.

4. Discussion and Conclusions

Based on publishing date, several models were briefly described in this narrative review trying to draw a comprehensive timeline of physically-based deterministic models for shallow landslides.

By analyzing the literature, purely deterministic approaches appear to have recently become less explored than they were in the last decades. However, it appears clear that the attention given to the stabilizing effects of vegetation has grown over years, as testified by the arising complexity involved when trying to assess these effects in a consistent way. Big challenges are still open about this topic, since most of the root reinforcement models are based on variables and parameters not easily derivable such as root architecture or root diameters. Most of the root reinforcement models have a strong empirical basis not easily applicable at large scale, particularly when different kinds of plant species cohabit [54]. New paradigms and expedients should be explored for roots’ effect quantification over large areas, and purely deterministic models can help in this task, as the processes and the related parameters can be quantified and studied singularly in a specific way. An example is provided by [55], where root-induced modifications of soil hydraulic properties (namely, the saturated hydraulic conductivity and the Soil Water Retention Curve) are included in a physically-based model. From this point of view, involving vegetation as a dynamic variable, for example accounting for growing in time and space of roots and canopy, based on real meteorological data (i.e., models based on heat accumulation of plants and thermal growth thresholds [56]), worth to be tested at different spatial scales. More studies such as those based on remote sensing, linking canopies development with root architectures, rooting depth and root spatial extensions, could help in applying root development models over large areas and vegetation growth can be more easily derived [57].

An important aspect that constitutes a strong point since the first PBDMs were published, is the interoperability with GIS environments. This aspect can help decision makers to cross different sources of information about territories. As it is difficult to have a unique, comprehensive PBDDM for all the processes that should be accounted for in landslide hazard assessment, overlaying different models outputs or different spatial information can provide more insights on a large scale. In fact, it should not be overlooked that PBDDMs pretend to give a single, specific output based on single, selected input parameters that may not represent the reality, especially when large periods are considered. Comparing different sources of information through pre- or post-processing techniques and procedures may help raise the reliability of landslide risk assessment analyses, maintaining an acceptable operational time.

In distributed models, sensitivity analysis to different DEM (or mesh) spatial resolutions should be included. It is known that accuracy changes according to the spatial resolution considered, either when computing hydrology or stability, especially when this latter is based on an infinite slope method. This aspect appears not properly considered when new models are developed, although the problem is crucial for a proper application of the model itself. In fact, spatial discretization of the domain can lead to different results [58].

Although rarely discussed when new models are developed, an important aspect for practical applications is the required running time of the different algorithms. This problem is crucial especially if the model should be used in early warning systems for civil protection purposes.

It is worth remembering that a very important expedient to overcome spatial uncertainty of parameters is including probabilistic approaches in PBDMs, thus leading to hybrid solutions. There are several valid examples in literature, such as SINMAP [59], GEOTOP-Fs [60], HIRESSS [61], SlideForMap [62] and FSLAM [63] that include these aspects at different complexity rates, but in this paper it was preferred to focus on physically-based distributed models that include only deterministic input parameters and output results.