Assessment of a Dynamic Physically Based Slope Stability Model to Evaluate Timing and Distribution of Rainfall-Induced Shallow Landslides

Abstract

Shallow landslides due to hydro-meteorological factors are one of the most common destructive geological processes, which have become more frequent in recent years due to changes in rainfall frequency and intensity. The present study assessed a dynamic, physically based slope stability model, Transient Rainfall Infiltration and Grid-Based Slope Stability Model (TRIGRS), in Idukki district, Kerala, Western Ghats. The study compared the impact of hydrogeomechanical parameters derived from two different data sets, FAO soil texture and regionally available soil texture, on the simulation of the distribution and timing of shallow landslides. For assessing the landslide distribution, 1913 landslides were compared and true positive rates (TPRs) of 68% and 60% were obtained with a nine-day rainfall period for the FAO- and regional-based data sets, respectively. However, a false positive rate (FPR) of 36% and 31% was also seen, respectively. The timing of occurrence of nine landslide events was assessed, which were triggered in the second week of June 2018. Even though the distribution of eight landslides was accurately simulated, the timing of only three events was found to be accurate. The study concludes that the model simulations using parameters derived from either of the soil texture data sets are able to identify the location of the event. However, there is a need for including a high-spatial-resolution hydrogeomechanical parameter data set to improve the timing of landslide event modeling.

1. Introduction

Shallow landslides due to hydro-meteorological factors are one of the most common destructive geological processes on the Earth’s surface and are responsible for the loss of lives, both human and livestock, every year, in mountainous regions. Although landslide occurrences are confined to hillsides and steep terrain, the devastating after-effects are not only limited to their origin but also affect the downstream areas. The cascading effects of frequent landslides, often with intense precipitation, create innumerable adversities for people and ecosystems in both highlands and lowlands [1,2,3]. Geologically, shallow, rainfall-induced landslides have an important role as an agent of shaping mountainous landscapes. While modifying the topography and renewing the terrain and ecosystems, both in the upstream and downstream areas, landslides cause the degradation of fertile land, damage to agricultural land, traffic disruption, and destruction of infrastructures, especially settlements, and have many other indirect social and economic implications [4,5,6,7,8,9,10].

Although landslides are triggered by geological (volcanic eruption, earthquakes, etc.) and anthropogenic (slope or toe excavation, slope loading, drawdown and irrigation, blast-vibrations, etc.) factors, hydro-meteorological factors still remain the major extrinsic factor [6,11]. In tectonically stable areas, e.g., southern peninsular India where this study was carried out, hydro-meteorological (intense rainfall) factors are the significant triggers, in addition to anthropogenic activities.

The current trends and projections of precipitation patterns show an increase in occurrence of high-intensity rainfall events globally, and heavy rainfalls in a short span of time will trigger a greater number of shallow landslides. As observed, a higher number of landslides have been reported globally because of the surge in the occurrence of frequent intense-localized precipitation events. The increased human presence in the mountainous regions makes the landslides and their immediate effects more severe and destructive [12,13]. With the alarming situation of climate change, population expansion, and uncontrolled exploitation of mountainous regions, tropical and subtropical mountainous areas such as Western Ghats, where this study was carried out, will be the worst affected areas by any typical hydro-geotechnical disaster such as landslides [6,14].

Shallow landslides are defined as slope failures that are generally confined to a depth of less than 3 m [15,16,17,18], and depth is determined as the distance to the failure plane from the surface. The failure plane, often termed as the slip surface, is between the soil column and bedrock or within the soil column where a sharp contrast in permeability is present [19,20]. The extent of a geographical area directly hit by landslides, especially shallow, rainfall-induced landslides, is often much smaller than the area affected by similar natural hazards such as floods, droughts, volcanic eruptions, and earthquakes. However, landslides are more frequent and recurring than many of the above-mentioned disasters and cause extensive destructions [6,8,21,22,23,24,25]. Although an individual rainfall-induced shallow landslide involves less runout volume, it has the potential to evolve into debris or earth flows and thereby gather more channel sediments and cause significant destruction along its path and in the downstream areas [3,26,27].

In India, landslides are one of the most common geohazards and are often triggered by hydrometeorological factors. India’s mountainous regions experience enormous pressure due to increased human activities along with adversities due to climate change. The impacts of geohazards such as rainfall-induced landslides have intensified. Since it is impossible to control the occurrence of rainfall-induced landslides, it is best advised to deploy comprehensive strategies to mitigate the effects. Additionally, it has to be noted that rainfall-induced landslides are site-specific and their destructive power is limited to the area of origin and the channels along the downstream areas. Adequate information on when and where a landslide may occur will help the policymakers and administrators with issuing site- or region-specific early warnings.

A number of methods are used for computing the information on when and where a landslide will occur [28]. Conventionally, the methods utilized for landslide hazard assessments or slope stability analysis to demarcate landslides can be classified into qualitative and quantitative methods [29,30,31]. The qualitative approaches can be further categorized into inventory-based mapping and heuristic methods. Inventory-based landslide assessments are one of the simplest forms of landslide hazard assessment, in which historical information of individual landslides are directly collected by visiting the affected area or through various remote sensing techniques [32,33]. Heuristic methods provide information on the degree and type of risks based on the knowledge and experience of an expert [34,35]. Because assigning weightage and ranking to causative and triggering parameters are highly subjective and localized, the expert with most experience and knowledge might produce the best results. Since each area is unique and complicated in its own ways, the ranking and weightage for one region might not produce desirable results for another region. The quantitative methods include two major categories: statistical- or stochastic- and physically based models. The former method employs mathematical models to assess slope stability based on the weighting assigned to the causative factors [36,37]. The weightages are relative, and a possibility of generalization is always present in the statistical methods. Statistical methods often work in combination with heuristic methods through bivariate, multivariate, or neural network analysis. Even though physically based models are relatively complicated compared with the rest of the methods, they are most suitable in local-scale and site-specific studies. Landslide hazard is determined using process-based slope stability models that consider the interplay of causal factors [38,39]. The results are often expressed in terms of a Factor of Safety (FoS), where lesser values indicate unstable areas.

