Articles | Volume 19, issue 4
Research article
09 Apr 2019
Research article |  | 09 Apr 2019

Development and validation of the terrain stability model for assessing landslide instability during heavy rain infiltration

Alfonso Gutiérrez-Martín, Miguel Ángel Herrada, José Ignacio Yenes, and Ricardo Castedo

Slope stability is a key topic, not only for engineers but also for politicians, due to the considerable monetary and human losses that landslides can cause every year. In fact, it is estimated that landslides have caused thousands of deaths and economic losses amounting to tens of billions of euros per year around the world. The geological stability of slopes is affected by several factors, such as climate, earthquakes, lithology and rock structures, among others. Climate is one of the main factors, especially when large amounts of rainwater are absorbed in short periods of time. Taking this issue into account, we developed an innovative analytical model using the limit equilibrium method supported by a geographic information system (GIS). This model is especially useful for predicting the risk of landslides in scenarios of heavy unpredictable rainfall. The model, hereafter named terrain stability (or TS) is a 2-D model, is programed in MATLAB and includes a steady-state hydrological term. Many variables measured in the field – topography, precipitation and type of soil – can be added, changed or updated using simple input parameters. To validate the model, we applied it to a real example – that of a landslide which resulted in human and material losses (collapse of a building) at Hundidero, La Viñuela (Málaga), Spain, in February 2010.

1 Introduction

Landslides, one of the natural disasters, have resulted in significant injury and loss to human life and damaged property and infrastructure throughout the world (Varnes, 1996; Parise and Jibson, 2000; Dai et al., 2002; Guha-Sapir et al., 2004; Crozier and Glade, 2005; Khan, 2005; Toya and Skidmore, 2007; Raghuvanshi et al., 2015; Girma et al., 2015). Normally, heavy rainfall, high relative relief and complex fragile geology with increased manmade activities have resulted in increased landslide (Gutiérrez-Martín, 2015). It is essential to identify, evaluate and delineate landslide hazard-prone areas for proper strategic planning and mitigation (Bisson et al., 2014). Therefore, to delineate landslide susceptible slopes over large areas, landslide hazard zonation (LHZ) techniques can be employed (Anbalagan, 1992; Guzzetti et al., 1999; Casagli et al., 2004; Fall et al., 2006).

Landslides are a result of intrinsic and external triggering factors. The intrinsic factors are mainly geological factors or the geometry of the slope (Hoek and Bray, 1981; Ayalew et al., 2004; Wang and Niu, 2009).

The external factor which generally triggers landslides is rainfall (Anderson and Howes, 1985; Collison et al., 2000; Dai and Lee, 2001). Several LHZ techniques have been developed in the past, and these can be broadly classified into three categories: expert evaluation, statistical methods and deterministic approaches (Wu and Sidle, 1995; Leroi, 1997; Guzzetti et al., 1999; Iverson, 2000; Crosta and Frattini, 2003; Casagli et al., 2004; Fall et al., 2006; Lu and Godt, 2008; Rossi et al., 2013; Raia et al., 2014; Canili et al., 2018; Zhang et al., 2018). Within these categories, we want to highlight the empirical models that are based on rainfall thresholds (Wilson and Jayko, 1997; Aleotti, 2004; Guzzetti et al., 2007; Martelloni et al., 2011). Each of these LHZ techniques has its own advantage and disadvantage owing to certain uncertainties on account of factors considered or methods by which factor data are derived (Carrara et al., 1995). Limit equilibrium types of analyses for assessing the stability of earth slopes have been in use in geotechnical engineering for many decades. The idea of discretizing a potential sliding mass into vertical slices was introduced in the 20th century. During the following few decades, Fellenius (1936) introduced the ordinary method of slices (Fellenius, 1936). In the mid-1950s Janbu and Bishop developed advances in the method (Janbu, 1954; Bishop, 1955). The advent of electronic computers in the 1960s made it possible to more readily handle the iterative procedures inherent in the method, which led to mathematically more rigorous formulations such as those developed by Morgenstern and Price and by Spencer (Morgenstern and Price, 1965; Spencer, 1967).

Until the 1980s, most stability analyses were performed by graphical methods or by using manual calculators. Nowadays, the quickest and most detailed analyses can be performed using any ordinary computer (Wilkinson et al., 2002). There are other types of software based on the modelling of the probability of occurrence of shallow landslides LHZ, in more extensive areas using geographic information system (GIS) technology and DEM (digital elevation model), as is the case of deterministic models like the software TRIGRS, SINMAP, R-SHALSTAB, GEOtop or GEOtop-FS, and r.slope.stability, among others (Montgomery and Dietrich, 1998; Pack et al., 2001; Rigon et al., 2006; Simoni et al., 2008; Baum et al., 2008; Mergili et al., 2014a, b; Michel et al., 2014; Reid et al., 2015; Alvioli and Baum, 2016; Tran et al., 2018). These are widely used models for calculating the time and location of the occurrence of shallow landslides caused by rainfall at the territorial level, some even in three dimensions, in order to obtain a probabilistic interpretation of the factor of safety. Currently other approaches and/or theoretical studies for landslide prediction are used (for triggering and/or propagation; Martelloni and Bagnoli, 2014; Martelloni et al., 2017). One of the achievements of the presented study is to discretize the potential slip mass in the critical profile of the slope, once unstable areas have been detected through the LHZ (landslide hazard zonation) programs. The terrain stability (TS) calculation tool is not limited to shallow landslides and debris flows but allows analysis of deep and rotational landslides, which other models often do not account for. We use the hydrological variable ru of Spencer in our algorithm to consider the infiltration of rainfall in the calculation of stability of the considered slope.

Limit equilibrium types of analyses for assessing the stability of earth slopes have been in use in geotechnical engineering for last years. Currently, the vast majority of stability analyses using this method of the equilibrium limit are performed with commercial software packages like SLIDE V5, SLOPE/W, Phase2, GEOSLOPE, GALENA, GSTABL7, GEO5 and GeoStudio, among others (González de Vallejo et al., 2002; Acharya et al., 2016a, b; Johari and Mousavi, 2018). Currently there are other slope stability models based on the theory of limit equilibrium that are still in analysis and testing, as is the case with the SSAP software package (Borselli, 2012), but in this case a general equilibrium method model is applied. Secondly, sometimes for commercial models, the introductions of parameters to perform calculations are not very interactive. For the stability analysis, different approaches can be used, such as the limit equilibrium methods (Cheng et al., 2007; Liu et al., 2015), the finite element method (Griffiths and Marquez, 2007; Tschuchnigg et al., 2015; Griffiths, 2015) and the dynamic method (Jia et al., 2008), among others. Limit equilibrium methods are well known, and their use is simple and quick. These methods allow us to analyse almost all types of landslides, such us translational, rotational, topple, creep and fall, among others (Zhou and Cheng, 2013). For the stability analysis, different approaches can be used, such as the limit equilibrium methods (Zhu et al., 2005; Cheng et al., 2007; Verruijt, 2010; Liu et al., 2015), the finite element method (Griffiths and Marquez, 2007; Tschuchnigg et al., 2015; Griffiths, 2015) and the dynamic method (Jia et al., 2008), among others (SSAP 2012, Slide V5, 2018). Also, limit equilibrium methods can be combined with probabilistic techniques (Stead et al., 2006) or with other models, like stability analysis of coastal erosion (Castedo et al., 2012). However, they are limited in general to 2-D planes and easy geometries. Numerical methods – finite element methods – give us the most detailed approach for analysing the stability conditions for the majority of evaluation cases, including complex geometries and 3-D cases. Nevertheless, they present some problems, such as their complexity, data introduction, the mesh size effect, and the time and resources they require (Ramos Vásquez, 2017).

