An Idealized Landslide Failure Surface and Its Impacts on the Traveling Paths

Numerical scenario simulation may serve as an efficient and powerful tool for hazard assessment, but it often suffers from the lack of a definite failure surface before the occurrence of failure. In the present study, an idealized curved surface (ICS) is proposed for mimicking the sliding surface in the numerical simulation. Different from the conventional Sloping Local Base Level (SLBL) or Scoops3D method, this idealized surface consists of two curvatures, which are defined in the down-slope and cross-slope directions, respectively. Applying this idealized surface to 45 historical landslides of sliding type in southern Taiwan, two specific relations of geometry (length, curvature radius, and depth) in the down-slope and cross-slope directions are figured out. These specific relations simplify the complexity of constructing the idealized curved surface for areas prone to landslides for the sake of hazard assessments. That is, once the area with landslide susceptibility is identified and the associated released volume is given, the idealized sliding surface can be uniquely determined with the help of these geometric relations. The proposed method is integrated with a two-phase grain-fluid model and numerically validated against a historical large-scale landslide for investigating its feasibility. Although the idealized failure surface may deviate from that determined by the post-event investigation, it is interesting to note that the major discrepancy is mainly found at the first stage. The differences (both flow thickness and paths) reduce over time, and only minor discrepancies can be identified at the deposition stage. These findings reveal the weak co-relation between the geometry of the failure surface and the flow paths. It also indicates the feasibility of the ICS for predicting the flow paths by scenario investigation, especially when the exact sliding surface of a landslide is not available.