Regardless of the complexity and data intensity of computing, physically based models are widely used to forecast rainfall-induced landslides because of the increased control over the hydrological and geotechnical aspects [40,41]. Physically based models in Geographic Information System (GIS) platforms have provided promising results in computing the timing and localization of shallow, rainfall-induced landslides at the regional scale [16,22,25,42,43,44,45,46]. These models can be employed to model slope stability ranging for areas from tens to thousands of kilometers squared. r.rotstab [47], High Resolution Slope Stability Simulator (HIRESS) [48], GEOtop-FS [49], Shallow Landslides Instability Prediction (SLIP) [50], Stability Index Mapping (SINMAP) [51], Shallow Slope Stability Model (SHALSTAB) [52], distributed Shallow Landslide Analysis Model (dSLAM) [53], and Transient Rainfall Infiltration and Grid-Based Regional Slope Stability (TRIGRS) [54] are some of the available physically based models. SINMAP, SHALSTAB, and TRIGRS are widely used in modeling the slope stability of hills. The models have been tested and validated in diverse geological, geomorphological, and climatic conditions globally. SINMAP, a probabilistic steady-state model, and SHALSTAB, a deterministic steady-state model, require geotechnical parameters and altitude information to compute slope stability and combine the Mohr–Coulomb infinite slope stability model with a steady-state hydrological model [55]. The illustration of slope failure mechanism is achieved mainly through the computation of FoS in physically based models. The FoS is the ratio of shear strength (resisting or stabilizing forces) to shear stress (driving or destabilizing forces) parallel to the surface. The SLIP model considers the saturation of part of the soil column as the major factor affecting slope failure. It assumes that the terrain is stable due to shear strength, and water flows easily into the soil column as the rainfall progresses. Large portions of the soil eventually become saturated, and the soil sliding process starts. SLIP assesses the stability of terrain using the FoS computed through the limit equilibrium method.

TRIGRS, a transient model, considers hydrological information in addition to geotechnical parameters and altitude information. Because TRIGRS considers the transient effects of rainfall on infiltration and computes the change in pore water pressure (PWP) and FoS [54], time varying analysis of slope stability is possible with TRIGRS.

TRIGRS has been previously utilized for regional-scale analysis in various studies and has been found to produce reliable results for transient rainfall events varying from hours to days [42,56,57]. In another study, importance was placed on the transient pore water pressure, which may result in the triggering of a slip [58]. In a study, based on flume tests, it was suggested that the initial porosity and saturation have an impact on the prediction capacity of the SLIP model [59]. It has been noted that another model, SHALSTAB, considers a single set of geotechnical values for a whole basin. The study emphasized the spatial discretization of these parameters to improve the predictive nature of the model [60]. One of the major limitations of SHALSTAB and SINMAP is the assumption of uniform thickness of soil, hydraulic conductivity, and steady-state shallow subsurface flow. Additionally, the slide is assumed to be translational, and the hydrological boundaries always follow the topography of the terrain [22]. While SHALSTAB and SINMAP are suggested for primary investigations of terrain stability, TRIGRS is often recommended to model specific events (such as the 2018 monsoon rainfall events in the study area) and early warning systems [45]. Although considering the soil as homogeneous has been found to over-estimate unstable pixels when using TRIGRS, it has been successfully utilized at the regional scale. As stated in one of the studies, the model parameters were set to be constant for the study area, which resulted in false positives [61]. These studies have stated the limitation of the model lies in the accuracy of the hydrogeomechanical parameters, which may need to be obtained through field investigations. These parameters are known to be spatially variable. Additionally, the lack of field data and laboratory analysis at the regional scale makes these parameters scarce. So, most of the studies are based on the inverse optimization of these parameters through the utilization of landslide inventories [56,61,62,63]. However, most of the areas lack landslide inventories with correct timing of the landslide occurrence, creating uncertainty in the optimization. The present study used the TRIGRS model [54,64,65] to compute the timing and distribution of shallow, rainfall-induced landslides in a catchment region in Western Ghats, India. In this study, the soil hydraulic parameters (SHPs) were based on two different soil texture maps available at different scales, and their influences on the occurrence of a landslide event were compared. Therefore, the impact of hydrogeomechanical parameters on the determination of FoS and PWP was quantified along hill slopes based on two soil texture data sets (FAO-based and regional), which differ in introducing soil heterogeneity to the model.

2. Materials and Methods

A catchment area was delineated to execute the TRIGRS, and the catchment area lies mostly in the district of Idukki, Kerala, India, being on the windward side of the southern part of Western Ghats. Figure 1 shows the geographical location of the study area along with the administrative boundaries and extent of Western Ghats. The present study examined the distribution of 1913 landslides as a result of rainfall events during the devastating 2018 monsoon period (Figure 2) along with timing of 9 selected shallow landslides in the study area (Figure 3).

The study area lies in one of the most landslide-prone areas in India, and both climate and topography play a significant role as causative and triggering factors [66] for the occurrences of landslides in the region. The tropical climate with intense rainfall and scorching summers have resulted in a thick overburden of poorly consolidated soil, with thickness often ranging up to 5 m (depending upon the slope) in the study area [67,68]. The combination of highly weathered bedrock, steep slopes, and heavy precipitation in the monsoon season makes the study area an ideal location for understanding, rainfall-induced shallow landslides.

The study area receives the majority of its rainfall in two monsoon periods (South West monsoon from June to September and North East monsoon from October to November). Additionally, the high-altitude and steep terrain of Western Ghats receives more rainfall than the low-lying areas due to orographic effect of Western Ghats [67,68].

Moreover, the devastating impacts of global climate changes have already been observed in the study area in the form of extreme rainfall, which are expected to increase in the coming years, which will in turn increase the frequency of slope failure in the study area [69,70,71,72,73,74]. According to the India Meteorological Department (IMD), the entire state of Kerala received 2346.6 mm of rainfall between 1 June 2018 and 19 August 2018 rather than the expected 1649.5 mm, which was about 42% above the normal rainfall rate. Furthermore, the district of Idukki received 3555.5 mm rainfall instead of the normal 1851.7 mm from 1 June 2018 to 22 August 2018 [75]. This is a 92% increase from the normal rainfall, and the IMD has identified this as large-excess deviation from normal rainfall. On 9 August, three stations in the catchment region delineated for this study received 255 mm, 254 mm, and 211 mm of rain. The intense and incessant rainfall received in the steep and undulating terrain of Western Ghats triggered multiple landslides in the study area on 9 and 11 June 2018 [76].

2.1. Data Sets

The determination of landslide occurrences through physically based models involves multiple levels of complexity. The output of models such as TRIGRS are significantly dependent on the resolution and precision of the input data sets. The following data sets, materials, tools, and software were used for understanding the timing and distribution of rainfall-induced shallow landslides in the study area. A detailed description of the data sets used for the present study is provided in the sections below.

2.1.1. Slope Failure Inventory