The above-mentioned software packages provide useful tools for determining the stability through the Fs (safety of factor) and for giving the most probable breakage (shearing) surfaces. This technique is fast and allows the field or emergency engineer to make timely decisions. Although this methodology is only available in some current software (Slide V 5.0, STB 2010, GEOSLOPE) and is based on limit equilibrium methods, it is highly recommended because of its reliability for representing real conditions in the field (Chugh and Smart, 1981). This rain infiltration produces a substantial reduction of cohesion (a key soil parameter for stability) that cannot be reproduced by actual software, and then several real situations cannot be predicted.

Delft University of Technology has developed a well-known and free software program to analyse landslides, the STB 2010 (Verruijt, 2010). This program is based on a limit equilibrium technique, using a modified version of Bishop's method to calculate the Fs only for circular failures. It is a user-friendly tool, but it does not allow the calculation of water infiltration on a hillside. This is a critical point, as it is well known that rainfall infiltration is one of the main causes of landslides worldwide (Michel et al., 2015). Reviewing these issues, a new solution must be developed for cases where landslides are linked to heavy rainfall. In this study, we developed a new model and programed it using MATLAB. The primary result of this model was a stability index, namely the minimum Fs, based on the limit equilibrium technique, in this case Bishop's method. The model also provides a possible failure curve and surface area, including the infiltration effects, which can be used to coincide with analysis of the actual event as tested with field data. Topographical data can also be introduced into the model from the DEM in a GIS.

2 Terrain stability model development

In the model we developed, the TS model, we used the limit equilibrium technique for its versatility, calculation speed and accuracy. An analysis can be done studying the whole length of the breakage (shearing) zone or just small slices. Starting with the original method of slides developed by Fellenius (1936), some methods are more accurate and complex (Spencer, 1967; Morgenstern and Price, 1965) than others (Bishop, 1955; Janbú, 1954). Using Spencer's method (Spencer, 1967; Duncan and Wright, 1980; Sharma and Moudud, 1992) here would mean dividing our slope into small slices that must be computed together. This method is divided into two equations, one related to the balance of forces and the other to momentum. Spencer's method imposes equilibrium not only for the forces but also for the momentum on the surface of the rupture. If the forces for the entire soil mass are in equilibrium, the sum of the forces between each slice must also be equal to zero. Therefore, the sum of the horizontal forces between slices must be zero as well as the sum of the vertical ones (Eqs. 1 and 2):


In this equation, Q is the resultant of the pair of forces between slices, and θ is the angle of the resultant (Fig. 1). From this, it can be stated that the sum of the moments of the forces between slices around the critical rotation centre is zero, conformed to Eq. (3):

(3) [ Q R cos ( α - θ ) = 0 ] .

When the R is the radius of the curvature, α is the angle of the slope referred to each slice. This takes into account that the sliding surface is considered circular, so the radius of the curvature is constant.

Figure 1Representation of the forces acting on a slice, considered in Spencer's method (Spencer, 1967). W is the external vertical loads; Zn and Zn + 1 are the forces acting on the left- and right-hand side of each slice, respectively, with their horizontal and vertical components. P and S are the normal and tangential forces at the base of the slice; α is the angle of the slope referred to each slice, b is the slice width and h is the mean height of slice (if the height is not constant).


These equations must be solved to get the Fs and tilt angles of the forces among the slices (θ). To solve these equations, an iterative method is required until a limiting error is reached. Once Fs and θ are calculated, the remaining forces are also obtained for each slice. Spencer's method is considered to be very accurate and suitable for almost all kinds of slope geometries and may be the most complete equilibrium procedure. It may also be the easiest method for obtaining the Fs (Duncan and Wright, 2005). Depending on the type of slope analysed, this model is able to establish the failure curve following the typical rotational circle, among other uses (Verruijt, 2010).

The Fs, classically defined as a ratio of stabilizing and destabilizing forces, determines the stability of a slope as follows:

(4) F S = ( Forces standing against / oppose sliding ) ( Forces that induce sliding t ) .

According to limit equilibrium methods, the two equilibrium conditions (forces and moments) must be satisfied. Taking into account these elements, the Fs is then obtained from the following expression (Spencer, 1967):

(5) F s = 1 W sin α [ c b sec α + tan ϕ ( W cos α - u b sec α ) ] ,

where ϕ is the friction angle at the fracture surface, u is the pore pressure at the fracture zone, c is the soil cohesion, α is the angle at the base of the slice, W is the external vertical forces and b is the width of the slice. According to Eqs. (4) and (5), the slope can be considered unstable if its value of the safety factor Fs is lower than 1 or stable if it is equal to or higher than 1. It should be noted that, when applying the factor in the engineering and architecture fields, the limiting value tends to be higher than 1, with common values being 1.2 or even up to 1.5 (Burbano et al., 2009) and security coefficients that include the European technical regulations, specifically the technical regulations of Spanish application (Table 2.1 of the DB-C of the CTE, or technical code of the building), among others. This is just a confidence measure for your calculations. The Fs can also be defined as the ratio between the shear strength (τ), based on the cohesion and the angle of friction values, and the shear stress, based on the cohesion and the internal friction angle required to maintain the equilibrium (τmb).

As mentioned, the minimum Fs for considering a slope stable is equal to 1. However, several authors (Yong et al., 1977; Van Westen and Terlien, 1996) suggest that the angle of a slope would have to be defined by a value of the Fs superior to the unity to take into account the exogenous factors of the slope. Following Jiménez Salas (1981), a value of FS≥1.3 can be considered stable by most standards.

To analyse the slope using Spencer's method, a set of equations must be solved to satisfy the forces and momentum equilibrium and to obtain the Fs. The values of Fs and θ are the unknowns that must be solved. Some authors suggest that the variation in θ can be arbitrary (Morgenstern and Price, 1965), although the effect of these variations in the final value of Fs is minimal. The variation in the angle depends on the soil's ability to withstand only a small intensity of the shear stress.

With that being said, if we assume that the forces between slices are parallel (in other words, that θ is constant), Eqs. (1) and (2) become the same, resulting in

(6) Q = 0 .

The assumption that the forces between slices are parallel gives optimal results for the calculation of the critical safety coefficients in Eq. (5) (Spencer, 1967). To solve these equations, we used the FSOLVE function of the MATLAB software, giving an initial Fs and angle. The FSOLVE function is a tool inside the optimization toolbox from MATLAB that solves systems of non-linear equations. When using this tool, an initial value must be provided to start the calculation.

When solving the normal and parallel forces at the base of the slice of the five acting forces, we obtain (Q), resulting from the forces between slices:

