1. Introduction

Human and economic losses generated by landslides occur every year in all countries; however, the impact of landslides varies considerably according to the local geological conditions and socio-economic vulnerability ¹⁾⁽². Although landslides do represent changes in terrain morphology within the natural and continuous geomorphological cycle ³ their occurrence in recent decades has been closely tied to world population growth and consequent urban expansion on susceptible slopes to this type of processes ⁴⁾⁽⁵.

Landslides are caused by a variety of hydrological, geological, and anthropological factors, compelling researchers to take an interdisciplinary approach to predicting their occurrence ⁶. The factors controlling the occurrence and distribution of landslides can be categorised as quasi-static variables or dynamic variables. The quasi-static variables contribute to slope susceptibility, which defines the spatial distribution of landslides, and dynamic variables control the triggering factors on landslide-prone-slopes and characterise landslides’ temporal patterns, such as rainfall ⁶⁾⁽⁷.

Regarding to landslides induced by rainfall, two different failure mechanisms have been discussed in the literature. The first mechanism is based on the idea of the rainfall rate exceeds the percolation rate, creating a perched subsurface water flow in the residual soil parallel to the slope. Failure occurs due to increased pore pressure caused by the rapid reduction of shear strength in undrained conditions, which has been called static liquefaction of the material by several authors ^{8)(9)(10)(11)(12}. The second failure mechanism proposes the development of an advancing wetting front from the slope surface in which the material is still in an unsaturated state when failure occurs due to reduced suction, as the mass behaves like a rigid body ^13)(14)(15.

The complexity of calculating the probability of obtaining an increased critical pore pressure or reduced suction and, consequently, predicting a rainfall-triggered landslide occurrence, is a function of many intimately related parameters. One of the most important recent advances that has allowed researchers to consider all of these variations is the use of physically based models ^{16)(17)(18)(19)(20)(21)(22}.

Physically based modelling of rainfall-triggered landslides proposes hillside hydrology as a subsurface flow in a steady state ²³ or vertical transient dynamic flows ²⁴. Other models based on the static state in kinematic wave hydrology have been used for saturated hillside slopes in many studies ^25)(26)(27. Several models are available in the literature. One of the first and most recognised physical models is called SHALSTAB ²³. SINMAP (Stability Index MAPping) ²⁸, and LISA (Level I Stability Analysis) ²⁹, are similar models using the hydrological SHALSTAB approach. The dSLAM (distributed Shallow Landslide Model) is a coupled hydrological-geotechnical model modified as the Integrated Dynamic Slope Stability Model (IDSSM) ¹⁶⁾⁽³⁰. CHASM (Combined Hydrological and Stability Model) is a distributed physically based two-dimensional model ³¹. TRIGRS is a Fortran program ³², based on the transient one-dimensional vertical infiltration model created by Iverson ²⁴. A similar approach, using the analytical approximate solution of the Richards equation, is called HIRESSS (High REsolution Slope Stability Simulator) ³³. Recently, a model called GEOtop-FS was proposed ³⁴; and tRIBS (TIN-based Real-Time Integrated Basin Simulator), an infinite slope and movement module in a physically based, distributed hydrological model ³⁵⁾⁽³⁶.

All these available models have been applied to very specific environmental conditions that do not allow for adjustment to particular rainfall and terrain complexities, such as those for tropical and mountainous terrains. A conceptual and physically based model called SHIA_Landslide for the prediction of rainfall-triggered shallow landslides in tropical environments and complex terrains was developed by ³⁷. The model is based on the Open and Distributed Hydrological Simulation methodology ³⁸. SHIA_Landslide is an original and significant contribution that offers a new perspective to analyse shallow landslide processes by incorporating a comprehensive distributed hydrological tank model that includes water storage in the soil and a geotechnical, classical analysis of infinite-slope stability under saturated conditions. The model was proposed to improve understanding of the mechanisms associated with slope instability and rainfall infiltration in mountainous areas located in rainy environments, where increased population pressure has led to the expansion of development into landslide-prone areas.

More recently, a number of studies has examined the forecasting performance of linear and non-linear hydrological models ³⁸⁾⁽³⁹. Linear models have been widely applied in hydrology simulation, and especially in models based on unit hydrograph. However, several authors do not consider it a safe practice because hydrological events of great magnitude underestimate the discharge flow when nonlinearities are not considered. Nonlinearity is more critical for small catchment in complex terrains during intense rainstorms ⁴⁰.

This paper compares the forecasting performance of linear and nonlinear SHIA_Landslide model, to evaluate the nonlinearity effects on the performance of a combined hydrology and slope stability model in tropical and mountainous terrains. The study simulates the September 21, 1990 rainstorm event in the La Arenosa catchment, where more than 600 of landslides were triggered during just three hours of rainfall.

2. The model SHIA_Landslide

SHIA_Landslide is a FORTRAN program for computing positive pore pressure changes as well as resulting changes in the factor of safety due to rainfall infiltration, using a physical and conceptual based, distributed hydrological and geotechnical coupled model to provide an assessment of slope-failure condition.

The hydrological module of SHIA_Landslide is formed by two fundamental components: a water balance that simulates the dominant hydrological processes in the catchment and a routing component that simulates the flow of water through the river network. Each grid cell corresponds to a system of five interconnected tanks that communicate with their respective tank in the downstream cell. The first four tanks represent the basin’s runoff production processes, while the last tank represents the transfer process runoff. A more detailed description of the SHIA_Landslide can be found in ³⁷.

The first tank (T1) is called static storage and represents interception and water detention in puddles and the capillary water storage in the soil-rooting zone, which is a function of field capacity and effective root depth. The only outflow from this storage is real evapotranspiration ( E 1 ). The second tank (T2) is called surface storage and represents water on the hill’s sloped surface that is flowing over the slope and has not infiltrated. The third tank (T3) represents the gravitational water storage in the residual soil between field capacity and saturation. This tank models the water column due to subsurface flow parallel to the slope surface through the soil layer and into the drainage system. The fourth tank (T4) corresponds to the aquifer, where vertical flow represents the system’s groundwater outflow and horizontal flow represents the base flow. The final tank (T5) represents the stream flow channel at the cell level, where each cell is connected to the downstream cell according to the drainage network.

The vertical connections between tanks describe the rainfall, evapotranspiration, infiltration and percolation processes. Based on the tropical weathering profile, each grid cell is formed by three soil horizons with three different saturated hydraulic conductivities: residual soil (𝐾𝑠), saprolite (𝐾𝑝), and rock (𝐾𝑝𝑝), where the residual soil is more permeable than the underlying saprolite soil and the impermeable rock.

The horizontal connections describe the overland flow, subsurface flow and base flow. The model considers the horizontal transfer of water to eight possible neighbouring cells by using a single flow direction sub-model. Horizontal connections are based on the topology of the basin built on a grid cell type.

To evaluate slope stability after the hydrological hillslope simulation, SHIA_Landslide uses an infinite slope stability analysis. Hillslopes, with steady subsurface flow parallel to the slope and the perched water level at distance 𝑍 𝑤 above the slide surface, have a pore water pressure of 𝜇= 𝛾 𝑤 𝑍 𝑤 𝐶𝑜𝑠 2 𝛽, ⁴¹, and therefore the factor of safety (FS) is shown in Eq. ¹:

(1)

where γ is the unit weight of soil, γ_w is the unit weight of water, Z is the soil thickness measured vertically, 𝛽 is the gradient of the hillslope, and 𝜙 is the friction angle.

For the calibration process, SHIA_Landslide includes the split-parameter procedure for distributing model proposed by ⁴². The effective parameters for the model at each cell are split into two components, a hydrological or geotechnical characteristic and a correction factor common to all cells that accounts for all modelling errors. Correction factors account for time and space scale effects together with the model and input errors, leaving the hydrological or geotechnical characteristics free of these problems while maintaining the parameter’s physical meaning.

To evaluate the model performance, and if observed stream flow data is available, SHIA_Landslide provides the root mean square error (RMSE) and Nash-Sutcliffe efficiency coefficient (NS). Moreover, the model uses the total rainfall over the catchment and the outflows along the different grid cells for water balance, and considers the water storage in the tanks at the end of the simulation.

The model outputs four products during the simulation process. The hydrological module gives the hydrograph and perched water table variation for any point in the catchment. In addition, the geotechnical module outputs two maps, a susceptibility map before the target area is exposed to rainfall, and a landslide hazard map showing the spatial distribution of landslides triggered by rainfall.

Figure 1 shows the flow diagram structure of the coupled model, SHIA_Landslide, delineating the steps used in the hydrological and slope stability calculations.

Figure 1 Flow chart of SHIA_Landslide program 

3. Linear and nonlinear approach

A great advantage of the Open and Distributed Hydrological Simulation methodology implemented in SHIA_Landslide is that permits to include different approaches to estimate the flows in the vertical and horizontal direction. SHIA_Landslide permits to estimate the overland flow and subsurface flow by a linear reservoir equation or a non-linear approach.

In the linear reservoir approach, in which flow velocity ( 𝐸 x ) is assumed constant and independent of flow magnitude, the discharge correlates water storage and a discharge coefficient, 𝛼, according to Eq. ²:

(2)

where the water level of the tank 𝑥 is represented by 𝑆 x and the discharge coefficient of the linear reservoir, 𝛼, is a function of the cell size 𝑑𝑥, the temporal discretization 𝑑𝑡 and a constant flow velocity, 𝑣 x . For the overland flow, the velocity is a constant value related to the hillslope velocity, several values have been proposed in the literature. For the subsurface flow, the constant velocity is given by the horizontal hydraulic conductivity of the upper part of the soil, this being mainly defined by its macro pore structure.

To estimate the overland flow and subsurface flow by a nonlinear approach, the flow velocity is a function of the water volume, specifically the water level of the tank and the flow section.

For the overland flow by a nonlinear approach, most authors have recommended a uniform flow such as the Manning equation, where the slope of the energy line is similar to the slope angle of the terrain ³⁸. In this way, the equation for overland runoff velocity ( 𝑉 2 ) is a function of the transversal flow section (A), slope angle (𝛽) and Manning coefficient (n) according to Eq. ³:

(3)

where ξ and e₁ are parameters associated with the surface type. For flows over natural terrains values of 0.038 and 0.315, respectively, were recommended by ⁴³.

For the subsurface flow by a nonlinear approach, the flow velocity ( 𝑉 3 ) is estimated as a lateral subsurface flow in mountainous terrains covered by forests according to Eq. ^{4 44}:

(4)

where K_s is the saturated hydraulic conductivity of the residual soil, and b is dependent on the soil type. A value of b = 2 was recommende by ⁴⁴ for a mountain basin covered by forest that represents a non-homogeneous hydraulic conductivity along the weathering profile.

For the base flow, SHIA_Landslide only considers a linear reservoir approach, where the constant velocity is related to the aquifer’s saturated hydraulic conductivity.

Finally, the routing along the channel network is carried out at a non-stationary velocity using the Geomorphological Kinematic Wave (GKW) ³⁸. The GKW is a simplification of the Saint Venant equations, in which inertial and pressure terms are neglected. The GKW uses the correlation proposed by ⁴⁵, which relates the cross-sectional geometry and velocity to the river’s discharge using potential equations, the Manning’s equation, which relates the flow velocity and flow cross section, and an empirical roughness coefficient. Using these correlations, the Eq. ⁵ is obtained for the velocity of the water ( 𝑉 5 ) in the channel as a function of the channel geometry and the terrain’s geomorphology:

(5)

The GKW requires nine independent exponents and coefficients, which can be obtained from a geomorphological regional study of hydrologically homogeneous zones. However, empirical studies have been performed by multiple authors that propose different values according to local conditions. Table 1 shows the constant regional parameter range values proposed by different authors ^{38)(39)(40)(41)(42}.

Table 1 Geomorphological Kinematic Wave parameter ranges proposed for the model 

4. Case study

A short-duration, high-intensity rainfall event affected the basin of La Arenosa on 21 September 1990. In less than 3 hours, a precipitation of 208 mm fell within the study area, triggering many landslides. The September 21, 1990 event is unique considering the huge number of failures that took place as a result

La Arenosa, with an area of 9.91 km², is located 160 km east of the Aburrá Valley, on the south-eastern side of the Central Cordillera in the Antioquia region (Figure 2). The area has a tropical humid climate with a mean annual precipitation of 4,300 mm and a mean annual temperature of 23° ⁴⁶. The geology of the study area consists of residual soils from granodiorite rocks covered in the gently sloping areas with slopes and fluviotorrential deposits ⁴⁷.

Figure 2 Location of the La Arenosa catchment, in the southeastern side of the humid tropical and complex terrains of Central Andean Cordillera in Colombia 

Analysis of post-event aerial photos and field investigations allowed for a partial reconstruction of the pattern and characteristics of the landslides in the La Arenosa catchment. A detailed landslide inventory and a comprehensive description of the landslides triggered during the event was carried out by ⁴⁷⁾⁽⁴⁸, and reported 699 landslides in the La Arenosa catchment. These landslide inventory maps cover only about 70% of the catchment.

SHIA_Landslide considers three morphometric input parameter maps, namely accumulated area, flow direction, and slope angle; and nine soil input parameter maps, namely soil thickness, cohesion, saturated unit soil weight, friction angle, saturated hydraulic conductivity of residual soil, saturated hydraulic conductivity of saprolitic soil, saturated hydraulic conductivity of rock, maximum static water storage, and maximum gravitational water storage

Morphometric parameters such as slope angle, flow direction and flow accumulation were calculated using ArcGIS 10.1 hydrological tools and a DTM provided by the Instituto Geográfico Agustin Codazzi (IGAC) with a raster size of 10 m. Soil parameters were obtained according to the official soil map elaborated by IGAC ⁴⁶, and soil descriptions, field tests and laboratory analyses on soil samples ⁴⁷⁾⁽⁴⁸ (Table 2).

Table 2 Morphometric and soil parameters of La Arenosa 

5. Calibration and validation

Model calibration for the hydrological component was conducted for the period from March to May, 2011. Figure 3 show the rainfall data for these periods.

Figure 3 Rainfall time series of La Arenosa rain gauge for the period March and May 2007 

The calibration results for La Arenosa discharge for nonlinear and linear version are shown in Figure 4, respectively. It illustrates the simulated discharge provided for linear and nonlinear SHIA_Landslide compared to the stream flow of La Arenosa for the calibration period. Overall, the simulated hourly flows resulting from the nonlinear and linear SHIA_Landslide agreed very well with La Arenosa stream flow. Both model approaches predicted flow peaks and slightly under-predicted low flows.

Figure 4 Results using nonlinear and linear SHIA_Landslide of simulated hourly discharges at the calibration flow gauge station La Arenosa during the calibration period compared with the discharges obtained for La Arenosa stream 

Comparison of predicted and discharge data of La Arenosa suggest that the nonlinear model gives good prediction of the landslide hydrology. Most of the peaks are simulated for the model with a good time precision. For the highest peaks the model slightly overestimates the stream flow. The root mean square error (RMSE = 0.292) and the Nash - Sutcliffe coefficient (NS=0.852) show good correlation between simulated and La Arenosa stream flows. The water balance shows a low negative difference (WB= -1.59 %).

For the linear version, the simulation is acceptable with values of RMSE (0.347) and NS (0.791), slightly less than for the nonlinear version. Moreover, the water balance shows similar negative values (WB = -1.61%). Although the linear model simulates adequately the event, it simulates recession curves substantially lower when compared to La Arenosa stream flows.

Figure 5 shows the simulated perched water table level using nonlinear and linear model for the calibration period evaluated in a slope grid cell with an accumulated drainage area of 800 m², it means 8 slope grid cells flow through this grid cell. During all the simulated period, perched water table conserve a level around 0.2 m, and the maximum value was 0.6 - 0.7 meters, correlated to a rainfall peak in the last week of April.

Figure 5 Nonlinear and linear simulated perched water table level for an slope grid cell (accumulated area = 800 m²) of La Arenosa for the period between March and May 2011 

For the linear model, the perched water table values decrease, the two remarkable peaks are lower, and the valleys get values of zero for short time periods. Evaluating the slope of the peaks, it can be seen that the linear version obtains a faster positive pore pressure response of the soils, increasing and reducing the peaks quicker than the nonlinear version.

Considering that validation provides a direct measure of the degree of uncertainty that may be expected when the model is applied to conditions outside of the calibration series, several verification tests were applied to evaluate the performance of the model. One representative validation test is presented, from September 2012 to November 2012.

Figures 6 and 7 show rainfall time series and the discharge simulated compared to the La Arenosa discharge.

Figure 6 Rainfall time series for La Arenosa rain gauges 

Figure 7 Simulated discharge from nonlinear and linear SHIA_Landslide compared to La Arenosa discharge for the period September to November 2012 

Overall, the simulated hourly flows agreed well with La Arenosa stream flows for both nonlinear and linear SHIA_Landslide. Despite occasional discrepancies for peak flows, most peaks are predicted correctly.

For the nonlinear version, the RMSE is 0.297 and NS efficiency 0.724. The water balance decreases compared to calibration value (-0.568). For the linear version, values decrease conserving good correlation values, RMSE is 0.322 and NS is 0.677. The water balance was -0.62%.

For the perched water table simulated for SHIA_Landslide, results are consistent with the validation (Figure 8). The small peaks show lower values for the nonlinear version, but the highest peak has similar maximum values for both versions. Similar to the calibration procedure, the nonlinear version conserve higher minimum values than the linear version.

Figure 8 Nonlinear and linear Simulated perched water table level slope grid cell (Accumulated area = 800 m²) 

6. September 21, 1990 rainstorm simulation

Figure 9 shows the rainfall time series for the September 21, 1990 rainstorm event. A period of three months was selected to assure that the initial conditions previous to the event are reached previously to the target time.

Figure 10 shows the simulated discharge for the September 1990 rainstorm for the nonlinear and linear version of the model. According to SHIA_Landslide, the channel discharge flow at the point of La Arenosa station was 55 m³/s for the nonlinear version, and a flow discharge of about 35 m³/s for the linear version. It means more than 30 times the base flow of La Arenosa. It is important to note that this discharge only considered the water flow; according to the description provided by ⁴⁹, the event was characterized by a huge volume of sediments, estimated by the authors in 1.5 Mm³. Considering this volume, the total water and sediments discharge increased considerably.

Figure 9 Rainfall from La Arenosa rain gauge for the period July to September 1990 

Figure 10 Stream flow simulated for nonlinear and linear SHIA_Landslide for the period July to September 1990. Observe discharge is not available to compare with simulated data 

Figure 11 shows the nonlinear and linear simulated perched water table level for the September 1990 event evaluated in the same slope grid cell used for the hydrological calibration, which has an accumulated drainage area of 900 m². Perched water tables reached a peak of 1.3 m during the most intense period of the rainstorm. The peak obtained for the linear version is higher than the peak reach for the nonlinear version. For the nonlinear version, there is a minimum level that is conserved for the perched water table close to 0.1m, even for periods with very low rainfall,

Figure 11 Perched water table level simulated for nonlinear and linear SHIA_Landslide for the period July to September 1990 for a slope grid cell (800 m²) 

Figure 12 shows the susceptibility map for La Arenosa catchment. Landslide susceptibility map does not have any change due to it only depends on the geotechnical parameters, such as slope, cohesion, friction angle, soil unit weight, which are assumed that do not change during the simulated rainstorm.

Yellow areas are the potential unstable areas for landslides triggered by rainfall, which correspond to 52% of the total area. 92% of the scars of the September 1990 event are located into these areas, the rest 8% corresponds to an induce error in the hazard map (Table 3), due to there is no way to predict this small percentage of landslides in the next step of the model, because these areas have been identified as unconditionally stable by the model.

Figures 13 and 14 shows the hazard map provided for the nonlinear and linear version of SHIA_Landslide respectively, providing the areas with landslide occurrence triggered by rainfall.

The areas prone to rainfall-induced landslides change considerably depending on the model used. For the nonlinear version (Figure 13) 23,456 grid cells that fail during the rainstorm were identified, which amount to 24% of the total catchment area. For the linear version (Figure 14) the model identified 34,268, which correspond to 35% of the total catchment area (Table 4). Comparing with the nonlinear version, it means an overestimation of 11%.

Figure 12 Susceptibility map obtain for SHIA_Landslide 

Table 3 Susceptibility grid cells classification by SHIA_Landslide 

Figure 13 Areas with landslide occurrence triggered by the September 1990 rainstorm simulated by nonlinear SHIA_Landslide 

Table 4 Hazard map: landslides triggered by rainfall for linear and nonlinear SHIA_Landslide for the period July to September 1990. Partial percentage is estimated according to the total number of potentially unstable grid cells (51,176), and total percentage is estimated according to the total number of grid cells (98,646). 

7. Results and discussion

A comparative analysis between linear and nonlinear SHIA_Landslide was performed to examine which model would provide better predictions. The results from the tests are summarized in Table 5.

In terms of the hydrological component, the nonlinear model showed the best performance, with the highest NS efficiency and lowest RMSE values. Concerning the geotechnical component, both approaches identified the same number of observed scar grid cells; however, the linear model increases the number of unstable grid cells. The difference resides on the total number of unstable grid cells identified. The nonlinear SHIA_Landslide model identified a number of 23,456 grid cells as unstable, 46% of potentially unstable cells, and the linear SHIA_Landslide identified a total number of 34,268 grid cells as unstable, 67% of potentially unstable cells. That means an increase of 21.13% of potentially grid cells identified as unstable. Furthermore, considering that both approaches identified similarly the actual scars grid cells, it means this additional percentage could be erroneous grid cells identified as unstable for the linear model, or unstable grid cells not included in the partial reconstruction of the landslide inventory maps.

Table 5 Comparison between nonlinear and linear SHIA_Landslide. Maximum rainfall intensity (MRI). 

One of the disadvantages of nonlinear model is that the parameter-fit process is iterative, consuming valuable time for an early warning system. To estimate the velocities of the overland runoff and subsurface flow the nonlinear process must start with a set or user-supplied starting values. SHIA_Landslide then tries to find the correct velocities by adjusting the values provided using velocities equations that are function of the transversal flow section. All this process takes additional time.

For the La Arenosa catchment, composed by 98,646 grid cells with a simulation period of three months, which means a number around 2,200 hours, the nonlinear model takes 5 minutes, whereas the linear model takes about 1 minute for the same conditions. The nonlinear model obviously requires more time to execute, due to the area-velocity iterative process to arrive at the velocity fitting parameter. The linear model is approximately 5 times faster than the nonlinear model. Although time is a valuable factor for early warning system, the time used for the nonlinear model is reasonable for an early warning system. The results suggest that non-linearity increases the predictions of the model, conserving appropriate time running.

Finally, it is well known in hydrology simulation that linear models overestimate high frequency floods and underestimate low frequency floods, provided both models have been adequately adopted for mean frequency floods ⁵⁰. Hydrologic systems are nonlinear and the implications of this non-linearity should be taken into account in the formulation and application of distributed models ⁵¹.

The results show than non-linearity on hillslope hydrology simulations affect the forecasting of landslide triggered by rainfall, then must be consider into the analysis. Non-linearity increases the predictions of the model, conserving appropriate time running; and linear models overestimate the potentially unstable areas.

8. Conclusions

SHIA_Landslide is a modelling program for computing positive pore pressure changes and associated changes in factor of safety due to rainfall infiltrations using a conceptual and distributed hydrological module coupled with a physically based and infinite stability slope geotechnical module. SHIA_Landslide represents an innovative approach that introduce new concepts for the simulation of shallow landslides triggered by rainfall. It uses a robust and simple conceptual hydrological model to simulate all the main components of the land phase of the hydrological cycle, and to establish the different components of hydrographs at multiple sites. These advantages of the model allow a better understanding of the hydrology of slopes, and consequently the stability of slopes.

SHIA_Landslide could be implemented using linear and nonlinear flow velocities for the different tanks levels. Comparisons between nonlinear and linear model indicate that, although a linear approach of the model takes much less time to run, a nonlinear SHIA_Landslide shows much better performance, decreasing false alarm rate. Considering this conclusion, SHIA_Landslide is an interesting tool that can be implemented in the development of early warning systems, aid in providing real-time rainfall monitoring and assist in the dissemination of alerts and communications.

The results show than consider the non-linearity on hillslope hydrology increases the predictions of physically based modelling of rainfall-triggered landslides; therefore, more effort should be put in the development of nonlinear hydrologic models.

9. Acknowledgements

The authors wish to thanks the Hans Wildsdorf Foundation and the Colombian Association of Petroleum Geologists and Geophysicist (ACGGP) for providing financial support for conducting this research