The entire state of Kerala received one of the most severe rainfalls in the century during June and August 2018, which triggered thousands of landslides, affecting millions of people [77,78,79,80]. The severity of the disaster led to studies on understanding the causes of slope failures in the region. The National Remote Sensing Center-Indian Space Research Organization (NRSC-ISRO), Geological Survey of India (GSI), and Kerala State Disaster Management Authority (KSDMA) generated various landslide inventories primarily by analysing pre- and post-disaster high-resolution satellite images. Later, Hao et al. improvised a landslide inventory by adding missing events and eliminating falsely detected landslides [66]. Further details on the elimination of erroneous entries and utilization of object-based image analysis for automated landslide detection can be found in [66]. The present study used the spatial information of the 1913 landslide events that occurred in the study area during 2018 monsoon rainfall based on the above-mentioned landslide inventory.

Along with the spatial information on the distribution of landslides in the study area, 9 slope failures in the study area that were triggered by the very intense rainfall between 1 and 19 June 2018 were validated for the timing of the occurrence of event. The events were documented by Geological Survey of India (GSI) under the program “Post disaster studies in Kerala” (FSP No: M4SI/NC/SR/SU-KRL/2018/21108).

Landslide 1 (Figure 3) happened in the Survey of India (SOI) toposheet number 58F/04 in Kallar Vettiyar, and the event happened on 9 June 2018. It had a dimension of 60 m × 30 m × 60 m (length × width × height), and the depth to the failure surface as observed was less than 4 m. The failure can be classified as a rainfall-induced shallow planar landslide. The area was an agricultural land with major crop being cardamom (low height plantation with no deep roots). Landslide 2 (Figure 3) was a shallow planar rainfall-induced landslide in a locality called Anachal (SOI toposheet number 58F/04) and occurred in the afternoon, around 2 PM, of 11 June, 2018. The slide had a length of 20 m, width of 30 m, and a height of 25 m. The failure surface was less than 3 m from the top, and the overlying material was very loose coarse-grained in situ soil. The slide was caused by a rise in the PWP as a result of a continuous downpour on the overlying material. A nearby four-story building collapsed as a result of the slide. Landslide 3 (Figure 3) happened on 9 June 2018 and was a shallow planar failure due to incessant rainfall. It happened in a locality called Eetticity, which lies in SOI toposheet 58F/04. The failure was a shallow planar rainfall-induced landslide with a 60 m length, 50 m width, and 45 m height. The overlying soil was planted with mixed plantations. On 9 June 2018, another landslide (landslide 4, Figure 3) happened in a locality called Kallimali (SOI toposheet number 58G/01). The landslide had a dimension of 45 m × 8 m × 40 m (length × width × height), and the overlying land was used for agriculture (pepper cultivation), and almost 1.5 acres of agricultural land on the hillslope was destroyed. According to the GSI reports, unplanned agricultural practices had caused higher infiltration of rainwater, which in turn caused a rise in the PWP and resulted in the failure. Landslide 5 (Figure 3) happened in Kallarkutty village in SOI toposheet 58G/01 on 9 June 2018, and was a shallow rotational failure triggered by rainfall. The failure was relatively small and had a length, width, and height of 8 m, 9 m, and 20 m, respectively. The same village had another landslide (Landslide 6, Figure 3) on the same day early in the morning (2 AM) in a pepper plantation. A shallow planar landslide happened in Cheruthoni on 9 June 2018, with a failure depth of 0.5 m deep and a very small run out distance of less than 2.5 m. The area lies in SOI toposheet 58B/16 (Landslide 7, Figure 3). Landslide 8 (Figure 3) occurred on 9 June 2018 in the village of Anaviratty (SOI toposheet number 58F/04). The failure was a shallow planar rainfall-induced landslide, and the overburden had thick vegetation. The depth to the failure was almost 1 m and had a very small, less than 5 m, run-out distance. Landslide 9 (Figure 3), chosen in the catchment area for the present study, happened on 9 June 2018, and was a shallow planar landslide less than 1 m deep to the failure plane. The top soil was characterized by moderate vegetation.

2.1.2. Strength and Hydraulic Parameters of the Study Area

The quantification of infiltration, runoff, flow routing, and thereby PWP and FoS was based on the strength and hydraulic parameters of the study area. The present study derived the strength and hydraulic parameters for the area from a soil texture map. The soil texture map of the area was obtained from the NASA Centre for Climate Simulation (NCCS) data portal. The Global Hybrid STATSGO/FAO Soil Texture [81,82] is provided by the NCCS data portal and is a 16-category soil texture map with 30 s resolution. Although the data are provided in 30 s spatial resolution, they were originally remapped data from the Food and Agricultural Organization of United Nations (below) 5 min data. The STATSGO/FAO data’s soil texture indices were set to the United States Geological Survey (USGS) index texture values, and further information on the original STATSGO and original FAO soil map can be found on the NCCS data portal. Three major designated soil textures were identified in the study according to the 30 s STATSGO/FAO soil texture map, and

Table 1. Hydrogeomechanical parameters (hydraulic and strength parameters) of the study area for three designated soil textures. C’ is cohesion for effective stress, ϕ’ is angle of internal friction, γ is unit weight of soil, D₀ is hydraulic diffusivity, and K_s is saturated hydraulic conductivity.

Soil Texture	Loam	Sandy Clay Loam	Clay Loam
FAO soil texture index	6	7	9
C’ (KPa)	10	29	35
ϕ’ (degree)	22.5	20	20
γ (Nm−3)	13,000	15,000	14,000
D0 (m2s−1)	0.0000094	0.0000062	0.000005
Ks (m1s−1)	0.00000453	0.00000659	0.00000272

summarizes the key hydrogeomechanical parameters (hydraulic and strength parameters) with respect to the three major corresponding soil textures of the terrain. A number of authors have examined the sensitivity of hydraulic and mechanical properties of terrain to compute the FoS and have found that physically based models perform better when the input parameters, hydraulic and mechanical, are high resolution [83]. The present study further used a detailed soil texture map from the Department of Soil Survey and Soil Conservation (DSSSC), Thiruvananthapuram, Kerala, India. The DSSSC is the state nodal agency for the conservation as well as management of soil resources for the state of Kerala. The department provides scientific databases to researchers and policymakers for the best practices on soil and land management. The detailed soil texture map was derived from the local soil texture associations provided by the DSSSC. The high-resolution hydrogeomechanical properties derived from the DSSSC soil texture maps were further utilized for deriving high-resolution FoS and then for sensitivity analyses (Section 2.2 and Section 4). Figure 4 illustrates the soil texture maps (FAO-STATSGO-derived and DSSSC-derived soil maps) of the catchment area selected for the present study. It can be seen that the FAO-STATSGO database demarcates only three soil types in the region, while eight soil types are identified in the field-based high-resolution DSSSC regional map.