(7) Q = c b F sec α + tan ϕ F W cos α - u b sec α - W sin α cos ( α - θ ) 1 + tan ϕ F tan ( α - θ ) .

In this expression, u is the pore pressure (permanent interstitial pressure) at the base of the slice and the weight of the slice is determined by W. If we assume that the soil is uniform and its density (γ) is also uniform, the weight of a slice of height h and width b can be written as follows:

(8) W = γ b h .

The application of a homogeneous pore-pressure distribution (permanent interstitial pressure) has been included in the model (Bishop and Morgenstern, 1960). In this case, the permanent interstitial pressure on the base of the slice was determined by the following expression:

(9) u = r u γ h .

In this expression, u is the pore pressure (permanent interstitial pressure) at the base of the slice, γ is the density of soil, h is the mean height of slice (if the height is not constant) and the weight of it affects the W evaluation.

The pore pressure will be hydrostatic, defined by the following: u=γw(h-hw), where γw is the saturated density of soil, and h and hw are the difference between saturated and dry height. The calculation of the infiltration factor is calculated with the following equation:

(10) r u = u γ h .

The factor ru is a coefficient of pore pressure (interstitial pressure coefficient), which determines the rain infiltration factor on the slopes. It is well known that the water that infiltrates the soil may produce a modification of the pore pressure, affecting its resistant capacity. This factor may vary from 0 (dry conditions) to 0.5 (saturated conditions). In the article of Spencer (Spencer, 1967), assuming a homogeneous pore-pressure distribution as proposed by Bishop and Morgenstern (1960), the mean pore pressure on the base of the slice can be written like Eq. (7).

This equation is used in our proposed model for calculating the safety factor (substituting the expression of u in Eq. 5).

3 Terrain stability (TS) algorithm and tests

Figure 2 shows the results of applying the terrain stability model to an irregular slope, including the initial and final points of the first failure circle (shown in yellow). This circle corresponds with the initial value introduced by the user into the FSOLVE function. The points of the slope are extracted from a DEM model in ArcGIS 10 (Glennon et al., 2008). The slope height is equal to 15 m, and the soil is considered to be uniform with the following nominal properties: γ=19 500 N m−3, φ=22, c=15 000 N m−2 and u=0 N m−2. For the application example of our algorithm in this section, we have used geotechnical data of a cohesive soil of the flysch type of Gibraltar (González de Vallejo et al., 2002).

Figure 2Idealized cross section of a slope. In this example, the centre coordinates are equal to xc=7 m, yc=14 m, and the lower cut has the slope coordinates (P1 point) equal to xt=0 m and t=0 m, data that the user introduces.


The code works as follows: the initial circular failure curve is plotted using the FPLOT tool, as shown in Fig. 2 (yellow line). In this example, the centre coordinates are equal to xc=7 m, yc=14 m, and the lower cut has slope coordinates (P1 point) equal to xt=0 m, yt=0 m. The Fs obtained was 1.6, which is, in principle, a stable slope. It must be taken into account that the mass susceptible to sliding must be divided into a sufficient number of slices. This value is entered into our code through the parameter N. In the application example of our algorithm, the sliding mass was divided into N=500 slices; this value of N is entered into the code by the user, who decides the value of that parameter. The greater the number of slices in which we divide the sliding mass, the more accurate the calculation will be. For the N=500 slice, we consider it to be a balanced value for an optimal calculation, which relates two fundamental parameters (computer calculation capacity and capacity accuracy).

The next step is to apply Spencer's method to the different breakage surfaces until the curve with the lowest Fs is found, and that will be the critical surface susceptible to a circular slip. To determine the minimal Fs using this model, the algorithm calculates the displacement of the lower cut-off point of the critical slip from the slope as well as the position of the centre of rotation of the critical failure curve. In addition, the user must enter a series of possible circular faults. Then, the user introduces the following constraints into the program: the initial or lower point of the failure curve (P1) in its intersection point with the slope, which may or may not match the origin of the slope analysed. Another restriction is the centre of the failure circle (Xc, Yc), which should initially cut the slope, i.e. the breaking curve must be within the feasible sliding region. With these data, the program automatically draws a first curve, in this case the yellow line in Fig. 3, and calculates the safety coefficient Fs for that initial curve.

Figure 3Results following the application of the software showing the slope profile and surface damage. The Fs and the clearest proof of circular failure are also provided (see the yellow line). P1 coordinates are (0, 0), and P2 coordinates are (38.85, 14.6) in metres.


On the basis of this first curve (yellow line in Fig. 2), the program enforces new restrictions:

  • The curve passes through the origin of slope P1=(0,0).

  • The centre of the possible circles of critical breakage is inside the rectangular box defined as (xboxmin.<xc<xboxmax.; ycboxmin.<yc<yboxmax). Note that the coordinates are entered with the 2-D expression (XY).

Both coordinates of the rotation centre position are free and can change for every circle. From the initial failure curve, characterized by the point x=(xc, yc), the MATLAB “fmincon” function is used to obtain a new critical point (xc*, yc*), where the Fs from the breakage curve is the minimum provided by fmincon. In this example, starting from the initial curve (yellow curve) with point x=(7,14), the TS model provides a new point x*=(4.4910,28.1091,0) with a new Fs, FS=1.45. In this case, the new search has been carried out with the following restrictions in the rectangular box, such as 2 m <xc< 8 m and 16 m <yc< 40 m. These restrictions are imposed in order to determine the critical circle. With all these restrictions, and because of the first calculated curve (the yellow curve), the developed model calculates the critical curve among the number of curves selected by the user (500 in this case), as well as the failure circle centre, by applying the fmincon (MATLAB function). This defines the curve with minimum Fs (Fmin) as the value of Fs (see green curve in Fig. 3). When solving this problem, a critical selection is the lower cut-off point of the slope. According to different authors, such as Verruijt (2010) and Castedo et al. (2012), the selected point is the same as the P1 point.

To complete the second phase in the TS model operation, the effect of rain infiltration must be introduced by the coefficient of the pore-pressure factor ru. In this example, the infiltration factor was introduced at the base of each slice to account for the infiltration and pore pressure at the base of the break surface of the slope. If ru increases, the cohesion of the soil mass of the slope decreases, directly affecting the reduction of the slope's Fs. The result is that a dry slope has FS=1.45, but if including the ru parameter equal to 0.3, the Fs decreases to a value of FS=0.95, which means that an Fs is below the unity, so an unstable circular failure appears (see Fig. 4). Entering the infiltration factor, ru, in Spencer's method to introduce the infiltration effects in slopes, the geotechnical cutting elements of the analysed soil are reduced, also reducing the values of the Fs, both for the initial yellow curve and the optimum green curve (Fig. 3). Note that the initial curve in the run shown in Fig. 4 is different from the one in Fig. 3, as it depends on the data introduced.

Figure 4Outcome of the TS model after the introduction of the infiltration factor, producing an unstable circular failure (Fs=0.95).


We can determine that if this infiltration factor value is small enough, taking into account the safety coefficients, the design may still be adequate, but critical information was missing to calculate this parameter.

To clarify the procedure employed in the suggested algorithm, the flow chart (block diagrams) presented in Fig. 5 demonstrates the calculation and iteration process as implemented in our software.

