Method to estimate the initial landslide failure surface and volumes using grid points and spline curves in MATLAB

This paper presents a new method to estimate a landslide failure surface and its volume based on a digital elevation model (DEM) using spline curves by assuming values of the tangent of the failure surface along with transversal vertical profiles. The model will give the depth of the probable failure surface plotted using the 2D grid function. We can easily visualize the results by running the codes; the display will show the step-by-step processes involved. The model requires fundamental data inputs, which are readily available online throughout the world, such as a DEM of the study area and a KML file of contour limits of the landslide as input. The calculation procedure is simple; it evaluates the depth of the failure surface at each point of the grid. The cross-section for each point is estimated perpendicular to the slope line joining the highest elevation point and the lowest elevation point within the contour limits. A cubic spline curve is calculated using two endpoints and the first derivative at these points. The Z value is calculated using these parameters at each grid point (centre of a grid), and finally, a 3D failure surface is generated, enabling the volume calculation. MATLAB code has been written to execute this process automatically. This model can easily use input data with very little fieldwork, and the processing can produce a result in a few minutes. We can enrich a real scenario by imposing more constraints on the variation in the cross-section angle based on the limit of the landslide using field measurements. This model is also helpful in predicting the area affected by the mass of the failure surface, and we can determine a predefined hazard zone to avoid any casualties in the future.