2.1.3. Digital Elevation Model (DEM)

Terrain parameters, such as elevation, slope, and flow direction, have been extensively used in every landslide study at the catchment scale. The present study derived the terrain parameters from the Shuttle Radar Topography Mission (SRTM) Global DEM (Digital Object Identifier number:/10.5066/F7PR7TFT). The void filled elevation data at 1 arc-second (30 m) with a worldwide coverage was used for generating a grid-wise distribution of the altitude, slope, and flow direction. SRTM, hosted by the Endeavour space shuttle, was the first mission to use C-band spaceborne imaging radar and X-band synthetic aperture radar (X-SAR) for the global acquisition of elevation. The void-filled DEM used for the present study is hosted and distributed by the Earth Resources Observation and Science (EROS) Archive. Even though the DEM is tagged as void-filled, there were some anomalously low-altitude pixels, often referred to as depressions or pits or sinks, entirely surrounded by high-altitude pixels. Because sinks significantly interfere with the routing and flow across terrain, they have to be eliminated by filling (increasing the altitude to an extent where it allows draining off to the downhill) or breaching (lowering the edges of the dead-ends or sinks blocking the flow to allow draining off downhill). The present study used a highly efficient utility program called TauDEM, Terrain Analysis Using Digital Elevation Models [84] to remove the sinks and generate a hydrologically conditioned DEM for seamlessly deriving the terrain parameters for topographic indexing and TRIGRS.

2.1.4. Precipitation

Time-varying rainfall intensity is one of the major parameters used as input to compute the FoS and PWP in the TRIGRS model. The present study used a carefully curated synthetic rainfall history for analyzing the distribution and timing of the landslides in the study area. The IMERG and IMD data sets and a number of reports and articles from the State Disaster Management Authority (SDMA) and various authors (Section 1) were used for generating the time-varying rainfall intensity curve for a nine-day study period. The precipitation data sets from the IMD have a spatial resolution of 0.250° × 0.250° and provide daily gridded data over the Indian region [85]. Additionally, the Integrated Multi-satellitE Retrievals for GPM (IMERG) daily “final” precipitation (Level 3, version 6) is a global precipitation data set with high spatial resolution (0.10^° × 10^°) that is available half-hourly, daily, and monthly [86]. The present study used IMERG’s “final” (~3.5 months after the observation) precipitation data sets, which are satellite-gauge products using both forward and backward morphing, and the data sets are available from web portals such as the Global Precipitation Measurement (GPM) portal (https://gpm.nasa.gov/data/sources, accessed on 24 November 2022). Because the routing scheme in TRIGRS is established through the mass balance of rainfall, infiltration, and runoff, the reliability of precipitation data plays a significant role in computing slope stability.

The synthetic rainfall history prepared for the present study is shown in Figure 5. The curve shows a gradual increase in the daily rainfall from the first to the seventh day and a drop for next two days; the highest daily accumulated daily rainfall was 250 mm. Because Central Water Commission (CWC) reported heavy rainfall on 9 August 2018, and three stations in the catchment region received 255 mm, 254 mm, and 211 mm (Peermade, Munnar, and Myladumpara in the Idukki district, respectively), the rainfall history taken for this study is a possible scenario in the catchment region.

2.1.5. Spatial Distribution of Soil Thickness

Soil thickness is one of the significant input parameters in TRIGRS and plays a crucial role in determining the FoS, PWP, depth to slip surface, and hydrological response of the terrain [44,87,88,89]. Although a uniform soil depth is accepted in the model, it might produce considerably less-realistic FoSs and lead to the erroneous timing and spatial distribution of landslides. The present study obtained the depth to slip surface at 9 nine landslide locations where field observations were available through GSI.

Because the present study was carried out for larger and complex terrain, it was rather difficult to obtain (or interpolate from the available sources) the spatial distribution of soil depth in high resolution. A number of studies have been carried out in different parts of the world, and a firm linear relationship between soil depth and slope has been established [42,44,90,91]. The present study therefore assumed a linear relationship between slope and soil thickness over the study area (Equation (1)) and derived the soil depth for the study area (Figure 6).

D = −0.3437s + 3.5

(1)

where D is the soil depth at each raster cell with slope s. The slope distribution for the study area was derived from the pit-filled SRTM 30 m DEM. The depth D is in meters, and the slope s is in degrees. Waterbodies in the study area were masked as pixels corresponding to waterbodies (i.e., have the least slope), might have erroneously resulted in thicker soil.

2.2. Model Description

The dynamic physically based slope stability model used in the present study, Transient Rainfall Infiltration and Grid-based Regional Slope-Stability Model (TRIGRS) is a commonly used slope stability model for computing the timing and distribution of rainfall-induced shallow landslides [54,63,64,92,93]. The TRIGRS model consists of a runoff routing component and a slope stability component. The former computes the infiltration and subsurface flow of storm water, and the latter models the grid-based slope stability over an area of interest. Many studies have been carried out globally for different terrain with varying topographic settings, and a number of researchers have shown that the open-source model is capable of modeling rainfall-induced shallow landslides over a large region. A brief description of infiltration, runoff, flow routing, and slope stability models is given below.

2.2.1. Modeling Infiltration, Runoff, and Flow Routing in TRIGRS

Since the study area is in a low-latitude area and experiences a warm, temperate climate, it can be safely assumed that rainwater is the only form of precipitation in the study area. It is hypothesized in the TRIGRS flow routing model that the rainwater from the upslope cells flows downslope, based on the infiltrability of the cells. The infiltrability or infiltration capacity is a function of the saturated hydraulic conductivity of the medium. Additionally, runoff from a particular cell occurs only when the sum of the direct precipitation received on a cell and the accumulated runoff from the upslope cells exceeds the infiltrability of the cell.

The model computes infiltration (I) at every cell as a sum of precipitation (P) and runoff (R_u) from cells with higher elevation with an assumption that infiltration cannot exceed saturated hydraulic conductivity (K_s).

I = P + R_u [Provided (I − K_s) is negative]

(2)

Furthermore, whenever the sum of P and R_u exceeds K_s, the excess quantity, P + R_u − K_s, is channeled as runoff and routed to adjacent downslope cells, and the runoff diverted to the adjacent downslope cells (R_d) is calculated as,

R_d = P + R_u − K_s [only when (P + R_u) > K_s]

(3)

Additionally, the input water is forced to run off from every cell where the ground water table reaches the surface (implying saturated grids). The saturated cells exfiltrate water received, either in the form of rain or run off from upslope cells, and is modeled to further run off to the downslope cells.

The entire process of flow routing after the computation of run off for each cell is based on the directions of the nearest and steepest confining downslope cell. D8 flow directions [94,95] were used in the present study for flow routing cell-by-cell in the mass balance calculations, and the D8 numbering scheme was set up according to the TopoIndex (Section 2.3.1). D8 flow directions were derived in TauDEM, which output gridded cells with a numbering scheme where 1 to 8 represent east, north east, north, north west, west, south west, south, and south east, respectively. Furthermore, adequate measures were taken to avoid any possible conflicts along the grid cells adjacent to the edges of the hydrologically conditioned DEM and along the flat areas in the basin [96].

2.2.2. Modeling Slope Stability

The slope stability component in the model computes the grid-based slope stability of each cell using an infinite-slope stability analysis. A parameter called FoS, which is the ratio of resisting basal Coulomb friction to gravitationally induced downslope basal driving stress, was employed to determine the failure of infinite slope as result of rainfall. Equation (4) governs the FoS for each grid for the selected time step.

$$FoS=\frac{tan{\varphi }^{\prime }}{tan\alpha }+\frac{c^{\prime }-\psi \left(Z,t\right){\gamma }^{w}tan{\varphi }^{\prime }}{{\gamma }_s{d}_{lb}sin\alpha cos\alpha }$$

(4)

where ϕ′ is the effective angle of internal friction, c′ is the effective cohesion, ψ is the pressure head as a function of depth Z and time t, d_lb is the depth of the lower impervious boundary, and γ^w and γ_s are the unit weights of water and soil, respectively. The model was modified to include a layer of an unsaturated zone. The modified equation for obtaining the FoS in this zone is,

$$FoS=\frac{tan{\varphi }^{\prime }}{tan\alpha }+\frac{c^{\prime }-\psi \left(Z,t\right){\gamma }^w\chi tan{\varphi }^{\prime }}{{\gamma }_s{d}_{lb}sin\alpha cos\alpha }$$

(5)

where,

$$\chi =\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}$$