Figure 5Sequential TS algorithm (block diagrams). Numbers in parentheses refer to numbers in the text.


Our algorithm (software) is more versatile compared to the STB 2010; the model developed here can analyse slope from right to left and vice versa, and the STB 2010 only allows the analysis from right to left. Other software programs, like the STB 2010, use a modified version of Bishop's method, a less accurate methodology than Spencer's method. A modified version of Bishop's method solves only the equilibrium in momentum, while the Spencer method also considers the equilibrium in forces.

Another improvement made by the TS code, in comparison with others, is that the use of Spencer's method allows us to analyse any type of slope and soil profile. In this procedure, we calculated the worst breaking curve by modifying the calculation points.

In the TS model, from the first slip rotational circle obtained in MATLAB, many circles were then calculated using the fmincon function, with some user restrictions. However, other models, like the STB 2010, require the definition of a quadrangular region (to look for the centres of rotational failures) and a point (namely five; see Fig. 9) to define the curve as where the failure must pass. Also, the number of circles that the STB 2010 model can analyse for their minimum value is limited to 100.

The TS model can detect relevant earth movements derived from rainfall infiltration, both translational and rotational types (Stead et al., 2006), such as those that usually occur in regions like India, the US, South America and the UK, among other places. The programs that do not contemplate this option will overestimate the Fs, potentially with great errors.

The TS model has an additional advantage: it also offers the opportunity to incorporate, in the same code, the stability analysis and the effect of the infiltration factor in the rainfall regime. This is a step forward from open-access programs, such as STB 2010, and also from alternative payment software, such as Slide.

4 Example of this application in the municipality of La Viñuela, Málaga, Spain

In 2010, La Viñuela, Málaga, (Spain) experienced torrential rainfall. The main consequence was a devastating landslide with serious personal and material losses, as shown in Fig. 6. The coordinates where this event occurred were in degrees (36.88371409801, −4.204982221126).

Figure 6(a) Spanish map with the location of La Viñuela (Google Maps). (b) Real images taken by the authors at La Viñuela in 2010.

4.1 Geological and hydrological environment

The study area is located in the county of La Viñuela, specifically in the Hundidero village, which is located immediately north of the swamp of La Viñuela (El Hundiero) and south of the Baetic System mountain ranges (southern Iberian Peninsula).

According to the Cruden and Varnes' (1996) classification, the slide corresponds to a rotational slide-like complex movement because it was generated in two sequences at different speeds. This type of mechanism is characteristic of homogeneous cohesive soils, as was the one analysed here (Cornforth, 2005; Rahardjo et al., 2007; Lu and Godt, 2008).

This event caused serious damage to different buildings. Regarding the damage caused, in the initial stretch of the slope (its head), a house was dragged and destroyed and another was seriously damaged. On the right bank of the mentioned house, another building was affected. In total, this event left a total of two buildings destroyed and one seriously compromised. Although 15 people lived in these houses, there were no fatalities. About 20 houses were to be constructed at the head of the slope; fortunately, the event happened before this construction. Figure 7 shows an aerial picture from 2006 before the disaster as well as the affected area and landslide in 2010.

Figure 7(a) An aerial photograph from before the event (2006). (b) An aerial photograph taken after the landslide (2010).


4.2 Event features and geometry

For this example, we used data of IGN, the Spanish National Geographic Institute (, last access: 11 December 2017), and a downloaded bitmap MTN25, which is a 1:25 000 topographic map in ETRS 89 coordinates and Universal Transverse Mercator (UTM) projection. The downloaded map is generated in a file by means of a geo-referenced digital rasterization (vector to raster conversion). Specifically, we downloaded page number 1039, which is the one corresponding to the landslide zone of the case study. Figure 8 shows the area of the case study.

Figure 8(a) Topographic map in a GIS map from page number 1039 of the IGN (Spanish National Geographic Institute).

From this map we obtained the topographic information to acquire all necessary profiles to study the landslide. Moreover, as our algorithm is a 2-D model, with this topographic map we study the critical curve of the slip in the most unfavourable profile of the landslide (Fig. 8).

It is well known that mass movements, such as landslides, are highly complex morphodynamic processes. We selected The Hundidero as our study area because it is prone to landslides. In order to analyse this case study using our model, we first calculated the initial displaced volume of the study area. According to the dimensions of the problem, the initial displaced volume was calculated, equivalent to the volume of half an ellipsoid (Varnes, 1978; Beyer, 1987; Cruden and Varnes, 1996) that has Vol =1∕6π (width × length × depth). In our particular case, the width was equal to 70 m, the length was equal to 235 m, and the depth was equal to 5 m, making up a total volume of 4.364 m3 (Fig. 9). Taking an average of 33% elongation, as proposed by Nicoletti and Sorriso-Valvo (1991) and Cruden and Varnes (1996), we determined that the total material displaced in this landslide had an approximate volume of 5.804 m3. In this mass displacement, it is also necessary to consider material added by erosion and dragged from the initial mass displaced. In Fig. 7, the straight line indicates the first rotational movement, and the zigzag line shows the planar drag and glide after the first rotational movement. The green region is the total area displaced or affected by mass movement. After the first circular movement, the mass moved rapidly, associated with a continuous rise in incremental pore pressure and the rapid reduction of shear strength, without allowing pressure dissipation.

Figure 9Characterization and longitudinal section of the rotational sliding (Geolen Engineering, 2010). The location of the dragged house is noted in red, which was analysed by the TS model.


The initial spit of land had an approximate size of 235 m in length by 70 m in width. Due to this initial displacement, there was a drag and a huge posterior planar displacement of about 550 m length, affecting a zone with several parcels of land and buildings. These sizes were confirmed using aerial photography and field data. The soil is basically composed of clays of variable thicknesses and of fine grain, with fluvial sediments and silty clay. The authors obtained these data by conducting a field survey as well as through the laboratory tests carried out by the laboratory Geolen S.A. (Geolen Engineering, 2010). From a geological and geotechnical point of view, according to a survey of those present as the laboratory extracted the materials, different lithological levels can be distinguished, as shown in Table 1.

Table 1Lithology of the area affected by the failure, according to the laboratory tests of Geolen S.A. No groundwater level was detected.

Download Print Version | Download XLSX

The laboratory tests included a sieve analysis (following UNE 103 101) in three of the samples extracted from the field at depths of 1.80–2.00 m, of which 70.3 % were composed of clay and silt; according to this, the sample is classified as cohesive. The liquid limit and the plastic limit were determined in two of the samples (following UNE 103 103 and UNE 103 104, respectively), yielding liquid limit values of 57.5 % and 64.2 %, respectively, and a plasticity index of 37 %. According to the lab results, the material can be classified as high plasticity material with the potential of having a high water content. The landslide analysed began in February 2010, ending in March of that same year. However, based on the field inspection and the analysis of the rainfall series in the La Viñuela region in 2010 (see Fig. 10), it can be inferred that the main causes of the event were the following:

  • the poor geomechanical parameters of the material that formed the affected hillside,

  • the hydrometeorological conditions in the days preceding and days after the event, according to the histogram.