Introduction
In hilly terrain with a high annual rate of rainfall, landslides are a common phenomenon (Clague and Stead 2012). To increase human development and create reliable transportation networks, more roads, tunnels, and other infrastructures such as dams and barrages are being built, which can induce landslides. This increase in human activities in hilly areas can trigger further landslides as slopes become unstable because of the removal of base materials, dynamite blasting for tunnels, and the increased pore pressure in surrounding rocks by modifying surface runoff etc. (Froude and Petley 2018;Jaboyedoff et al. 2016). These situations can worsen during rainy seasons. Some natural phenomena, such as earthquakes and other tectonic activities, further increase the chance of landslides (Keefer 1984;Lacroix 2016). To reduce and mitigate the losses due to landslides, necessary steps are required. To achieve this outcome, we need a detailed understanding of landslide behaviour. The characterization of landslide volume is a key element for characterizing their behaviours and assessing their potential impact on objects at risk. The landslide volume, defined by its failure surface, is an important parameter that controls the affected surface area or the propagation distance of a failed mass and thus the potential damages related to this phenomenon (Corominas 1996;Jakob 2005;Jakob et al. 2012;von et al. 2016). Volume calculation is a difficult task because (Jaboyedoff et al. 2020), in most cases, the failed mass generally rests over a potential failure surface, and its propagation is triggered by an event such as an increase in pore pressure because of rainfall, snowmelt, earthquakes, and human activities leading to a decrease in safety factor below 1 (Terzaghi 1950;Jaboyedoff et al. 2016). Without in situ investigations, we cannot acquire reliable information about failure surface geometry. To estimate initial volumes, we typically multiply the area by the average depth of landslides, which is a very rough approximation (Guzzetti et al. 2009). We can make further improvements by dividing the area into small cross-sections and measuring the depth using boreholes, and then multiplying the cross-sectional area by the depth of that portion; additionally, we can make further improvements by using techniques based on surface information. Calculating the volume involved when there are slope movements is a difficult assessment task.
There have been many methods developed to calculate these details; some are based on ellipsoids (Guzzetti et al. 2009;Larsen et al. 2010), and some are based on balanced cross-sections (Hutchinson, 1983;Carter and Bentley 1985), but most of these methods require expertise in that field, and they are based on the available data, such as surface morphological features, borehole data, and geophysical interpretations. Few methods are supported by a systematic approach. With new technologies and the increasing availability of new satellite data, such as high-resolution digital elevation models, ASTER Global DEMs (G-DEMs), and enhanced computer processors, we can improve our older methods and create new methods. A recent method called the sloping local base level (SLBL; Jaboyedoff and Derron 2013) uses the DEM file as input and tries to predict the depth of failure. This paper presents a new method to estimate the landslide failure surface (mainly in the case when whole debris is resting on the failed surface and the surface is not visible) and volume based on a DEM using vertical profiles made of spline curves evaluated at each point of the DEM grid. The model gives the depth of probable failure surface plotted using the 2D grid function. We have provided the MATLAB scripts and the data online on the GitHub platform (https:// github. com/ Gauta m5673/ Failu re_ surfa ce_ estim ator). Users can easily follow the results from starting to end. The necessary data include the limit of the landslide and the slope of the failure surface at the limit of the landslide, along with the profile. This method requires fundamental data inputs, which are readily available online throughout the world, such as a DEM file of the study area and a KML file of contour limits of the landslide as input. The procedure calculates the altitude of the failure surface at each point of the grid (Fig. 1).
The cross-section for each point is created perpendicular to the line joining the highest elevation point and the lowest contour limits. "A cubic spline curve requires four parameters provided by each cross-section endpoint" (Bartels et al. 1987), i.e. the endpoints and their first derivative. There are several ways to define the first derivative, whether constant (method 1) or varying (methods 2 and 3), from top to bottom on each side of the contour. The Z value is calculated at each grid point (here, in this case, the grid points refer to the centre of a grid) using the spline equation. The end of the process leads to a 3D failure surface DEM. The volume is calculated using the elevation difference of the failure surface. To execute this process automatically, MAT-LAB code has been written. We have also checked the computational limits; after testing on different data, we have indicated the conditions in the code. Based on the calculation limits, for normal calculation, the number of cells is 3000; i.e. if we have a landslide with the dimension of 3000 m 2 , it can calculate the data up to 1-m resolution. For the landslides with dimensions in the range of 3000-17,000 m 2 , it can calculate the data up to the 2-m resolution. For landslides with dimensions greater than 17,000 m 2 , it can calculate the data up to 3-m resolution. DEM resolutions lower than these values can be easily calculated without facing any computational load. Furthermore, this method provides an estimation of the failure surface topography; it does not need too high resolution because such an approach includes rather large uncertainty. We have proposed three methods for varying the angles between the failure surface and the cross-section lines.

Theoretical background
"A spline is a flexible strip that is used to produce a smooth curve through a delegated set of points" (Bartels et al. 1987). The benefit of a spline is that if we increase the parameters, we can create a more accurate curve by increasing the degree of the polynomial. In the case of a landslide, if we have some spatially distributed borehole data (here to calculate the depth at a particular grid point, the spline curve used is a three-degree polynomial. If, in some cases, borehole data are available at some points, then the model with multiple splines segments must be used instead of the one segment spline used here. Mathematically, a spline is described as "a curve with a piecewise cubic polynomial (the degree can be increased) function whose first and second derivatives are continuous across the various curve sections" (Guéziec 1996). Here, we apply spline estimation to calculate the failure surface depth. In this section, we describe the methodology and the codes written to execute the process automatically. All the diagrams are illustrated using the example of the Kotropi landslide, Himanchal Pradesh, India. The overall procedure consists of these eight steps (Figs. 1, 2, and 3): S1. Draw the contour limits of the landslide failure. S2. Create the steepest slope joining the highest and lowest elevation points within contour limits.
S3. Take a point from the grid of the DEM. S4. Create a cross-sectional plane perpendicular to the steepest slope (S2) and pass-through grid point (S3).
S5. The slope of the failure surface must be fixed at least on points of intersections.
S6. Create a one segment spline curve. S7. Estimate the depth at the location grid point (S3) using spline value, topography, and coordinate changes.
S8. Iterate for all the parts of the grid inside contour limits (S1). The data inputs and the MATLAB codes and functions are described below to complete all the steps mentioned above. There are many functions required to compile the codes together, including -deg2utm (Palacios 2021), kml2struct (Slegers 2021), gridfit (D'Errico 2021), poly3n, and splin_param_new. The user needs to have the MATLAB Mapping Toolbox (MATLAB 2021) installed as an extension to run the code. This is compatible with organizing geographic data and allows for the interpolation, trimming, resampling, and transformation of coordinates.

Data required
The following data are required: 1. KML file of the contour lines tracing the boundary of the failure surface 2. DEM of the area in TIFF format covering the contour limits The imported data are converted into a suitable projection system based on the zone of the study area (WGS 1984 UTM projection) by using the function deg2utm. Additionally, the KML file needs to be converted into the structure using the function kml-2struct. These two data sources need to be overlapped to obtain the grid points lying inside and on the contour limits so that we can select the required points for calculating the Z values ( Fig. 1). Once the data are in grid form, we need to convert them into the columns of X, Y, and initial Z values so that we can apply the iteration for each point; for more details, read the script of the written code.

Create the steepest slope
We need to select the grid point with the highest and lowest elevations to draw the steepest slope of the landslide surface. A line is drawn passing through these points, which has been projected on the horizontal surface and named as elevation line. All the crosssection plans for each point are perpendicular to the elevation line.

Coordinate transformation
There are as many cross-sections as the number of grid points lying inside or on the contour line (Fig. 2). However, for ease of calculation, we need to rotate the axis so that its Y-axis becomes parallel to the line with the steepest slope/elevation line and the X-axis parallel to the cross-section lines. Therefore, we rotate the axis with angle θ (the angle formed between the initial Y-axis and the elevation line). Thus, each grid point has been transformed according to the new coordinate system.

Intersection of cross-sections with curves and sorting the data
The cross-section planes intersect with contour limits (Fig. 2), so we need to obtain the coordinates of the intersection to use these values as input during spline calculations. To achieve this, the data are sorted based on transformed Y values, and all the X coordinates are indexed according to new sorted columns of transformed X and Y coordinates. Using the polyxpoly function, the intersection points are calculated, and then these points are converted into the original coordinate system to calculate the Z values for each intersection point.

Calculating the elevation (Z values) for intersection points
For spline parameters, we need the X and Z values for each point, so we calculate the Z values using a function called gridfit. It requires the grid data as input and calculates the elevation at required points, which are also called as quarry points. For further processing, the calculated elevation data are combined into a single matrix.

Methods for calculating the variation of the angles
Once we determine the end coordinates, the slope angles at the endpoint of each cross-section need to be calculated. The result will be more accurate if we feed these values by measuring the dip and strike angles in the field, but it is impossible to measure for each cross-section. Because this is a hectic task and continuous measurement is not possible, we must interpolate this variation. Here, we propose three methods for varying the angles between the failure surface and the cross-section planes so that they can be used to calculate elevation using the spline method in the next step. These are as follows: 1. Constant average dip angle for the entire failure surface 2. Varying dip angle (nonlinear variation) 3. Varying dip angle (linear interpolation)

Method 1: constant average dip angle for the entire failure surface
In this case, it is assumed that the measured dip angle is constant for the entire failure surface. The angle formed between the crosssection and failure surface (Fig. 2) is also assumed to be constant for both profiles (the left profile is the portion lying left to the steepest slope line, and the right profile is the portion lying right to the steepest slope line (Fig. 6a, b)). This assumption can be made if the failure surface depicts geological structures, such as joints, faults, and bedding planes. In some cases, when structures (faults, bedding, etc.) are present, the development is highly similar to a classic geological cross-section (Chigira et al. 2013). In this method, there is no required presumption on the nature of the contour limits. We can easily feed the dip and strike angles for the right and left profiles and determine the calculated results. Only two inputs are required: left profile average dip angles and right profile average dip angles, which can be measured in the field.

Method 2: varying dip angle (nonlinear variation)
In this case, the dip and strike angles vary from top to bottom for the left and right profiles. The weight for the angle variation is given based on the distance of intersection points from the elevation line (Fig. 4B); see section Create the steepest slope. This is inversely proportional to the distance, i.e. IDW method. The higher the distance from the elevation line, the lower is the weight for the dip and strike angle variations. In Fig. 4, AB is the elevation line (horizontal projection of the steepest slope; see section Create the steepest slope), and OP and O'P' are the distances Fig. 4 Weightage for the variation of the dip in method 2, which is inversely proportional to distances OP and O'P' from this line to the intersection points. Generally, the change in curvature of the boundaries is higher at the upper and lower sides, where the distance from the elevation line is less. As in the case of an ellipse, the change in curvature is higher near the ends of the major axis, where the distance of a point on the ellipse is less from the major axis. Here, OP < O'P'; hence, the weight will be greater for variations in the strike and dip angle at intersection point P than at point P'. The same step will be true for the right profile; for more details, read the script of the written code. The code is adjustable depending on the user's requirement, and modifications can be made if required. It is better to use the 1st and 3rd methods. Through this method, we give more weight to the variation in dip angle to the top and bottom sides of the landslide because these are the areas where the distance from the elevation line to the intersection point will be less than that in the central part. Here, it should be clear that in the central part, the distances from the elevation line have less variation, so that we will have small changes in the dip variation in the centre for a large range. In general, if we are defining a circle regarding a point, the change in curvature will be higher for a smaller circle as compared to the larger circle, if we're moving on the circle. So, in this case, also the change in curvatures for the nearer points will be high as compared to the farthest. With this assumption, we are distributing the variation of dip and strike from the top profile to the bottom profile. The weight distribution based on the actual situation will be more accurate, but it will become a case-by-case calculation and not the generalized one. The codes are simple; if the user wants to modify the code, they can easily modify it to give the weightage based on the actual situation.
This method can be applied if the contour limits are not in a particular shape. There will be input for the left and right profiles, and each will require four parameters, which are as follows: 1. Top dip and bottom dip angle for the right profile 2. Top dip and bottom dip angle for left profile The values for the angles between the failure surface and each cross-section will vary nonlinearly between these extreme dip values, which are provided as inputs.

Method 3: varying dip angle (linear interpolation)
In this method, it is assumed that the left and right profiles are a part/arc of a large circle. Then, the dip angle is linearly distributed to the intersection points, creating a linear interpolation. This method can be applied if the contour limits are elliptical, circular, or consist of other convex shapes. Here, again, there are inputs required for the left and right profiles. We need two parameters for each profile: 1. Top dip and bottom dip angle for the right profile 2. Top dip and bottom dip angle for left profile For calculation purposes, the left and right profiles are defined with respect to the elevation line (horizontal projection of the steepest slope or line joining the highest and lowest elevation points), i.e. the portion left to the elevation line is the left profile, and the right side of it is the right profile. This is because, during calculation, it is difficult to distinguish the sign convention. If the elevation line is in the north direction (i.e. strike = 0°), then the situation is normal. If the strike of the elevation line is not in the north direction and has some angle (strike = θ°), then there may be some confusion regarding the measurements. To avoid this problem, it is better to define these values with respect to the elevation line. In the field, we can easily identify the peak point of a landslide, so for the top side (Fig. 5), for the top-right dip, we measure the dip value at a location in a range of 25-30° in a clockwise direction from the peak point, and for top-left dip measure, we measure the dip at a location in a range of 25-30° in a counter-clockwise direction from the peak (a rough approximation for locations). These will be the inputs for top dip angles. In Fig. 5, we can see the locations where we should measure the dip values of the failure surface (regions TL and TR for measurement on the top side, BL and BR for measurement on the bottom side).
The strike values 125 N and 235 N are default values used in the calculation, and we must change them if selecting different ranges for the measurement of dip. Similarly, for the bottom side, we can easily find the lowest elevation point in the field. Therefore, for the bottom right dip, only the dip value of the failure surface at a location 25-30° in the counter-clockwise direction from the lowest point is measured, and for the bottom left dip, the dip of the failure surface at a location 25-30° in the clockwise direction from the lowest elevation point is measured. Here, we are given a range of locations for measuring the dip angles so that we can create a linear distribution for the variation of these dip values. It is not necessary to measure the dip only in the range of the mentioned location; users can set their own range to measure the dip values, but they will have to feed the strike values of these ranges into the code script. Here, the strike of the abovementioned locations is set as default values in the script. Therefore, if other locations are used, these default values should also be changed; for more details, read the code script. In this case, the dip angle will linearly decrease  (Fig. 6a) and for the right profiles (Fig. 6b).

Calculation of Z values for the failure surface using splines
We calculate all four parameters required for the spline curves, which are the left endpoint and its first derivative (slope of the left profile measured along the cross-section line) similarly and the right endpoint and its first derivative (slope of the left profile measured along the cross-section line). The input data are stored in a Spm.parm structure, and the curve for each cross-section is obtained using the splin_param_new and poly3n functions. Once we have the curve, we can easily obtain the Z value for each grid point on a unique cross-section line. The corresponding spline curves are plotted in Fig. 7. The main concept of predicting the failure depth is based on a single spline curve method along multiple profiles. To calculate a value at a point, we need the curve/ Fig. 6 a Left profile and the intersection of the normals for the approximated arc. b Right profile and the intersection of the normal for the approximated arc Fig. 7 Spline curves for each cross-section. Each curve is used to calculate the Z value for a unique grid point lying within or on the contour line function or relation between dependent and independent variables. Here, in this case, we can get the required parameters to define the curve with the help of dip angles at the end of a cross-section in the form of slope at both ends and the coordinates of the endpoints of cross-sections. Then with the help of the dip angles and the end coordinate points of a particular cross-section, we are estimating the failure depth at a particular grid point. For each grid point, we are having this curve and estimated depth. Finally, we are plotting the data as a mesh of refined grids in 3D form; in this way, we are predicting the failure surface.
We know that the elevation for each grid point will be the same for any coordinate system because we are rotating in the horizontal plane, so the elevation for the original grid data and the transformed grid data will be the same. Therefore, we can convert all the X coordinates and Y coordinates into their original latitude and longitude settings to plot the data by creating the grid.

Plotting the 3D failure surface and stem
Once we have the Z values for the failure surface, we can create the refined grid and X_mesh, Y_mesh, and Z_mesh to plot the data in 3D and make it easy for visualization (Figs. 8 and 10). Similar grids are created for the original Z values for each grid point and a 3D surface for the failure surface, and the original surface is plotted together. Using plot3 and stem3, the point-wise Z values for the original and failure surfaces are also plotted to make them easier to interpret. The failure surface stem is a representation of Z values with a particular referenced X-Y plan. In this plot, Z is represented as a vector with a magnitude of Z values.

Volume calculation
When we import the standard data, e.g. in Tiff format, it is a Geo-graphicCellsReference with various information, such as cell size, raster size, number of columns and rows, latitudinal limits, longitudinal limits, and coordinate system type (Table 1). From this information, we can obtain the cell size of the grids multiplied by the required conversion factor to obtain ΔX and ΔY in the X and Y coordinates, respectively (suitable UTM projection): • Here, ΔX (in metres) = cell_extent_in_longitude (in the corresponding UTM projection) (Table 1

) • ΔY (in metres) = cell_extent_in_latitude (in the corresponding UTM projection)
Alternatively, we can divide the difference between the maximum and minimum longitude with the number of columns to obtain the cell size in the X direction, and similarly, the difference between the maximum and minimum latitudes can be divided by the number of rows to obtain the cell size in the Y direction. Fig. 8 Stem of the failure surface created after calculating the Z values using spline curves. Each grid point represents a unique Z value with respect to the X-Y plane Table 1 The values stored in geographic reference cells.
The code will automatically detect the projection system for the given area and will convert it into units of measurement. The multiplication of ΔX and ΔY gives the area of a particular grid cell. For the volume, we can multiply each grid area by the elevation difference. The elevation difference is the difference between the height of the original surface and the calculated depth of failure. The total sum for each grid point is the volume of the failure surface: where Z topo (i,j) and Z surf-fail (i,j) are the elevations of the topography and of the failure surface calculated using the spline; here, three situations arise. In the 1st case, if the mass of failure lies over the surface, we can use the post-landslide DEM file in our calculation of the volume because the mass is resting over the surface, so the DEM of the post-landslide will have the topography, and the failure topography is calculated using a spline. In the second case, the material is partially remaining on the failure surface, implying that the lower limit of the scar is buried. In that case, this lower limit must be extrapolated, and the method must be used only in the perimeter defined by the top of the deposit lying in the scar (the exposed surface must not include in the calculation) and the lower and lateral limits. In some cases, the calculation must be divided into two steps, the volume that moved out the scar on the pre-existing topography. The volume is obtained by subtracting the resulting DEM from the post-DEM. In the 3rd case, if the whole mass has been transported from the surface, then there is no need to calculate the failure surface, as it has already been exposed. For volume calculation, we can trace the boundary of landslides in Google Earth, and we need the DEM results of post-and pre-landslide failure events. From these two DEMs, we can generate the grid and obtain the elevation difference of failure depth and the original surface for each grid point stored in these two DEMs. Then, a summation of each cell volume can be made for the total volume calculation. If the pre-DEM does not exist, it can reconstruct by using the inverse SLBL (Jaboyedoff et al. 2020).

Application
To illustrate this method, we used the DEM and KML files for two known landslides. The model has been tested on the Kotropi landslide that occurred in Himachal Pradesh in India in August 2017 (Pradhan et al. 2017) and the La Frasse Landslide in Vaud, Switzerland (Tacher et al. 2005).
The Kotropi landslide (Fig. 9) has a latitude of N31°54′39.4″ and a longitude of E76°53′27.4″ and is located on the Indian Toposheet Number 53A/13, Himachal Pradesh, India. There is a trench on the other side of the highway, and a large amount of the failure mass of the landslide has been dumped in the trench because of the large displacement and favourable slope movement. As the road is busy and the main transport line, when the failure occurred, two-state government buses and a few other vehicles were hit and dumped into the trench.
According to the data released by the Himachal Government, at least 46 people died because of this disaster, and at least 500 m of the highway was also destroyed, which disrupted state transportation. In Fig. 10, we show a 3D model of the calculated failure surface and the original surface together. Method 1 was used with a right profile constant dip = 42° and a left profile with a constant dip = 59°. The volume has been calculated to be approximately 5.9 million m 3 . The calculated volume can be adjusted by changing the dip angles for the left and right profiles. The accuracy can be increased by measuring the dip and strike in the field and then taking the average. The third method can also be used by feeding the required inputs from the field. We have selected the Kotropi landslide as an example because currently, there is no available data on the volume of this landslide, demonstrating that the volume estimation by this method can be obtained very easily.
By analysing the time series data available on Google Earth Pro, we can easily visualise that the boundaries of failure limits were already demarcated; this is because it was not the first time this slide was activated. Two earlier incidents defined the boundaries of failure limits as per the available info, but no required actions were taken related to safety issues. According to Pradhan et al. (2017), along with rainfall, there was a pre-existing fault zone present in the area, a weathering rate of debris and soil layers, and the presence of low-strength materials. These were the primary causes of the slope instability. As mentioned above, this location is prone to landslides; people had reported issues since 1997 at approximately 20-year intervals.
The La Frasse Landslide has a latitude of N46°21′22.89″ and a longitude of E7°2′22.88″, located in the canton of Vaud (Switzerland). According to Tacher et al. (2005), it has a 2 km length and 500 m width, maximum depths of failure ranging between 40 and 80 m, and an active volume estimated by Tacher et al. (2005) of 41 to 42 million m 3 . The La Frasse Landslide is in a landslide-prone area. It is a compound translational slide following Hungr et al. (2014) classification. It is composed of tertiary flysch material. This is lying over flysch and limestone bedrock. The landslide is eroded at its toe by the Grande Eau River, which increases its activity. The area is covered by forest along with some habitation (Fig. 11).
In a similar manner, the grid point is generated from a DEM, and the contour limit is traced (Fig. 11a, b). Using the coordinates of intersection points and slope at these points, spline curves ( Fig. 12a, b) are generated, and then the required Z values are generated for each grid point. In Fig. 13, we show the 3D model of the calculated failure surface and the original surface together. In this case, method two is used because the contour limit is not The volume is calculated to be approximately 46.4 million m 3 . The calculated volume can be adjusted by changing the dip angles for the left and right profiles. This volume lies between the estimation of Tacher et al. (2005) of 42 million m 3 and the one calculated using the SLBL method (~ 49 million m 3 ) on the whole surface area of the landslide, which gives the active landslide volume for the most likely failure surface (Jaboyedoff and Derron 2013). This approach can be made more accurate by measuring the dip and strike in the field and feeding the input data.

Discussion and conclusion
These methods provide a rapid tool to estimate depths to a surface, which can be useful for slip surface location and borehole planning. This assessment indicates that using a spline calculation provides an accurate estimate of surface depth (Jaboyedoff and Derron 2013). This method can be made more powerful with additional field measurement parameters, such as whether the surface depth is known at a borehole location, and the surface can be constrained by using an estimate of the tangent and modifying the limits used to define the spline. Alternatively, if known parameters such as borehole data increase, we can use a multiple segment spline.
There was a limitation regarding the intersection points of the cross-section lines. If the cross-section lines intersect at more than  two points, then there will be an issue related to the number of grid points on that line. In the calculation, it is possible that we selected the wrong pair of points to define a spline, while the grid point lies in some other set of intersection points. This problem has been resolved by applying suitable conditional operations. However, it would be better if we could assign new grid points for each set of intersection points, but this needs further developments to be able to integrate several additional cross-sections automatically.
The implication of this model is straightforward. We can easily obtain the required data inputs with very little fieldwork, and the processing takes a few minutes to produce the result. Its advantages are that it allows the use of data sources available all over the globe that can be used based on a DEM and images such as tiff files, ASC files (with suitable conversion into a tiff format), KML, and other standard data inputs. The MATLAB script is available online on the GitHub platform. There are numerous parameters that can be manipulated accordingly to obtain better results. At the end, the model is storing the results in.tiff file format, which can be easily exported in ArcGIS and other software. The code is written on the MATLAB platform, but the results are provided in a standard format. Running the code will save the results (original grid, refined original grid, refined grids of the failure surface, and refined grid of the whole surface) in tiff format, which can be easily exported in ArcGIS and other software. Also, we will launch the model in the Python platform to avoid the work from being too site-specific. We can enrich the real scenario by imposing more constraints on the variation of the cross-section angle based on field measurements. This model is also helpful in predicting the area that can be affected by the mass of the failure surface, and we can make a predefined risk zone to avoid any casualties in the future.

Funding
Open access funding provided by University of Lausanne.

Available codes and examples
The codes, the example files and their descriptions are available online on the GitHub platform (https:// github. com/ Gauta m5673/ Failu re_ surfa ce_ estim ator).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.