(6)

where the χ approximation was given by Vanapalli and Fredlund [97] and is known as the effective stress parameter. θ is the soil moisture content in the unsaturated soil zone, θ_r is the residual soil water content of the soil, and θ_s is the saturated soil water content. Figure 7 (modified from [54]) is a schematic cross-section of the hillslope where the TRIGRS model was applied. The infinite slope or the particular cell is considered stable when FoS is above 1; when the FoS falls below 1, the slope fails. Thus, the depth Z, where FoS first falls below 1, is considered the depth of landslide initiation [92].

The model considers the catchment region as digital terrain, which can be divided into zones, based on the soil texture, and three-dimensional cell grids. The user or modeler can introduce more control by feeding each zone topographic parameters (elevation and slope), hydraulic properties (flow direction, diffusivity, infiltration rate, depth to ground water, saturated water content, residual water content, saturated hydraulic conductivity, and inverse of the height of capillary rise), and strength parameters (cohesion, friction angle, and weight of the soil). Furthermore, control can be achieved by providing elevation, slope, flow direction, depth to ground water, and depth to bedrock grid-wise. The TRIGRS version used for the present study did not have the versatility to provide strength and hydraulic parameters grid-wise. The present study employed version 2.1 [54] of TRIGRS, which is based on version 1.0 [64], but with fewer restrictive hypotheses. The model is capable of performing on saturated as well as unsaturated media and is capable of capturing unstable cells, which are effectively landslide locations, in the study area by considering precipitation and the corresponding changes in the hydrogeomechanical parameters.

The execution time required for TRIGRS has an exponential negative relationship with the number of pixels, grids, and cells in the input data sets. A larger number of pixels can be expected whenever the model is executed with high-resolution input data sets or over a large area. Because the present study was carried out for a relatively larger catchment area with a higher number of grids in the input data sets, multiple model executions became time consuming. In order to overcome this issue, the study employed TRIGRS MPI [65] for the parallel implementation of TRIGRS by utilizing the maximum computational capacity. The TRIGRS MPI significantly reduced the time consumed for each run, provided more flexibility and freedom while changing the input parameter sets for multiple executions of the model.

2.3. Application of the Model

The present study assessed the distribution of slope failures in the catchment region at 1913 landslide sites along with computing the fall in FoS and rise in PWP at 9 landslide sites in the study area with respect to a 9-day rainfall period. The open-source model, TRIGRS, is capable of considering topographic parameters (elevation and slope), hydraulic properties (flow direction, permeability, diffusivity, conductivity, infiltration rate, and depth to ground water), and strength parameters (cohesion, friction angle, and weight of the soil) while computing the FoS.

A brief description of the significant steps involved in this work is described in the following sections.

2.3.1. Topographic Indexing for Flow Routing

The execution of the model TRIGRS is carried out through two individual programs: TopoIndex and the main TRIGRS program. The former, as the name suggests, is used for topographic indexing; the latter, the main program, is used for flow routing and runoff calculations. TopoIndex prepares a list and a grid of downslope receptor cells (TIdsneiList and TIdscelGrid), a list of index numbers corresponding to each cell number (TIcelindxList), a list of downslope cells for which nonzero weighting factors have been computed (TIdscelList), and a list of weighting factors for downslope receptor cells (TIwfactorList). The outputs of the TopoIndex are generated in ASCII, and text formats and each output’s name start with what is given above in the brackets.

A hydrologically conditioned DEM and flow direction are the input files required for generating the essentials data sets for TRIGRS. Section 2.1.3 briefly describes the hydrological conditioning of the DEM. The present study used the D8 algorithm to derive the flow direction through TauDEM (Section 2.1.3). Although TopoIndex uses the D8 flow direction, the numeric coding is different than that of the TauDEM-derived flow direction or ESRI flow direction. Figure 8 shows the designated codes in both the flow direction along with the geographical direction. Although the flow direction for TopoIndex can be generated through with D8 and D-infinity methods, the present study used the TauDEM-derived D8 flow direction scheme for generating the outputs.

Figure 9 shows a schematic diagram of the study, which depicts the derivation of the input data sets, application of the model, and validation of the output (grid-based FoS) from the hydrogeomechanical properties through FAO-STATSGO and high-resolution hydrogeomechanical properties through the field-based DSSSC against the landslide inventory. The landslide inventory is based on the landslides that occurred in the catchment region during the 2008 monsoon season.

2.3.2. TRIGRS Initialization for Slope Stability