Figure 10Rainfall histogram at La Viñuela from August 2009 to April 2010. Rainfall data have been provided by the Spanish Meteorological Agency (station of Viñuela).


Most of the landslides observed during these days occurred as a consequence of exceptionally intense rainfalls. The precipitation data were provided by the meteorological station of La Viñuela (Fig. 10). It can be observed that large amounts of precipitation fell during the months of December, January, February and March of 2010, with peaks being, at the most, 60 L m−2 in a single day (January and February). In total, 890 L m−2 fell in the 2009–2010 hydrological cycle, which ended at the end of April 2010. This is a key point in slope stability to consider when dealing with areas capable of having high infiltration rates.

The rotational slide analysed had occurred between level 2 and level 3, when the water content reached that depth, as confirmed by the infiltration calculations in the terrain (see graphs in Fig. 9, reaching depths of up to 5 m). Two direct shear tests (consolidated and drained) were conducted in unaltered samples extracted from the boreholes at 3.00–3.60 and 4.00–4.60 m deep. The cut-off values of the soil are specified in Table 2. Those values were used in the developed software to obtain the safety coefficient and the theoretical failure curve.

Table 2Summary chart of the characteristics of the soil analysed at the Geolen S.A. laboratory: φ is the angle of internal friction, c is the cohesion, γSat is the saturated specific gravity and γa is the apparent specific gravity.

Download Print Version | Download XLSX

The dynamic and continuous tests were carried out by the Geolen S.A. laboratory with an automatic penetrometer of the ROLATEC ML-60-A type. The data obtained were transcribed by the number of strokes to advance the 20 cm tip, which is called the “penetration number” (N20).

This test is included in the ISO 22476-2:2005 standard as a dynamic probing super heavy and consists of penetrating the ground with a conical tip of standard dimensions. The depth of the failed mass can be estimated as well as the theoretical failure curve for an increase in the soil consistency (see data in Table 3).

Table 3Summary chart of the soil analysed at the Geolen S.A. laboratory. Bold values show, according to the data of the field penetrometers, the depth mobilized by the rotational sliding.

Download Print Version | Download XLSX

The change in the geomechanical response of the soil takes place at a depth of 4–5 m, according to the results of N20 and US (samples without changes) taken along the analysed column. In this case, the sloped ground mass showed a characteristic striking relationship of a displaced terrain (Gonzalez de Vallejo et al., 2002). This differs from the underlying or unmoved terrain, which indicated a more consistent striking relationship that was taken within the area of the landslide behind the house drawn in accordance with the analysis of the hits N20 from Table 3.

Figure 11(a) Hydraulic potential. (b) Volumetric water content. Both have been plotted as a function of the depth (mm) at different times (d).


4.3 Input data

To analyse the topography of the critical section, we obtained the DEM data from the ArcGIS 10 software program (Environmental System Research, 2017), with a scale of 1:1000, through IGN raster maps, with adequate accuracy. These data were interpolated to a 2 m grid using a triangulated network interpolation methodology. Orthophotos proved to be very useful in locating the landslide with accuracy and to validate the field survey. The model developed here applies to failure in an initiation zone, in addition to predicting landslides, including those induced by the infiltration of critical rains.

To complete the input data, we plotted the hydraulic potential and the volumetric water content, as a function of depth in the ground for different time steps, using a previously developed infiltration model, as shown in Fig. 11 (Herrada et al., 2014). The figure shows the evolution of how the wetting front advances can be observed. These reached almost 5 m depth at the end of April 2010.

4.4 Analytical results

We applied the TS model using topographic data obtained from the ArcGIS 10 software program. We did so to obtain the degree of stability of the sliding land based on the angle of internal friction, the cohesion, the density and the angle of the slope we analysed. Figure 9 shows the analytical results from the real slope by studying and analysing the most unfavourable profile of the landslide studied. In addition we compared the results given by the developed TS model and the results given by STB 2010 model, using free surfaces in both cases. In our model the worst curve (shown in green) was calculated automatically from the initial curve (shown in blue), resulting in FS=2.300, in the dry state (Fig. 12).

Figure 12(a) TS model with a critical failure of FS=2.300. (b) Results are from the STB 2010 model with an Fs of 2.063.


As can be noted, the failure curves are similar, and the safety coefficients Fs only differ by 0.237. In both cases, the results indicated are conservative estimates, resulting in a stable slope that was not realistic, as was the case in La Viñuela. In order to get the most unfavourable curve, which would match the analysis of the actual event, the pore coefficient must be introduced. At the first runs of the model, the ru was equal to zero (dry soil in Fig. 9), but if this value is changed to ru=0.35, the results are quite different (Fig. 13). The resulting failure was near the surface and the top cut with the slope found relatively near to the houses. Taking into account the infiltration of rainwater, the slope analysed in the TS model showed a value of FS=0.98, in other words showing that it was unstable.

Figure 13A new calculation including the pore coefficient ru, showing the worst curve in green. The circles show the houses dragged by the landslide.


This calculation and the theoretical failure curve provided by our model was able to reproduce, in a realistic way, the landslide which occurred in La Viñuela. Our model found that the critical surface area that corresponded with the profile of the terrain was 12 927.45 m2, which closely matches the real situation. In the STB 2010 program, it was 7825.35 m2; therefore, our prediction was more accurate.

As mentioned, the STB 2010 model does not allow stability calculations to apply to rainfall infiltration on a hillside. Hence, it is not capable of predicting a hillside's instability in a critical rainfall scenario, which was critical in the slope analysed. The STB 2010 model found that the hillside studied had an Fs of FS=2.063; that means that it was a very stable slope. Consequently, our original algorithm TS model appears to be more efficient and accurate.

If we compare the results of the penetrometric tests (Table 3) and the laboratory tests (Geolen Engineering, 2010) summarized in the actual critical surface in the most unfavourable profile of landslide (Fig. 9) with those offered by our TS algorithm (Fig. 13) to which we apply the infiltration factor ru=0.35 (high interstitial pressure), we can check the similarity between the two critical surfaces of the landslide.

A value of ru=0.35 was introduced in the calculation, and the code gave us a value of the slope safety factor of Fs=0.95 (unstable); when in the dry state, the code calculated a safety factor of Fs=2.300 (stable). The calculation of the safety factor in the STB 2010 program that lacks the analysis of infiltration in the calculation offered a result of Fs=2.063 (stable).

Using the STB 2010 program, we would not have been able to previously detect the landslide of the case study of the paper, calculation that is not normally done in the stability calculations; with the calculation with our code we could have avoided the collapse of the building.

With these results, the terrain stability analysis performed using the developed model defines the slip-breaking curve that intuitively appears to be susceptible to failure fairly well, especially when heavy rains occur. For example, the landslides which occurred in the La Viñuela area could only have been predicted if the infiltration had been taken into account. Even then, it could not have been done with other available software programs, which were not able to consider it.

5 Conclusion

The terrain stability (TS) analysis defines the critical surface to landslides in 2-D of each profile of the analysed slope and the safety slip factor (Fs) fairly well. We developed this model due to the need for a useful tool to predict landslides, especially when heavy rains occur.