INTRODUCTION
Landslides, slope failures and the sequential mass movements caused by heavy or abnormal rainfalls often take place in mountainous areas and sometimes cause huge damages to property and people Lin et al., 2011;Iverson et al., 2015;Zhang et al., 2018). The volume of the released mass and its impact area generally determine the scale of the disaster. However, the estimation of the released volume and the prediction of the sequential flow paths of the moving mass for a landslide-prone area are still challenging tasks because of their high degree of uncertainty. Numerical simulations, founded on physics-based mechanical models (Kuo et al., 2009(Kuo et al., , 2013Mergili et al., 2020), discrete numerical approaches (Wu, 2010;Lo et al., 2011;Song et al., 2017;Yu et al., 2019) or statistics-based empirical models (Scheidegger, 1973;Manzella and Labiouse, 2013;Chen et al., 2014;Zhan et al., 2017), are common tools for risk assessment or for evaluating mitigation or remediation measures.
In terms of scenario investigation, the identification of the area prone to landslide and the estimation of the associated volume of moving mass are crucial for performing any numerical simulation. Owing to the rapid development of modern remote sensing techniques (LiDAR, SAR, InSAR, UAVSAR, etc.), even the tiniest topographic features (micro-topography) or structure patterns can be recognized, making it possible to identify the areas susceptible to landslides or to investigate the landslide dynamics (Lin et al., 2013(Lin et al., , 2014Stumpf et al., 2013;Delbridge et al., 2016). For example, Lin et al. (2014) suggested a methodological approach intended for detecting the large-scale landslides of potential in southern Taiwan by using LiDAR data. That is, with the hill-shade imagines given in different azimuth based on a high-resolution LiDAR-derived Digital Elevation Map (DEM) and with the help of aerial photos and geological field investigation, one may be able to identify the outline/boundaries of an unstable area that is prone to movement, i.e., identifying the scarps of potential landslides. Applying the technique of Uninhabited Aerial Vehicle Synthetic Aperture Radar (UAVSAR), Delbridge et al. (2016) developed a method for characterizing the 3-D surface deformation. Hence, in addition to the initiation mechanism, the estimation of the failure surface for a landslide-prone area has become an issue of high priority. Nevertheless, the geometry of the failure surface is tightly linked with the local geological structure and hydrological conditions. The complex composition of the material and the infiltration pattern of groundwater add to the complexity of the problem.
In general, the failure surface can be of different shapes, but the vertical section view along the flow direction can be roughly approximated by an arc of a circle as the most common case (Briaud, 2013). Hence, in additional to 2-D cases, a spherical-shaped failure surface is often used when investigating the performance of methods for 3-D slope stability analysis (see Lam and Fredlund, 1993;Huang and Tsai, 2000;Huang et al., 2002). As the high-resolution DEMs have become popular, a geometric interpretation of the failure surface based on a DEM, the Sloping Local Base Level (SLBL) method, was proposed for studying the failure mechanism as well as for a preliminary analysis prior to field investigation (Jaboyedoff et al., 2004(Jaboyedoff et al., , 2009). The SLBL-defined sliding surface possesses a constant second derivative, so that the surface exhibits a parabolic curve along the down-slope direction. Taken into account the curvature in the transverse direction, Reid et al. (2000Reid et al. ( , 2001 suggested a spherical failure surface together with extension of Bishop's (Bishop, 1955) limit-equilibrium analysis for quantitatively assessing the three-dimensional gravitational slope stability, explicitly for the stratovolcano edifice failure with large volumes without considering the internal/local structures. This searching scheme has been extended in the software Scoops3D for slope stability assessment in the framework of a digital landscape, where the local material properties in the subsurface and the hydrological conditions (groundwater configuration) can be taken into account Alvioli et al., 2018;Zhang and Wang, 2019).
The sliding surface defined by means of Scoops3D is spherical, so that there exists only one curvature. One may benefit from its simplicity for constructing the surface in the three-dimensional configuration. However, this single curvature approach is sometimes not sufficient for describing a plausible failure surface with respect to a complex topography. Recently, Kuo et al. (2020) proposed a smooth surface enclosing a given volume and with minimal surface area (volume-constrained smooth minimal surface) for mimicking the fracture interface of a deep-seated landslide, of which the boundary on the topographic surface fits that identified by the specific micro-topography through the high-resolution DEM and from a geological field survey. In Kuo et al. (2020), the smooth minimal surface is constructed under the assumption that the slope material is homogeneous and isotropic, i.e., the stratigraphy and geological structures are omitted. Nevertheless, the verification against 24 landslides triggered by excessive rainfall during the 2009 Typhoon Morakot confirms its applicability of providing acceptable predictions. Instead of the minimization of a smooth surface area, an idealized surface with two distinct curvatures in the downslope and cross-slope directions, respectively, is suggested for mimicking the sliding surface in the present study. This concept has been motivated by the investigation on the source scars of historical events, and it provides the freedom of matching the identified surface boundary of a landslide-prone area. However, it also requests additional constraints to define a unique failure surface in this regard.
In terms of the released mass volume, this method was used to assess the failure surface of 45 landslide events (source scars) that occurred in 2009 in southern Taiwan during Typhoon Morakot (Lin et al., 2013;Tai and Jan, 2019). The results indicate that there are two specific relations for the size (length and width) of the source scar area and the failure depth, in the down-slope and cross-slope directions, respectively. These specific relations may significantly simplify the complicated process of determining surface curvature. Together with the empirical formula expressing the relation between landslide area and volume, e.g., the regressed formula in the present study or an appropriate one, such as the volume-area relation proposed by Guzzetti et al. (2009), one is able to mimic the plausible failure surface, when the area susceptible to landsliding is given. We have applied the hereby proposed idealized curved surface (ICS) to the back-calculation of the 2009 Hsiaolin Landslide event for investigating its feasibility, where a two-phase model and the associated numerical code used in  are employed. We have then compared our results to the ones computed with the failure surface (FS) identified by post-event investigations. Not only the geometry of the surface is compared, but also the flow paths (flow-covered area) and the differences in flow thickness are examined at different time intervals, from movement initiation to deposition.
In section 2, we present the construction of the idealized curved surface, which consists of the formulation of the ICS (section 2.1), the regressed geometrical characteristics (section 2.2), and the procedure of determining the unique ICS (section 2.3). The training data of the 45 landslides are collected and listed in Appendix A, while their locations in Taiwan are depicted in Appendix B. The application of the method to the 2009 Hsiaolin Landslide is given in section 3, where the differences of the results computed with the ICS and FS are detailed. The data of the discrepancy evolutions are given in Appendix C. At last, the great application potential of the proposed ICS approach for stability assessment in variant scenarios is highlighted in the concluding remarks.

Idealized Curved Surface
In Geographic Information System (GIS), the terrain data are generally provided in a Cartesian coordinate system. In the present study, we set the xy-plane on the horizontal plane and the z-axis accounting for the elevation. Consider a specific surface, of which the curvature in the principal direction (x-axis) reads κ x and the curvature in the lateral direction (y-axis) is κ y . According to the differential geometry (Bronstein et al., 1996), the surface elevation z and these curvatures are related by where C x = 1 + (∂z/∂x) 2 −1/2 and C y = 1 + (∂z/∂y) 2 −1/2 , respectively.
As shown in Figure 1A, line p 1 p 2 stands for the main (downslope) axis of the failure area on the slope with p 1 and p 2 the upper and lower points, respectively. Line p 3 p 4 is perpendicular to the vertical section of p 1 p 2 and vertically above the midpoint of p 1 p 2 . The idealized curve ⌢ p 1 p 2 below p 1 p 2 is indicated by the thick line in Figure 1B, in which L is the length of p 1 p 2 , angle θ represents the inclination angle of p 1 p 2 to the horizon. In Figure 1B, p m denotes the midpoint of p 1 p 2 and d m accounts for the vertical distance (depth) from the midpoint of the idealized curve to p 1 p 2 . With the geometric relation, one can easily arrive at with R the curvature radius and θ arc the half expanding angle. Hence, it is clear that the central curve of the idealized surface ⌢ p 1 p 2 along the main axis is a function of the five parameters R, L, θ , θ arc , d m or R, L, θ , φ, d m , where φ = θ + θ arc is the slope angle of the surface along the main axis. In general, (L, θ ) is available once the failure area is identified and the down-slope direction is defined, but the angle θ arc is a non-linear function of d m and R. The curvature radius R and angle θ arc can be determined by means of iteration, once L, θ , d m are given.

Geometric Characteristic of the Idealized Surface
The proposed ICS is applied to mimic the failure surface of historical events. Comparing the difference between the DEMs, before and after the events, it is possible to identify the landslide area and the associated volume of the landslide mass. Frontiers in Earth Science | www.frontiersin.org and Jan, 2019). Among them, more than 470 victims have been identified in the Hsiaolin Landslide event. Being impressed by its large scale, the huge volume of released mass and the breach of the short-lived landslide dam, it has attracted particular attention of the government and engineering/scientific interests from research communities for natural hazard mitigation Kuo et al., 2011Kuo et al., , 2013Lo et al., 2011;Tsou et al., 2011;Lin et al., 2014), and the Hsiaolin case will serve as the application example in section 3 in the present study. Because these landslides were triggered during the typhoon, for each source scar only the area with failure depth more than 0.5 m are taken into account in investigating the geometric characteristics of the ICSs for the sake of isolating the erosion at the scar boundaries during the heavy rainfall. In the classification of landslide types, proposed by Hungr et al. (2014), they are mainly debris avalanches (and potentially debris flows in some highly channelized cases). For an overview, the 45 source scars are depicted in Supplementary Figures B1-B3, in which the red boxes in the inset panel indicate their corresponding locations in Taiwan. For each scar area, 10 points at the top and 10 points on the lowest margin are selected, yielding 100 combinations.
That is, we have produced 100 idealized surfaces for each source scar, of which the five parameters R, L, θ , θ arc , d m are computed, where d m is determined with respect to the post-event surface and the system (2) is solvable. In terms of the volume of computed landslide mass, the top 10 combinations are figured out, in which the mean values of the geometric parameters, such as the curvature radius R, depth d, and length L, are taken into account. Together with the failure area A L and the associated released volume V L , these mean values are summarized and listed in Supplementary Tables A1, A2 for reference. It is worth noting that, in both the main down-slope and transverse directions, there are specific relations among the mean curvature radius R, depth d, and length L of the 10 best-fitted curved surfaces. They can be obtained by regression and given by Frontiers in Earth Science | www.frontiersin.org regressed specific relations are illustrated in Figures 2A,B, where panel (A) shows the data in the down-slope direction, and panel (B) is for the data in the cross-slope direction. In Figures 2A,B,D, the solid red circles and green squares, respectively, stand for the landslides taking place in the catchments of DF081 and DF054, the magenta triangle markers represent the source scars for the Hsiaolin event, and the blue diamond markers denote the ones in the mountainous areas of the Kao-Ping River basin. Although the two regressed lines are close to each other in the plot, it is inappropriate to merge them. With the help of Figure 1B, one can find that formula (3) is close to the ideal curve with inclination angle θ = 32 • , and relation (4) approximately follows the curve with θ = 10 • (cf. Figure 2C). Hence, relation (3) indicates a representative inclination angle (around 32 • ) of p 1 p 2 in the down-slope direction, and relation (4) reflects the fact that the transverse slope (cross-slope) of the source scar is slightly tilted in general.
The specific relations, (3) and (4), may simplify the complicated iteration process of solving (2) without knowing the magnitude of θ . In using them, one should pay special attention to the definitions of d 1 and d 2 in (3) and (4). As shown in Figure 1A, p 1 p 2 and p 3 p 4 may be on different levels so that d 1 and d 2 would have distinct values with respect to a given failure depth. For an identified area of a slope failure of potential, p 1 is chosen to be at the top of the area with p 2 at the lowest point, where the line p 3 p 4 is perpendicular to the vertical section of p 1 p 2 and vertically above the midpoint p m , see Figure 1. With respect to a given failure depth/position below the midpoint p m (cf. Figure 1B), the values of d 1 and d 2 can be easily determined. Because L 1 and L 2 are given by the lengths of p 1 p 2 and p 3 p 4 , R 1 and R 2 can, therefore, be determined with the help of (3) and (4), and the idealized failure surface is then uniquely determined.

Construction of the ICS for a Landslide-Prone Area
Until today, the ways to determine the failure depth are still open to debate. Conventionally, the failure depth is determined by the analysis of the safety factor with respect to a pre-defined failure surface (Bishop, 1955;Bishop and Morgenstern, 1960;Zheng, 2012;Das and Sivakugan, 2016). In the present study, we consider an alternative method, i.e. through the empirical relation between the source area of a landslide and its associated released volume. Since the 60's, the volume-area relation has been investigated based on individual landslides of variant types in different areas (e.g., Simonett, 1967;Rice et al., 1969;Imaizumi and Sidle, 2007;Parker et al., 2011). Guzzetti et al. (2009) analyzed 677 landslides of sliding type in Italy and suggested the representative empirical formula, with R 2 = 0.9707, where A L (in m 2 ) stands for the area and V L (in m 3 ) denotes the associated landslide volume, see the blue dash-dotted line in Figure 2D. In the present study, a similar relation is proposed (see the black solid line in Figure 2D), appropriate for the analyzed 45 landslide events (source scars) in southern Taiwan. As shown in Figure 2D, these two formulas, (5) and (6), are very close, while their reference cases lie in distinct areas. Although it has been elaborated in Larsen et al. (2010), in which 4, 231 cases were collected, that the exponent ranges from 1.1 to 1.3 for soil-based landslides and 1.3 to 1.6 for ones involving the failure of bed-rock, the similarity between (5) and (6) reveals more or less their applicability in the corresponding local area. With the help of an appropriate volume-area relation, such as (6) or (5), one may estimate the landslide volume for constructing the associated ICS, once the area that is prone to landslide is identified. As elaborated, each failure depth (d m in Figure 1B) uniquely defines its corresponding ICS for the identified area. Hence, the failure depth can be determined according to the estimated volume, and the method of iteration can be employed once the landslide volume is given. Figure 2D shows the volume of released mass against the failure area for the analyzed 45 landslide events (source scars), where the blue dash-dotted line indicates formula (5) and the black solid line represents relation (6). The sound agreement indicates that the empirical relation (6) or (5) together with the specific relations, (3) and (4), can be employed for determining the failure volume as well as constructing the idealized curved surface.

APPLICATION TO HSIAOLIN LANDSLIDE IN TAIWAN, 2009
The proposed idealized failure surface is then applied to retracing (back-calculating) the movement of a large-scale landslide event, the Hsiaolin Landslide that took place in Taiwan in 2009. The numerical simulation is performed and compared against the results computed with the post-event failure surface. In the computation, the two-phase debris model for non-trivial topography  is employed, where all the parameters are set identical to the ones used in . The computational domain was 3, 700 × 2, 210 m with mesh size x = y = 10 m, and three source scars (HL-1, HL-2, and HL-3, cf. Supplementary Figure B1C) were taken into account in the computation. All the sliding masses were released simultaneously for simplicity (Kuo et al., , 2013. Although this scenario (simultaneous release) might be highly unlikely in reality, it is a compromise for the lack of definite evidence of the release sequence, and the influence is rather minor. For the sake of brevity, we are not showing hereunder the model equations.
For details on the model equations and applied parameters, the readers may refer to .
As the idealized curved surface is a smooth surface, the calculated geometry of the initially released mass is certainly different from that calculated with respect to the post-event failure surface. In Figure 3, the released landslide mass calculated by considering the idealized curved surface (ICS) is shown in panel (A), whilst the initial landslide mass calculated with respect to the post-event failure surface (FS) is given in panel (B). The difference between these two masses is significant: the thickness discrepancy ranges from −42.58 to 59.34 m (cf. panel C or  Table C1). Moreover, the resultant scar areas are not identical and a deviation around 23% has been identified. Figure 3D illustrates the discrepancy of areas covered by the idealized curved surface (ICS) and the post-event failure surface (FS): the area in blue-green indicates the one covered both by ICS and FS; the area only occupied by the ICS is marked in golden yellow; and the blue color represents the scar area not covered by the ICS (i.e., FS only). In order to point out the size of the landslide, the Hsiaolin village is shown for scale, highlighted by the red line in Figure 3C. Figure 4 illustrates the evolution of the thickness difference computed with the ICS and the FS, where the data of flow thickness <1 cm are isolated. Together with Figure 3C, the left panels of Figure 4 illustrate the differences in flow thickness at various time levels. The discrepancy of the initial thickness is significant (Figure 3C), but the difference reduces as time passes (Figures 4A,C,E). The range of difference drops from [−43.4 m, 61.9 m] at t = 0 s to [−16.7 m, 11.7 m] at t = 181.83 s (cf. the blue shade in Figure 5), when the flowing body is nearly at the state of rest. Figures 4B,D,F depict the flow paths at different moments in time. It is found that the discrepancy dramatically decreases, especially for the flow paths (flow covered area). Figure 5 summarizes the evolutions of the volume difference (green line), the discrepancy of the covered area (red line) and the depth difference (blue shaded) from t = 0 s to t = 181.83 s. As shown in Figure 5, the variation of the volume difference (green line) slightly increases, from 2.83 to 3.74%, and this increase is suspected to be caused by the outflow from the computational domain. The discrepancy of the flow paths (flow-covered area, red line) is 26.82% at t = 0.0 s and reduces to 9.48% at t = 181.83 s. The notable variation takes place at the first stage (t < 80 s), and it is supposed to be caused by the difference between the ICS and FS. For t > 75 s, the discrepancy monotonically decreases and tends to be around 9.5% at the deposition stage. Regarding the thickness difference, a significant variation of the depth difference (blue shaded) is also noted when t < 42.4 s, while the range asymptotically reduces to [−16.7, 11.7 m] at t = 181.83 s. For detailed information, the data of the discrepancy evolutions are collected in Supplementary Table C1.  Figure 6A is an orthophoto (satellite image) taken ca. 7 months after the Hsiaolin Landslide event, in which the bare zones approximately indicate the areas covered by the moving mass (flow paths, source scars, and deposit areas). In Figure 6B, the computed flow paths during the period of computation (from t = 0.0 s to t = 181.83 s) have been marked with respect to the satellite image, where the idealized curved surface (ICS) is employed in the computation and the area with flow thickness <1 cm is isolated. The computed flow paths are composed of 61 sets of results, i.e., with an increment of 3.0305 s, which are in sound agreement with the bare areas in the satellite image. As shown in Supplementary Table C1, the discrepancy of the flow paths between the results computed with the ICS and FS is 26.82% at the stage of initiation, causing the overflows at the northern and southern flanks of the source scar. This discrepancy reduces, as time marches, to ca. 9.5% at the stage of rest. These findings  concerning the reductions of the discrepancy in flow paths and of the range of thickness difference indicate the feasibility of employing the ICS for representing the FS by retracing or predicting the plausible traveling paths, especially before the failure takes place and the slope-failure area of potential has been identified.

CONCLUDING REMARKS
In the present study, an idealized curved surface (ICS) is proposed to mimic the failure surface of a landslide of sliding type for the sake of performing numerical simulations. Through the investigation on historical events by means of the proposed ICS, specific relations are suggested for R/L and d/L in the down-slope and cross-slope directions. Another interesting finding is that, for the 45 source scars in southern Taiwan, the specific relation between landslide volume and area is very close to that suggested by Guzzetti et al. (2009) for landslides in Italy. These relations provide useful information for constructing the ICS. Since each failure depth d is associated with one unique ICS as well as the volume to be released, one may determine the corresponding ICS through iteration, once the landslide volume is given. For a landslide-prone area whose boundary is identified, the plausible released volume can be estimated by the empirical volumearea formula, either (5) or (6), so that the idealized sliding surface can be determined by iteration. Based on the constructed ICS, the associated plausible released volume is determined and numerical simulations can therefore be performed for the sake of hazard assessment or evaluating remediation measures.
The feasibility of the employment of the ICS is investigated through the back-calculation of a large-scale landslide with a twophase grain-fluid model . The ICS is constructed with respect to the source scar and released volume given in . Because the ICS is smooth but the FS is rather rugged, significant deviations of the flow thickness and flow paths are observed at the first stage. However, it is interesting that the discrepancy decreases as time passes and the deviations for the final shapes of the deposit are rather minor. The difference of flow-covered area drops from 26.82% (initial condition) to 9.54% at the deposition stage. In the computational period, the range of the depth difference has shrunk 73%, i.e., from [−43.4, 61.9 m] to [−16.7, 11.7 m]. In comparison with the post-event orthophoto, the results computed with the ICS allow to approximately retrace the flow paths as indicated in the photograph. This finding reveals that the flowing paths of a landslide are dominated by the topographic condition, such as channelization, and less sensitive to the geometry of the failure surface, especially after a long traveling distance. That is, in case the main part of the moving mass flows into a channelized topography at the first stage, the geometry of the failure surface does not play a crucial role.
It is worth pointing out that the proposed ICS does not take into account the spatial geological and hydrological variations, weathering effects, or the complex and inhomogeneous (nonisotropic) composition of the site. The current ICS is a compromise, which is not suggested for replacing those determined by the stability analysis. It only provides a plausible failure surface for a quick hazard assessment under the circumstance that no detailed local geological and hydrological condition is available. On the other hand, the current proposed ICS can provide the candidate surface for stability analysis as what has been done by Scoops3D , so that the hydrological condition or the local geological structure can be taken into account. That is, the ICS provides the possibility of investigating the potential area of slope failure by integrating the above mentioned hydrological condition as well as the local geological structure, which are generally of high uncertainty. With respect to the ICS, variant hydrological conditions (such as the ground-water level, distribution pattern, or seepage directions, etc.) can be applied for evaluating the associated safety factors, where failure is supposed to take place at the lowest value of the safety factor. That is to say, it enables us to retrace the plausible hydrological condition at the initiation stage for the historical landslide events. Moreover, variant scenario investigations can be performed for the purpose of hazard assessment or evaluating the countermeasure of disaster mitigation. With these, a high potential of the proposed ICS in engineering applications can be expected.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
Y-CT: provide the idea, design the work, and construct the MS. C-YK, R-FC, and C-WL: concept/opinion exchange and ensure the contents. C-JK, K-DL, and Y-CW: perform the data analysis and numerical simulation. All authors contributed to the article and approved the submitted version.