The present study divided the entire study area into three zones based on the FAO soil texture and eight zones based on the regional soil texture map (Figure 4). The hydrological and strength parameters were based on both soil texture maps. The hydrological and strength parameters, such as cohesion (c′), internal friction (ϕ′), unit weight of soil (γ), hydraulic diffusivity (D₀), saturated hydraulic conductivity (K_s), and saturated and residual water content, were derived solely based on the soil texture map. The values of these parameters were obtained from the literature based on the soil type. The present study divided the entire catchment region into 30 m × 30 m cells, and the FoS was computed for each grid with a time step of 1 day for a total period of nine days. An FoS value of 1 was considered as the threshold for slope failure in the catchment region. Thereby, a value greater than 1 was taken as a measure to classify the pixel as stable (non-landslide) and a value less than 1 was taken as a measure to classify the pixel as unstable (landslide) in the catchment region. The computed FoS for every pixel in the study area for the entire nine-day period was taken as an indicator of the initiation of slope failure.

3. Results

As discussed in Section 2.1.1, the present study examined the landslides in a landslide-prone catchment region in Western Ghats. The shallow landslides were triggered by widespread heavy and incessant rainfall that occurred in the first weeks of June 2018. Details of the hydrogeomechanical parameters and terrain parameters used for computing the FoS and PWP in the physically based model are provided in Section 2.1.2 and Section 2.1.3. The rainfall rate fed into the TRIGRS model is described Section 2.1.4. The distribution of the landslide pixels and the changes as the rainfall continued are described in Section 3.2, and the timing of the triggering of the nine landslides as a result of the rainfall in the study are is described in Section 3.3 and Section 3.4.

3.1. Impact of Varying the Number of Soil Zones

Even though the number of soil zones in the studied catchment varied for the two runs, being three and eight, the properties were still based on soil texture. It can be seen from Figure 10 that the true identifications of the landslide events in both runs were comparable. The model in both the cases could identify 1524 and 1323 landslides of the 1913 landslides. Additionally, a true positive rate (TPR) of 68% and 60% was obtained with a nine-day rainfall period for the two derived soil texture data sets, respectively. A false positive rate (FPR) of 36% and 31% was also seen, respectively. The total accuracy was 68%. Although previous studies [83] showed that the hydrogeomechanical parameters have an impact on the simulations of FoS and PWP, it was seen here that the properties just based on soil texture did not impact the model result. It was also suggested that these properties may need to be derived from the other sources such as Earth Observation (EO) data sets, which can result in a better simulation of the PWP [98].

3.2. Change in FoS and PWP in the Catchment Region

In the study area, the FoS corresponding to majority of the grids showed a significant decrease during the nine-day rainfall period. Furthermore, over time, a higher number of grids started showing an FoS below the threshold value of one, and the corresponding time was taken as the initiation of the landslide. Along with the FoS, the PWP was computed for each grid to further illustrate the interconnection among the FoS, landslide initialization, and rainfall history. Figure 11 shows the fall in the FoS and rise in the PWP for the study area for the period of nine days. It is depicted in the figure that the FoS was decreasing for the majority of the pixels, and new areas could be demarcated as landslide pixels as time progressed. The FoS was calculated for each pixel at the end of a 24 h rainfall period for nine days in Figure 11 (left). The same was repeated for the PWP and is plotted in Figure 11 (right). It is clear from the FoS–PWP through the rainfall rate plot that a number of pixels fell below the threshold FoS value and turned into unstable areas.

Furthermore, Figure 11 shows the inverse relationship between the FoS and PWP in the study area. It can be observed from the figure that when PWP increased as rainfall continued, the FoS significantly decreased. The blue region in the plots represents the relatively stable areas (lower PWP and higher FoS), and the red region represents relatively lesser stable areas (higher PWP and lower FoS). For better understanding, the usage of model with regard to landslide initiation timing, the nine landslides were chosen for the case study. The nine landslides chosen for the study are located in the black rectangle in the plots.

The plot in Figure 8 is a detailed illustration of the FoS–PWP in and around the nine landslide sites chosen in the study area. As discussed above, the rise in the PWP and fall in the FoS were computed and plotted for the nine-day rainfall period. The black dots represent the locations of the nine landslides chosen for the study area. It can be seen from Figure 11 and Figure 12 that there are false positive regions, where the FoS drops below one. This shows that the geohydrological parameters that were governing the hydrological simulation needed to be further improved so that the partition of water among the soil layers could be better simulated. The χ approximation from Equation (6), utilized in Equation (5) of the FoS, is based on the SHPs, saturated water content, and residual water content. It can be seen that the properties impact both the PWP and FoS. These parameters were literature-based and dependent on the soil texture. One of the limitations is that these properties were not available at a 30 m resolution. Again, this shows that these properties may have to be identified at finer resolution through satellite data that have higher spatial coverage compared with the field analysis, which may be not possible in remote and inaccessible areas. However, the landslide events were identified, which can still act as a warning for susceptible zones.

3.3. Fall in FoS

Figure 11 and Figure 12 show the distribution of the landslides in the catchment region and their changes as the rainfall progresses. The 8 landslides were triggered on 9 June 2018 and one landslide was triggered on 11 June 2018. The rainfall event used for computing the eight landslides which were triggered on 9 June 2018 started on 3 June and ended on 11 June 2018. Seventh day of the rainfall event (9 June) was the day on which eight landslides happened. Similarly, a rainfall event started on 5 June and ended on 11 June 2018 was used for computing the FoS at one landslide which was triggered on 11 June 2018.

Figure 13 shows the gradual fall in FoS at nine landslide sites with respect to the nine days rainfall event. A vertical dotted line is added to each plot for easy representation of day on which the landslide occurred (9 June in the case of Landslide 1, 3, 4, 5, 6, 7, and 8 and 11 June 2018 in the case of Landslide 2). As discussed above, the present study considered 1 as a threshold value for FoS and whenever a pixel’s FoS fall below 1 is considered as a landslide and the corresponding time is recorded as the landslide initiation time.

Landslide 1 (Kallar Vattiyar), Landslide 3 (Eetticity), Landslide 4 (Kallimali), Landslide 5 (Kallarkutty I), Landslide 6 (Kallarkutty II), Landslide 7 (Cheruthoni), Landslide 8 (Anaviratty) and Landslide 9 (Ambazhachal) were triggered on 9 June 2018. The FoS computed at these landslide locations on 9 June are 1.15, 6.87, 0.90, 0.89, 1.08, 1.08, 1.11, and 1.22 (FoS corresponding to dotted line on Figure 13). Additionally, Landslide 2 (Anachal), was triggered on 11 June and the FoS computed as per TRIGRS on the same day is 0.93 (FoS corresponding to dotted line on Figure 13). According to the FoS computed Figure 13, the model could indicate the slope failures at three landslide sites out of nine chosen landslides in the study area. Additionally, as the time progresses it can be observed that the model rightly shows the eight out of 9 landslides failures in the study area.