The TS model we developed uses Spencer's method, which is more precise than the modified Bishop method; this model is used by other software such as the case of the STB 2010, so it differs in the results it provides for the Fs. It also takes into account the factor of water infiltration due to critical rains, which other software programs do not consider. A failure surface can be determined by constraints using the MATLAB function fmincon. The data needed to run the model include soil and climate properties that may vary in space and time. The exit indices of the analysis (Fs) should be interpreted in terms of relative risk. The methods implemented in the TS model are based on data structures, which are based on the data entry of the elevation model (DEM), so we obtain a topographic map, a key element in obtaining the topographic profile to be studied with our algorithm.

In the case study analysed, the slope was initially stable and was determined by the analysis performed with the STB 2010 model. However, the slope became unstable due to the heavy rains of that hydrological period, which called for the application of the pore-pressure coefficient ru. For analysing cases of heavy rain, this model is a powerful tool for determining slope stability. In addition, thanks to the great versatility of this model, it is applicable to any analysis in other parts of the world, based on the methods of limit equilibrium (Spencer, 1967). The TS model can also be used in combination with GIS software, SINMAP, the TRIGRS model and aerial photographic analysis, and mapping techniques, or even as a part of other models like the coastal recession models (Castedo et al., 2012).

Data availability

The data that can be publicly accessed have been described within the text; however, raw data belong to the La Viñuela Municipality, which contracted the first author to conduct a field survey and analysis. Authors are not authorized to publicly share the data directly. However, some of them can be shared upon request.

Author contributions

Each author has made substantial contributions to the work. AGM contributed to the conception of the work, to the applied methodology, to acquisition, to formal analysis, to data elaboration and to writing the original draft of the paper. MAH contributed to the code development. JIY and RC contributed to the field survey, to formal analysis, to validation of the work and to writing the original draft of the paper. Each author has approved the submitted version and agrees to be personally accountable for their own contributions and for ensuring that questions related to the accuracy or integrity of any part of the work, even ones in which the author was not personally involved, are appropriately investigated, resolved and documented in the literature.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “Advances in computational modelling of natural hazards and geohazards”. It is a result of the Geoprocesses, geohazards – CSDMS 2018, Boulder, USA, 22–24 May 2018.

Review statement

This paper was edited by Albert J. Kettner and reviewed by three anonymous referees.


Acharya, K. P., Bhandary, N. P., Dahal, R. K., and Yatabe, R.: Seepage and slope stability modelling of rainfall-induced slope failures in topographic hollows, Geomat. Nat. Hazards Risk, 7, 721–746,, 2016a. 

Acharya, K. P., Yatabe, R., Bhandary, N. P., and Dahal, R. K.: Deterministic slope failure hazard assessment in a model catchment and its replication in neighbourhood terrain, Geomat. Nat. Hazards Risk, 7, 156–185,, 2016b. 

Aleotti, P.: A warning system for rainfall-induced shallow failures, Eng. Geol., 73, 247–265, 2004. 

Alvioli, M. and Baum, R. L.: Parallelization of the TRIGRS model for rainfall-induced landslides using the message passing interface, Environ. Model. Softw. 81, 122–135,, 2016. 

Anbalagan, R.: Landslide hazard evaluation and zonation mapping in mountainous terrain, Eng. Geol., 32, 269–277, 1992. 

Anderson, M. G. and Howes, S.: Development and application of a combined soil water-slope stability model, Q. J. Eng. Geol. Lond., 18, 225–236, 1985. 

Ayalew, L., Yamagishi, H., and Ugawa, N.: Landslide susceptibility mapping using GIS-based weighted linear combination, the case in Tsugawa area of Agano River, Niigata Prefecture, Japan, Landslides 1, 73–81, 2004. 

Baum, R. L., Savage, W. Z., and Godt, J. W.: TRIGRS-A Fortran program for transient rainfall infiltration and grid-based regional slope-stability analysis, Version 2.0, US geological survey open-file report 424, US Geological Survey, available at: (last access: 14 May 2018), 2008. 

Beyer, W. H.: Handbook of Mathematical Sciences, 6th Edn., CRC Press, Boca Raton, Florida, 1987. 

Bishop, A. W.: The Use of the slip circle in the Stability Analysis of Slope, Geotechnique, 5, 7–16, 1955. 

Bishop, A. W. and Morgenstern, N. R.: Stability coefficients for earth slope, Geotechnique, 10, 129–150, 1960. 

Bisson, M., Spinetti, C., and Sulpizio, R.: Volcaniclastic flow hazard zonation in the sub-apennine vesuvian area using GIS and remote sensing, Geosphere, 10, 1419–1431,, 2014. 

Borselli, L.: SAPP 4.2.0.: Advanced 2D Slope stability Analysis by LEM by SSAP software, SSAP code Manual, version 4.2.0, available at: http://www.Ssap.Eu/Manualessap2010.Pdf, last access: 20 September 2012. 

Burbano, G., del Cañizo, L., Gutiérrez, J. M., Fort, L., Llorens, M., Martínez, M., Paramio, J. R., and Simic, D.: Guía de cimentaciones en obras de carretera, Ministerio Fomento, Madrid, 2009. 

Canili, E., Mergili, M., Thiebes, B., and Glade, T.: Probabilistic landslide ensemble prediction systems: lessons to be learned from hydrology, Nat. Hazards Earth Syst. Sci., 18, 2183–2202,, 2018. 

Carrara, A., Cardinali, M., Guzzetti, F., aand Reichenbach, P.: GIS technology in mapping landslide hazard, in: Geographical Information System in Assessing Natural Hazard, edited by: Carrara, A. and Guzzetti, F., Kluwer Academic Publisher, the Netherlands, 135–175, 1995. 

Casagli, N., Catani, F., Puglisi, C., Delmonaco, G., Ermini, L., and Margottini, C.: An inventory-based approach to landslide susceptibility assessment and its application to the Virginio River Basin, Italy, Environ. Eng. Geosci., 3, 203–216, 2004. 

Castedo, R., Murphy, W., Lawrence, J., and Paredes, C.: A new process–response coastal recession model of soft rock cliffs, Geomorphology, 177, 128–143, 2012. 

Cheng, Y. M., Lansivaara, T., and Wei, W. B.: Two-dimensional slope stability analysis by limit equilibrium and strength reduction methods, Comput. Geotech., 34, 137–150, 2007. 

Chugh, A. K. and Smart, J. D.: Suggestions for slope stability calculations, Comput. Struct., 14, 43–50, 1981. 

Collison, A., Wade, S., Griffiths, J., and Dehn, M.: Modelling the impact of predicted climate change on landslide frequency and magnitude in SE England, Eng. Geol., 55, 205–218, 2000. 

Cornforth, D. H.: Landslides in Practice: Investigation, Analysis, and Remedial/Preventative Options in Soils, John Wiley, Hoboken, NJ, 2005. 

Crosta, G. B. and Frattini, P.: Distributed modelling of shallow landslides triggered by intense rainfall, Nat. Hazards Earth Syst. Sci., 3, 81–93,, 2003. 

Crozier, M. J. and Glade, T.: Landslide hazard and risk: issues, concepts, and approach, in: Landslide Hazard and Risk, edited by: Glade, T., Anderson, M., and Crozier, M., Wiley, Chichester, 1–40, 2005. 

