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

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.

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).

722
A. Gutiérrez-Martín et al.: Development and validation of the terrain stability model tain 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 TRI-GRS, 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 r u 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 F s (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 wellknown 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 F s 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 de-veloped a new model and programed it using MATLAB. The primary result of this model was a stability index, namely the minimum F s , 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.

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): 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. These equations must be solved to get the F s 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 F s 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 F s (Duncan and Wright, 2005). Depending on the type Figure 1. Representation 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).
of slope analysed, this model is able to establish the failure curve following the typical rotational circle, among other uses (Verruijt, 2010).
The F s , classically defined as a ratio of stabilizing and destabilizing forces, determines the stability of a slope as follows: According to limit equilibrium methods, the two equilibrium conditions (forces and moments) must be satisfied. Taking into account these elements, the F s is then obtained from the following expression (Spencer, 1967): 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 F s 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 F s 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 F s 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 F s superior to the unity to take into account the exogenous factors of the slope. Following Jiménez Salas (1981), a value of F S ≥ 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 F s . The values of F s 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 F s 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 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 F s 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: 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: 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: 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−h w ), where γ w is the saturated density of soil, and h and h w are the difference between saturated and dry height. The calculation of the infiltration factor is calculated with the following equation: The factor r u 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).
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 x c = 7 m, y c = 14 m, and the lower cut has slope coordinates (P1 point) equal to x t = 0 m, y t = 0 m. The F s 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 F s is found, and that will be the critical surface susceptible to a circular slip. To determine the minimal F s 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 (P 1 ) 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 (X c , Y c ), 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 F s for that initial curve.
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 P 1 = (0, 0).
-The centre of the possible circles of critical breakage is inside the rectangular box defined as (x box min. < x c < x box max. ; y c box min. < y c < y box max ). Note that the coordinates are entered with the 2-D expression (X, Y ).
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 = (x c , y c ), the MATLAB "fmincon" function is used to obtain a new critical point (x * c , y * c ), where the F s 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 F s , F S = 1.45. In this case, the new search has been carried out with the following restrictions in the rectangular box, such as 2 m < x c < 8 m and 16 m < y c < 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 F s (F min ) as the value of F s (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 P 1 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 r u . 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 r u increases, the cohesion of the soil mass of the slope decreases, directly affecting the reduction of the slope's F s . The result is that a dry slope has F S = 1.45, but if including the r u parameter equal to 0.3, the F s decreases to a value of F S = 0.95, which means that an F s is below the unity, so an unstable circular failure appears (see Fig. 4). Entering the infiltration factor, r u , 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 F s , 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.
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.
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 F s , 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.

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).

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.

Event features and geometry
For this example, we used data of IGN, the Spanish National Geographic Institute (http://centrodedescargas.cnig.es/ CentroDescargas, 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 correspond-ing to the landslide zone of the case study. Figure 8 shows the area of the case study. 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 m 3 (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 m 3 . 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.
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.
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.
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.
A. Gutiérrez-Martín et al.: Development and validation of the terrain stability model Figure 9. Characterization 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.  Table 2. Summary 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.

Soil
Result Units parameter 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" (N 20 ).
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). The change in the geomechanical response of the soil takes place at a depth of 4-5 m, according to the results of N 20 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 land-slide behind the house drawn in accordance with the analysis of the hits N 20 from Table 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.

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 F S = 2.300, in the dry state (Fig. 12).
As can be noted, the failure curves are similar, and the safety coefficients F s 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 r u was equal to zero (dry soil in Fig. 9), but if this value is changed to r u = 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 F S = 0.98, in other words showing that it was unstable.
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 m 2 , which closely matches the real situation. In the STB 2010 program, it was 7825.35 m 2 ; 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 F s of F S = 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 r u = 0.35 (high interstitial pressure), we can check the similarity between the two critical surfaces of the landslide.
A value of r u = 0.35 was introduced in the calculation, and the code gave us a value of the slope safety factor of F s = 0.95 (unstable); when in the dry state, the code calculated a safety factor of F s = 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 F s = 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.

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 (F s ) 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 F s . 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 (F s ) 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 r u . 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 TRI-GRS 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.