3.4. Rise in PWP

As part of the study, the PWP was calculated for the catchment region. Figure 14 plots the PWP extracted at the nine landslide sites in the study area against the time-varying rainfall rate. The plot shows a near linear rise in the PWP with respect to the rainfall. Furthermore, the PWP and FoS showed an inverse relationship at the nine landslide sites (Figure 13 and Figure 14). As in the case of rainfall-induced shallow landslides anywhere else, an increase in the PWP indeed reduced the FoS at the nine landslide sites in the study area.

4. Discussion

The timing and distribution of shallow landslides were modeled in a catchment region (Figure 1) in Western Ghats through TRIGRS. The hydrogeomechanical properties used in TRIGRS to model the FoS and PWP were derived through two different sources: FAO-STATSGO and a regional soil texture map from DSSSC. Although TRIGRS computes the FoS cellwise (i.e., gridded), the hydrogeomechanical parameters are provided zone-wise. The available versions of TRIGRS do not have the ability to provide cell-wise hydrogeomechanical properties. Because the hydro-geomechanical parameters were derived from the soil texture information, the zones in TRIGRS approximately followed the soil texture map of the study area. The present study incorporated FAO-STATSGO soil information and high-resolution regional soil information from DSSSC, which resulted in three and eight soil types in the study area, respectively (Figure 4). Therefore, the study area was categorized into three and eight zones. Although the spatial heterogeneity was better represented with the DSSSC soil texture, the modeled FoS from both soil textures were comparable (Section 3.1), leading to a lack of improvisation. This may have been due to the use of SHPs based on soil texture. This was further clarified by carrying out a sensitivity analysis to quantify the impact of the SHP on the analysis of slope stability and to further understand the absence of significant changes in modeled FoS when the two soil texture maps were used to derive the SHPs.

Because SHPs are important in determining the stability of hillslopes, saturated water content (θ_s); residual water content (θ_r); and a fitting parameter for soil size distribution (α), which is approximately equivalent to the inverse height of capillary rise [54], were selected for sensitivity analysis. In order to understand the impact of each parameter on the final grid-based FoS, one-parameter-at-a-time (OAT) analyses were carried out. OAT analyses are often used in sensitivity analyses, parameter optimization, and calibrating physically based models [42,99,100,101].

The current study considered a case of a finite and unsaturated profile. The finite depth was based on Equation (1). Based on Equations (4) and (5), the FoS was calculated for both the saturated and unsaturated zones, where each cell was treated as a unique sliding unit. Thus, failure analysis with pressure head determination provided the FoS for the depth profile for each cell. It is noted that the depth at which the FoS first reduced below one was considered the landslide initiation depth, which could vary depending on the soil properties as well as the rainfall intensity [92,102]. The depth may vary between the water table depth and the basal boundary [103]. The pressure head, utilized in Equations (4) and (5), was calculated based on the below equations based on Richard’s equation (Equation (7)).

$$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}\left[K\left(\psi \right)\left(\frac{1}{co{s}^2\delta }\frac{\partial \psi }{\partial Z}-1\right)\right]$$

(7)

where,

$$K\left(\psi \right)=K_s\exp \left(\alpha {\psi }^{*}\right)$$

(8)

The pressure head was obtained using Equation (9).

$$\psi \left(Z,t\right)=\frac{cos\delta }{{\alpha }_1}ln\left[\frac{K\left(Z,t\right)}{Ks}\right]+{\psi }_0$$

(9)

A realistic range of θ_s, ranging from 0.280 to 0.450, with a median of 0.365, was chosen to understand the impact of the saturated water content on the stability of terrain and to determine the landslide initiation zone (Figure 15). It can be seen that with the same rainfall intensities, a change in θ_s impacted the FoS of the depth profile (Figure 15) as well as the time taken for the rising of the groundwater table (

Table 2. Comparative analysis for landslide initiation based on variation in θ_s.

Sl. No.	θs	FoS < 1 (Landslide Initiation Day)	Depth to Groundwater Table	Depth of Weak Zone
1	0.28	Day 5	0.50 m	1–2 m
2	0.32	Day 6	0.48 m	1.7–2 m
3	0.36	Day 7	0.46 m	1.8–2 m
4	0.45	Day 9 (Stable)	0.76 m	Stable

). The FoS corresponding to θ_s = 0.365 reached < 1 on the seventh day. At the minimum θ_s, an FoS < 1 was seen on the fifth day, while the cell remained stable for θ_s = 0.450 (Figure 15). Figure 15 further quantifies the change in FoS (%) with respect to changes in θ_s (%), where a 30% increase in θ_s reflects a 15% change in the FoS. The OAT analyses carried out as part of the present study for θ_s showed an approximate linear trend for the study area. Additionally, it can be seen from Table 2 that the landslide was initiated as soon as the groundwater table rose to depth of 0.50 m. The water table rise mechanism was already previously explained (Baum et al., 2010). Thus, as shown in Take et al., 2004, with the rise in pore water pressure to zero, the shear resistance fell. This led to an unstable zone. It is illustrated from both the depth profile in Figure 15 and Table 2 that the zone of weakness increased with the decrease in θ_s, which represents the porosity of the soil. It is deciphered from Figure 15 that θ_s had a significant role in estimating FoS and, thereby, the initiation time and spatial distribution of the landslides in the study area.

As described above, an OAT analysis was carried out for θ_r. The θ_r values were considered from 0.01 to 0.1, with a median of 0.55, and the corresponding FoSs were modeled for a 9-day period and are plotted in Figure 16. The median θ_r (0.055) reached the threshold value before the seventh day, and the soils with lower θ_r values needed more time, up to 36 more hours, to become unstable. This change could be attributed to the increase in the water-holding capacity of the soil. On the other hand, higher θ_r values started crossing the threshold FoS of one starting on the fifth day. It can be seen that a 100% change in θ_r produced only a 4% change in the FoS. It might be because of the fact that the θ_r alone might have a lesser impact on the modeled FoS. In reality, landslide is a very complex phenomenon, and the interplay of multiple factors determine the mechanism behind its initiation. Furthermore, in order to understand the impact of θ_r together with θ_s and other factors, many-parameter-at-a-time (MAT) analyses should be carried out by simultaneously varying multiple parameters in the model. Even though MAT analyses will further shed light on the interconnection among various hydrogeomechanical parameters, the present study was limited to OAT because of its simplicity and the significantly less computational effort required.