Cruden, D. M. and Varnes, D. J.: Landslides types and processes, in: Landslides investigation and mitigations, Transportation Research Board Special report 24, edited by: Schuster, R. L. and Turner, A. K., National Academy Press, Washington, 36–75, 1996. 

Dai, F. C. and Lee, C. F.: Terrain-based mapping of landslide susceptibility using a geographical information system: a case study, Can. Geotech. J., 38, 911–923, 2001. 

Dai, F. C., Lee, C. F., and Ngai, Y. Y.: Landslide risk assessment and management: an overview, Eng. Geol., 64, 65–87, 2002. 

Duncan, J. M. and Wright, S. G.: The accuracy of equilibrium methods of slope stability analysis, Eng. Geol., 16, 5–17, 1980. 

Duncan, J. M. and Wright, S. G.: Soil Strength and Slope Stability, John Wiley, Hoboken, NJ, 2005. 

Environmental Systems Research: Institute Geographic information system (platform and resources ArcGIS), California, EEUU, available at:, last access: 4 September 2017. 

Fall, M., Azzam, R., and Noubactep, C.: A multi-method approach to study the stability of natural slopes and landslide susceptibility mapping, Eng. Geol., 82, 241–263, 2006. 

Fellenius, W.: Calculation of the stability of earth dams, in: Vol. 4, Transactions of 2nd Congress on Large Dams, 7–12 September 1936, Washington, D.C., 445–462, 1936. 

Geolen Engineering: Geotechnical study in the Viñuela, Sevilla, Spain, available at:, last access: 18 May 2010. 

Girma, F., Raghuvanshi, T. K., Ayenew, T., and Hailemariam, T.: Landslide hazard zonation in Ada Berga District, Central Ethiopia – a GIS based statistical approach, J. Geomat., 90, 25–38, 2015. 

Glennon, R., Harlow, M., Minami, M., and Booth, B: ArcGis 9., ArcMap Tutorial, Esri, North Carolina, USA, 2008. 

González de Vallejo, L., Ferrer, M., Ortuno, L., Oteo, C.: Ingeniería Geológica. Madrid: Prentice Hall, 2002. 

Griffiths, D. V.: Slope stability analysis by finite elements. A guide to the use of Program slope 64, Geomechanics Research Center Colorado School of Mines, available at:, last access: 20 January 2015. 

Griffiths, D. V. and Marquez, R. M.: Three-dimensional slope stability analysis by elasto-plastic finite elements, Géotechnique, 57, 537–546, 2007. 

Guha-Sapir, D., Hargitt, D., and Hoyois, G.: Thirty years of natural Disasters 1974–2003: The Numbers, Centre for Research on the Epidemiology of Disasters, available at: (last access: 5 June 2015), 2004. 

Gutiérrez-Martín, A.: El agua de infiltración de lluvia, agente desestabilizador de taludes en la provincia de Málaga. Modelos constitutivos, Doctoral Thesis, University of Granada, Granada, 2015. 

Guzzetti, F., Carrara, A., Cardinali, M., Reichenbach, P.: Landslide hazard evaluation: a review of current techniques and their application in a multi-scale study, central Italy, Geomorphology, 31, 181–216, 1999. 

Guzzetti, F., Peruccacci, S., Rossi, M., and Stark, C. P.: Rainfall thresholds for the initiation of landslides in central and southern Europe, Meteorol. Atmos. Phys., 98, 239–267, 2007. 

Herrada, M. A., Gutiérrez-Martin, A., and Montanero, J. M.: Modeling infiltration rates in a saturated/unsaturated soil under the free draining condition, J. Hydrol., 515, 10–15, 2014. 

Hoek, E. andBray, J. W.: Rock Slope Engineering (revised thirded.), Inst. of Mining and Metallurgy, London, 1981. 

Iverson, R. M.: Landslide triggering by rain infiltration, Water Resour. Res., 36, 1897–1910, 2000. 

Janbu, N.: Stability analysis of slopes with dimensionless parameters, in: vol. 46, Harvard University soil mechanics series, Harvard University, Cambridge, Massachusetts, 1954. 

Jia, G. J., Tian, Y., Liu, Y., and Zhang, Y.: A static and dynamic factors-coupled forecasting model of regional rainfall-induced landslides: A case study of Shenzhen, Sci. China: Technol. Sci., 51, 164–175, 2008. 

Jiménez Salas, J. A. and Justo Alpañes, J. L.: Geotecnia y Cimientos II, Rueda, Madrid, 1981. 

Johari, A. and Mousavi, S.: An analytical probabilistic analysis of slopes based on limit equilibrium methods, Bull. Eng. Geol. Environ., 15, 1–15,, 2018. 

Khan, M. E.: The Death Toll from Natural Disasters: The Role of Income, Geography, and Institutions, Rev. Econ. Stat., 87, 271–284, 2005. 

Leroi, E.: Landslide risk mapping: problems, limitation and developments, in: Landslide Risk Assessment, edited by: Cruden, F., Balkema, Rotterdam, 239–250, 1997. 

Liu, S. Y., Shao, L. T., and Li, H. J.: Slope stability analysis using the limit equilibrium method and two finite element methods, Comput. Geotech., 63, 291–298, 2015. 

Lu, N. and Godt, J.: Infinite slope stability under steady unsaturated seepage conditions, Water Resour. Res., 44, W11404,, 2008. 

Martelloni, G. and Bagnoli, F.: Infiltration effects on a two-dimensional molecular dynamics model of landslides, Nat. Hazards, 73, 37–62, 2014. 

Martelloni, G., Segoni, S, Fanti, R., andCatani, F.: Rainfall thresholds for the forecasting of landslide occurrence at regional scale, Landslides, 9, 485–495,, 2011. 

Martelloni, G., Bagnoli, F., and Guarino, A.: A 3D model for rain-induced landslides based on molecular dynamics with fractal and fractional water diffusion, Commun. Nonlin. Sci. Numer. Simul., 50, 311–329, 2017. 

Mergili, M., Marchesini, I., Rossi, M., Guzzetti, F., and Fellin, W.: Spatially distributed three dimensional slope stability modelling in a raster GIS, Geomorphology, 206, 178–195,, 2014a. 

Mergili, M., Marchesini, I., Alvioli, M., Metz, M., Schneider-Muntau, B., Rossi, M., and Guzzetti, F.: A strategy for GIS-based 3-D slope stability modelling over large areas, Geosci. Model Dev., 7, 2969–2982,, 2014b. 

Michel, G. P., Kobiyama, M., and Fabris, R.: Comparative analysis of SHALSTAB and SINMAP for landslide susceptibility mapping in the Cunha River basin, southern Brazil, J. Soils Sediments, 7, 1266–1277, 2014. 

Michel, G. P., Kobiyama, M., and Fabris, R.: Critical rainfall to trigger landslides in Cunha River basin, southern Brazil, Nat. Hazards, 75, 2369–2384, 2015. 

Montgomery, D. and Dietrich, W.: R-SHALSTAB: A digital terrain model for mapping shallow landslide potential, to be published as a technical report by NCASI, available at:, (last access: 14t May 2018), 1998. 