Along with θ_s and θ_r, α (a fitting parameter for soil size distribution and is approximately equivalent to the inverse height of capillary rise) was analyzed to comprehend its impact on FoS and thereby the initiation of landslides. The α values were varied from 0 to 10, and we modeled the FoS for a 9-day period. The experiments were carried out as in the cases for θ_s and θ_r and are plotted in Figure 16. The fall in FoS for 9 days with the α values 0.35, 1.00, 1.80, 2.00, and 10.00 are highlighted in the plot. It is observed from Figure 16 that α values above two were less significant in determining the FoS than values from zero to two. Because any α value above two had less significance in determining the FoS in the study area, a range of zero to two with a median of one was considered for analyzing the percentage change in the FoS with respect to the percentage change in α (Figure 16). It has to be noted that any α value less than zero acts as a flag in TRIGRS to treat the terrain as tension-saturated.

Two major observations were made as part the present study while selecting the input parameters for modeling the slope stability in complex terrain. Primarily, even though the regional DSSSC soil map had better soil information representing the study area, it still considered a mean θ_s, θ_r, and α for a zone based on the field observations and laboratory measurements of the soil samples. In reality, the SHPs vary for the same soil type based on the porosity, grain size, degree of compaction, etc., and complex terrain is more likely to reflect this trend. The subzonal variation must have strongly affected the modeled distribution and the computed initiation time of the landslides in the study area. Secondly, because the uncertainties present in spatially varying hydrological properties have an impact on the initiation time and spatial distribution in complex terrain, a robust approach has to be implemented to (1) accurately derive the hydrological properties through EO data sets in finer resolution and (2) provide the derived finer-resolution hydrological properties cell-wise rather than zone-wise as in the present versions of TRIGRS. An obvious solution for deriving hydrological properties that have a significant impact on the FoS is through inverse modeling. Although in situ collected soil samples and laboratory analyses are rather more intuitive than modeling approaches, it has to be noted that in situ data collection, sampling in uniform intervals (spatial and temporal), and analyses are very costly and often impossibly laborious for complex and inaccessible terrain.

Along with the results and findings of the present study, the significance of the inherent uncertainties present in the hydraulic parameters in determining the model outputs was also put forward in a previous study on a landslide-prone area in Brazil [104]. The study further emphasized the need for developing novel strategies to generate hydraulic parameters that represent spatial variation in finer resolution. A recent review [45] on TRIGRS and its performance compared with that of other physical slope stability models also concluded that refined hydrological parameters can generate more realistic results. Another study carried out in Norway focused on the initial conditions of the terrain and concluded that wetter initial conditions lead to early instability and overestimation of landslide pixels in the study area [105]. A data-sparse region in the northern part of Kerala was studied for landslide hazard zonation using TRIGRS and showed that providing spatially varying SHPs significantly reduces the overestimation of landslide pixels [62]. The study further raised concerns regarding the lack of flexibility to provide local variations in hydrogeomechanical parameters to the current versions of TRIGRS. The study suggested the modification of TRIGRS to receive input parameters with local variations rather than generalizing the properties zone-wise. A number of methods to inversely derive SHPs through EO data sets have been put forward to overcome the limitations of regression functions [106,107,108]. The high-resolution SHPs derived through inverse modeling [108] can be a possible alternative, especially in areas where field data are difficult to retrieve.

5. Conclusions

The present study modeled the distribution and timing of rainfall-induced shallow landslides in a catchment in Western Ghats. It was shown how a physically based model can be used for understanding the timing and distribution of landslides with minimal ground-based data sets. This study used a set of satellite-derived terrain parameters and hydrogeomechanical parameters with ground-measured rainfall data to compute the FoS and PWP in a landslide-prone area. The model was able to predict the landslide distribution in the spatial region with a TPR of 68%. Additionally, for nine landslides, the model predicted the occurrence of eight landslides, with precise timing of three landslides. In addition, the present study analyzed the impact of the SHPs in computing the FoS and thereby predicting the timing and distribution of landslides in the study area. It could be seen that with the increase in the storage capacity of soil, it became more stable, and a prominent shift was seen in the initiation of the landslide events. Thus, the precise timing of a landslide event was dependent on the SHPs. Although the model could compute the change in the FoS and PWP in accordance with the rainfall received in the study area, it is worth mentioning the challenges and room of improvement when using the model, TRIGRS. Primarily, the lack of high-resolution SHPs was one of the major challenges faced when executing the slope stability model. Although a field-based high-resolution soil texture map from DSSSC was used to generate high-resolution FoS and could reduce the number of false positives, it needs further improvement, especially in computing the timing of landslide initiation. A high-resolution map of SHPs can provide more control and heterogeneity and thereby improve the modeling of the FoS and PWP. Because the hydrologic response of terrain is determined by the spatial distribution of the thickness of soil columns, any physically based slope stability model heavily depends on the resolution and accuracy of the soil depth. In TRIGRS, the depth to slip surface and time taken to reach complete saturation of each cell are computed as a function of soil depth. It further emphasizes the requirement for the high-resolution spatial distribution of soil depth. Although it is extremely difficult to obtain in situ observations of soil thickness for larger and complicated terrain, a higher number of in situ and uniformly distributed soil thickness sampling observations would resolve this problem to a certain extent. Although necessary steps, such as pit removal, were taken prior to the analyses, the present study heavily relied upon SRTM 30 m DEM to derive the terrain and hydrological parameters, such as slope, flow direction, and TopoIndex parameters used in TRIGRS. A high-resolution DEM from airborne surveys, such as drone-based surveys, could have significantly improved TRIGRS’ outputs.

Furthermore, the dynamic physically based slope stability model, TRIGRS (version 2.1), does not provide the flexibility to input the hydrogeomechanical parameters grid- or pixel-wise. The grid-wise provision of input data sets with the inclusion of the impact of local geological and anthropogenic features and hydrogeomechanical properties can significantly improve the model outputs. Although the grid-wise provision of the hydrogeomechanical properties will be computationally expensive, it may significantly enhances modeling capabilities. On the basis of the results obtained in the present study and the sensitivity analysis carried out, we further plan to generate SHPs in finer resolution and provide SHPs cell-wise, rather than grid-wise, in TRIGRS.

Despite of the limitations of this study and the unavailability of high-resolution terrain and hydrogeomechanical information, the model could fairly demarcate the distribution and capture the timing of rainfall-induced shallow landslides in the study area. It can be further concluded that the model and method can be used as a measure to assess the landslide vulnerabilities in the Western Ghats area or anywhere else where the input data sets are available. Moreover, the physically based slope stability model used in the present study can be used as a primary approach to comprehensively understand the stability of hillslopes where input data sets are limited or only remotely sensed EO data sets are readily available.