Morgenstern, N. R. and Price, V. E.: The analysis of the stability of general slip surfaces, Geotechnique, 15, 79–93, 1965. 

Nicoletti, P. G. and Sorriso-Valvo, M.: Geomorphic controls of the shape and mobility of rock avalanches, GSA Bulletin, 103, 1365–1373, 1991. 

Pack, R. T., Tarboton, D. G., and Goodwin, C. N.: Assessing Terrain Stability in a GIS using SINMAP, in: 15th annual GIS conference, GIS 2001, 19–22 February 2001, Vancouver, British Columbia, 2001. 

Parise, M. and Jibson, R. W.: A seismic landslide susceptibility rating of geologic units based on analysis of characteristics of landslides triggered by the 17 January, 1994 Northridge, California earthquake, Eng. Geol., 58, 251–270, 2000. 

Raghuvanshi, T. K., Negassa, L., and Kala, P. M.: GIS based grid overlay method versus modeling approach – a comparative study for Landslide Hazard Zonation (LHZ) in Meta Robi District of West Showa Zone in Ethiopia, Egypt, J. Remote Sens. Space Sci., 18, 235–250, 2015. 

Rahardjo, H., Ong, T. H., Rezaur, R. B., and Leong, E. C.: Factors controlling instability of homogeneous soil slopes under rainfall, J. Geotech. Geoenviron. Eng., 133, 1532–1543, 2007. 

Raia, S., Alvioli, M., Rossi, M., Baum, R. L., Godt, J. W., and Guzzetti, F.: Improving predictive power of physically based rainfall-induced shallow landslide models: a probabilistic approach, Geosci. Model Dev., 7, 495–514,, 2014. 

Ramos Vásquez, A. A.: Análisis de estabilidad de taludes en rocas, Simulación con LS-DYNA y comparación con Slide, Trabajo Fin de Máster, Máster Universitario en Ingeniería Geológica, ETSI Minas y Energía, Universidad Politécnica de Madrid, Madrid, 2017. 

Reid, M. E., Christian, S. B., Brien, D. L., and Henderson, S. T.: Scoops-3D – Software to analyze Three-Dimensional Slope Stability Throughout a Digital Landscape, Version 1.0, US Geological Survey, Virginia, 2015. 

Rigon, R., Bertoldi, G., and Over, T. M.: 2GEOtop: A distributed hydrological model with coupled water and energy budgets, J. Hydrometeorol., 7, 371–388,, 2006. 

Rossi, G., Catani, F., Leoni, L., Segoni, S., and Tofani, V.: HIRESSS: a physically based slope stability simulator for HPC applications, Nat. Hazards Earth Syst. Sci., 13, 151–166,, 2013. 

Sharma, S. and Moudud, A.: Interactive slope analysis using Spencer's method, Stability and Performance of Slopes and Embankments II: Proceeding of a Specialty Conference Sponsored by ASCE, Geotechnical Special Publication 31, ASCE, 2, 506–520, 1992. 

Simoni, S., Zanotti, F., Bertoldi, G., and Rigon, R.: Modelling the probability of occurrence of shallow landslides and channelized debris flows using GEOtop-FS, Hydrol. Process., 22, 532–545,, 2008. 

SLIDE V5: 2D limit equilibrium slope stability for soil and rock slopes, user's guide, available at:, last access: 1 May 2018. 

Spencer, E.: A method of analysis of analysis of the stability of embankments assuming parallel interslice forces, Géotechnique, 17, 11–26, 1967. 

Stead, D., Eberhardt, E., and Coggan, J. S.: Developments in the characterization of complex rock slope deformation and failure using numerical modelling techniques, Eng. Geol., 83, 217–235, 2006. 

Toya, H. and Skidmore, M.: Economic development and the impacts of natural disasters, Econ. Lett., 94, 20–25,, 2007. 

Tran, T. V., Alvioli, M., Lee, G., and An, H. U.: Three-dimensional, time-dependent modelling of rainfall-induced landslides over a digital landscape: a case study, Landslides, 15, 1071–1084,, 2018. 

Tschuchnigg, F., Schweiger, H. F., and Sloan, S. W.: Slope stability analysis by means of finite element limit analysis and finite element strength reduction techniques. Part II: Back analyses of a case history, Comput. Geotech., 70, 178–189, 2015. 

Van Westen, C. J. and Terlien, M. J. T.: An approach towards deterministic landslide hazard analysis in gis. A case study from Manizales (Colombia), Earth Surf. Proc. Land., 21, 853–868, 1996. 

Varnes, D. J.: Slope movement types and processes, in: Landslides: analysis and control, Transportation Research Board. Special report 176, edited by: Schuster, R. L. and Krizek, R. J., Transportation Research Board, National Academy of Sciences, Washington, D.C., 11–33, 1978. 

Varnes, D. J.: Landslide Types and Processes, in: Landslides: Investigation and Mitigation, Transportation Research Board Special Report 247, edited by: Turner, A. K. and Schuster, R. L., National Academy Press, National Research Council, Washington, D.C., 1996. 

Verruijt, A.: STB – SLOPE: Stability Analysis Program, Delft University, available at: (last access: 2 June 2017), 2010. 

Wang, X. and Niu, R.: Spatial forecast of landslides in three gorges based on spatial data mining, Sensors, 9, 2035–2061, 2009. 

Wilkinson, P. L., Anderson, M. G., Lloyd, D. M., and Renaud, P. N.: Landslide hazard and bioengineering: towards providing improved decision support through integrated numerical model development, Environ. Model. Softw., 17, 333–344, 2002.  

Wilson, R. C. and Jayko, A. S.: Preliminary maps showing rainfall thresholds for debris-flow activity, San Francisco Bay Region, California, US Geological Survey Open-File Report 97-745 F, US Geological Survey, Menlo Park, California, 1997. 

Wu, W. and Sidle, R. C.: A Distributed Slope Stability Model for Steep Forested Basins, Water Resour. Res., 31, 2097–2110,, 1995. 

Yong, R. N., Alonso, E., Tabba, M. M., and Fransham, P. B.: Application of Risk Analysis to the Prediction of Slope Stability, Can. Geotech. J. 14, 540–553, 1977. 

Zhang, S., Zhao, L., Delgado-Tellez, R., and Bao, H.: A physics-based probabilistic forecasting model for rainfall-induced shallow landslides at regional scale, Nat. Hazards Earth Syst. Sci., 18, 969–982,, 2018. 

Zhou, X. P. and Cheng, H.: Analysis of stability of three-dimensional slopes using the rigorous limit equilibrium method, Eng. Geol., 160, 21–33, 2013. 

Zhu, D. Y., Lee, C. F., Qian, Q. H., and Chen, G. R.: A concise algorithm for computing the factor of safety using the Morgenstern–Price method, Can. Geotech. J., 42, 272–278,, 2005. 

Short summary
This research work of this paper completes the authors' activity in the Military Emergencies Unit of Spain in responding to the problem of slope instabilities and how to predict them in the case of heavy rains. This work completes the work of stabilization and stresses made by torrential rains in the south of Spain. The results have been satisfactory.
Final-revised paper