Application Assessments of Using Scarp Boundary-Fitted, Volume Constrained, Smooth Minimal Surfaces as Failure Interfaces of Deep-Seated Landslides

More than 9,000 potential deep-seated landslide sites in the mountain ranges of Taiwan have been identified by a series of renewed governmental hazard mitigation initiatives after the 2009 Morakot typhoon. Among these sites, 186 sites have protection targets where thorough mitigation strategies are to be implemented. One of the important tasks in the hazard mitigation initiative is to estimate the volume, failure interface and related quantities of each landslide site. In addition, with this number of sites, an automated tool is needed to generate predictions at low operational costs. We propose to use volume-constrained smooth minimal surfaces to approximate the landslide failure interfaces. A volume-constrained smooth minimal surface in the current context is defined as a differentiable surface that encloses a given landslide volume with the minimal surface area. Although the stratigraphy and geological structures are omitted, the smooth minimal surface method is verified with 24 known landslides and is shown to be able to generate acceptable, approximated failure interfaces. A collection of assessment indices is employed to measure the fitness of the predictions. Finally, the prediction fitness vs. the landslide scarp geometry is investigated.


INTRODUCTION
Deep-seated landslides pose severe threats to human lives and property. Typhoon Morakot struck Taiwan in 2009 which brought approximately 2,500 mm of precipitation in 4 days to the southern parts of the island and triggered numerous landslides, debris flows and vast flooded areas. This catastrophic extreme climatic event caused more than 22,000 landslides. Among these mountainous hazards, more than 320 landslides were found in scarp areas larger than 10 ha (Lin et al., 2011). These large-volume landslides often lead to composite casualties. For example, the Hsiaolin landslide, after sweeping through the village in its course, formed a short-lived blockage dam, and the continuous inflow triggered follow-up dambreak debris flows Li et al., 2011;Tsou et al., 2011). Similar rainfall triggered large scale deep-seated landslide examples and related research include those reported by Cardinali et al. (2002), Roering et al. (2005), Baroň et al. (2011), Chigira (2011), Xu et al. (2015), Vallet et al. (2015), and Lee et al. (2018).
Having learned from the Morakot typhoon landslides, Taiwan's government authorities officially defined the deepseated landslides for administrative purposes according to their geometric measurements: volumes larger than 1 × 10 6 m 3 , areas larger than 100 ha or depths deeper than 10 m (Lin et al., 2011;Chen et al., 2017). Using geometric measurements to classify the deep-seated landslides simplifies the administrative process and has been suggested in the literature (Roering et al., 2005;Lo, 2017). There are in fact other definitions of deep-seated landslides in different research contexts. For example, geologists may refer the term to slow moving large-scale landslides with failure surfaces occurring deep in a rock bed, or geotechnical engineers may refer it to landslides with failure surfaces below the underground water table. Nevertheless, as literally suggested, deep-seated landslides are usually associated with large slide volumes such that, in the present paper, we adopt the simple geometric definition for this type of landslides. Because of the large scale of the potential deep-seated landslides, the spatial geological variations, weathering effects and orogenic activities contribute to forming these sites. During long-term evolutionary processes, topographic features, such as crowns, bulges, and trenches, develop, and currently, these features can be detected by modern remote sensing techniques (Varnes, 1978;Chigira and Kiho, 1994;Chigira et al., 2013;Crosta et al., 2013;Lo, 2017). In a series of hazard mitigation initiatives implemented by the government starting from 2010 (Central Geological Survey, 2010) and continuing into the present (Soil Water Conserv. Bureau, 2017;Soil Water Conserv. Bureau, 2018), in which other sequential projects can be checked, airborne light detection and ranging (LiDAR), various satellite synthetic-aperture radar (SAR), and field investigations have been combined to survey the topographic surface activities (Lin et al., , 2014Tseng et al., 2013;Chen et al., 2015;Wu et al., 2017). Through these efforts, the scarp boundaries of the deep-seated landslides are identified according to the surface features. More than 9,000 potential deepseated landslide sites have been found, and among them, 186 sites have protection targets (Figure 1). In the figure, a close-up view of ten sites in Ren'-ai Township, Nantou County, is shown. Along this direction, assessments of the landslide volume and influence area are to follow hazard mitigation planning.
In traditional geological engineering approaches, the factor of safety of the potential landslide site will then be calculated by using slope stability tools. The most well-known methods are the those based on the limit equilibrium concept, published by Bishop (1955), Mogenstern and Price (1965), and Janbu (1973). In these methods, the slope body above the prescribed failure interface is discretized into a number of vertical slices (free bodies), and a system of algebraic equations is derived according to the force and moment equilibrium, with an additional proposition of internal forces as closure conditions. Because the failure interface is a presumed input parameter in these limit equilibrium methods, there are a few empirical strategies for proposing the interfaces. These strategies include circular, piecewise linear, and other methods, and the choice is made according to the geological condition of the investigated site. The circular failure interface is the most widely employed because this shape is easily parameterizable and agrees with a large portion of common observations. Consequently, automatic iterative search schemes for the least stable interface are implemented in many analysis computation tools (e.g., Siegel, 1978). Once the least stable failure interface is obtained, the landslide depth and 2D volume can be estimated.
Although the analysis is convenient and straightforward, 2D slope stability can be performed on only one or a few heuristically selected representative profiles for each landslide site. Because of this limitation, the analysis may not provide sufficient information for deep-seated landslides in estimating landslide-influenced areas. For this purpose, the 3D scarp depth distribution of the slide mass, for example, is one of the essential quantities required to estimate the runout and spread of the landslide, and it is particularly important for rapid large-scale avalanches, whose motion can be largely influenced by the topographic conditions of the terrain Luca et al., 2016;Tai et al., 2019), and the references therein.
The concept of the factor of safety has long been the key concept in slope stabilization plans. Since the wide deployment of the aforementioned 2D methods, there have also been attempts to extend the slope stability analysis to general 3D terrains. Similar to the 2D limit equilibrium formulation, most of these methods require a priori failure surface; the derived system of equations is a statically indeterminate system, and additional conditions are needed to solve it. For example, in a study by Leshchinsky and Huang (1992a) and Leshchinsky and Huang (1992b), the landslide mass is assumed symmetric and variational optimization is adopted for searching the minimum factor of safety. Ignoring the shear stresses in the internal landslide mass (Ugai, 1985) showed that the failure profile is circular along the sliding direction. A two-directional moment equilibrium method was proposed by Huang and Tsai (2000) and Huang et al. (2002), which can directly determine the sliding direction, instead of presuming a sliding direction. Further improvements are made by including the complete momentum and moment equilibrium conditions (Zheng, 2012;Jiang and Zhou, 2018), in which examples of practical 3D applications are demonstrated. In the group of slope stability analysis methods, Hungr noted that the unbalanced force in the transverse direction of the potential landslide mass is responsible for the errors in the factor of safety (Hungr, 1987;Hungr et al., 1989). In contrast to the limit equilibrium concept, elasto-plastic approach is a full determinant method (Griffiths and Lane, 1999). The adoption of the full determinant method requires more complete underground material parameters and stratigraphic details, such that fitting the failure surface obtained by this type of method to the identified landslide scarp outlines may involve a series of procedural justifications of the material parameters and geological conditions. It has a different mechanical complexity beyond the present consideration and is not pursued further here.  '-ai Township, Nantou County (Central Geological Survey, 2010). The dots in the scarp areas represent the line-of-sight deformation obtained by using the temporarily coherent point interferometry SAR (TCPInSAR) technique . Note that the length scale varies according to the 3D perspective.
Because the limit equilibrium approach for slope stability analysis yields a statically indeterminant system, additional empirical conditions, the failure surface and internal forces are needed to calculate the factor of safety. Unlike the 2D cases, there is a lack of systematic, convenient and commonly applied ways to form the 3D failure interfaces that fit the observed landslide scarps. However, simple spherical-shaped failure interfaces are often used as illustrative applications in developing 3D slope stability methods (Xing, 1988;Lam and Fredlund, 1993;Huang et al., 2002). Treating spherical surface sections as failure interface elements, a searching scheme has been proposed to construct the regional distribution of the factor of safety (Reid et al., 2000(Reid et al., , 2001. Although the method does not focus on finding precise 3D matches of landslide scarps for individual sites, it has become increasingly popular in creating landslide susceptibility maps (Reid et al., 2015;Tran et al., 2018;Zhang and Wang, 2019).
In the examples of the identified potential deep-seated landslide sites in Figure 1, we can see that the scarp areas are identified by the closed polygons and the line-of-sight deformation obtained by using the temporarily coherent point interferometry SAR (TCPInSAR) technique somewhat indicates the landslide activity . In addition due to the large area that can cover topographic heterogeneous landscapes, each landslide site may contain multiple secondary failure structures, e.g., the cracks and minor scarps in Figure 1 (Central Geological Survey, 2010). In landslide hazard mitigation plans with these types of information, a tool that it is able to generate 3D failure interfaces, conform with the surface scarps and provide related depth information for the landslide mass is required. Such a tool also provides great benefits for landslide sites, of which multiple secondary hazard scenarios are to be studied and managed. Furthermore, with the massive number of 186 potential landslide sites, the tool should be automated so that it can be operated at low costs.
As a tradeoff, the tool is aimed at providing predictions with satisfactory accuracy at small operational efforts but not necessarily with high degrees of precision. In this regard, we propose a simple method for computing the failure interfaces of deep-seated landslides, that is based on the minimization of a smooth surface that encloses a given landslide volume with the specified scarp boundary. The concept arises from the observations that (1), the scarp boundaries of potential deep-seated landslide sites often develop over time and become observable on the surface, for example crown fissures and flank scarps, and (2), landslide volume-area relations can be used to determine the volume from the scarp area. As a result, the average landslide depth is well-constrained, and the 3D failure interface and slide volume can be obtained.
In the proposed method, the geological settings and hydrogeological conditions of each landslide site are neglected. Geological settings include the rock texture, lithological stratas, rock cleavage or joint orientations, etc. These factors affect the intrinsic structure of the failure surfaces and landslide failure pattern. On the other hand, the hydrogeolocal conditions include the rainfall and underground water hydrology, which alter the force balance condition of the slide mass and trigger the landslide. Both geological settings and hydrogeological conditions are site specific details which lead to the deviations of the actual landslide failure surfaces from the predicted. However, to construct these details even for a single site may require extensive resources and investigations such that full coverage surveys and hazard monitoring for a population of landslides become virtually impossible. Under these circumstances, we defer the adoptation of any site specific detail to the present method.
Hence, it is an important task to establish the accuracy baseline of the method in the present paper, such that the paper is organized as follows: A brief introduction to the landslides used in the assessments is presented in section 2; The smooth minimal surface method and statistical assessment indices are described in sections 3 and 4, respectively; the applications of the method to two landslides and two conceptual examples are detailed in section 5 and Appendix 1, respectively, which includes the comparison of the predicted landslide failure interfaces to the actual failure interfaces and calculations of the assessment indices. Then, the application set is extended to include a total of 24 landslides, and their assessment indices are tabulated in section 6. With this amount of data, the accuracy bounds can be inferred.

LANDSLIDES IN THE STUDY AREA
In total, 24 deep-seated landslides were selected in the present study (Lin et al., 2011). They were all triggered by the excessive rainfall of the Morakot typhoon and are distributed in three different catchment areas in the southwest mountainous range of Taiwan. The first group, containing 11 sites, is located in the Cishan River catchment area, Jiaxian District, Kaohsiung City, as shown in Figure 2. All these landslides have areas larger than 10 ha, maximum depths over 10 m and, according to the definition from the Taiwanese government, are classified as deepseated landslides. Among them, the sites labeled with HLIN prefixes are the scarps associated with the Hsiaolin landslides Tsou et al., 2011;Tai et al., 2019). The landslides are defined by using the 2005 and 2010 LiDAR 1 m resolution digital elevation maps (DEMs). The second and third groups of landslides are identified by the DF054 and DF081 prefixes, respectively (Figure 3). DF054 is in the Longjiao River catchment, Dapu Township, Chiayi County, and DF081 is in the Laonung River catchment, Maolin District, Kaohsiung City.
The geological structure and stratigraphy of these catchment areas are very briefly reviewed here as background information. Because of space limits, the geological maps of the three catchment areas are relegated to the Supplementary Materials of the paper and are made by referencing (Fei and Chen, 2013). To summarize, the landslide sites of the first group are distributed in three types of surface strata in the Cishan catchment area: Hunghuatzu Formation, Changchihkeng Formation, and Tangenshan Sandstone. These units are arranged chronologically, with the oldest formed in the late Miocene. The landslide sites of DF054 are in the Tangenshan Sandstone and Ailiaochiao Formation (early Pliocene). They consist of sedimentary rocks and are mainly composed of sandstones and shales, in which marine microfossils commonly occur. The third group, the DF081 landslides, is within Chaochou Formation, which is middle Miocene in age and occasionally interlaces with quartzite.
The landslide volumes (V) and areas ( p ) 1 of the landslides are tabulated in Table 1 and the data are plotted in Figure 4A. In the figure, the regression model is drawn and compared to the well-known equation from Guzzetti et al. (2009) and Klar et al. (2011). It is found that the present 24 landslides agree excellently with the fitted regression line. With these two quantities, we can calculate the equivalent radius: R = p /π 2 . The characteristic length scale 2R, which is the equivalent diameter, will be used as the normalization length factor for statistical quantities varying in the horizontal or slopewise directions, such as the distance of the gravity centers between the actual and predicted landslide volume. In addition, two additional non-dimensional geometric parameters, the roundness r r and the sphericity r s , can be defined with these measurements and the scarp boundaries.
We use standard mathematical definitions to define the two parameters. The roundness is defined as the ratio between the radius of the maximum inscribed circle of the scarp and that of the minimum circumscribed circle. Its values range between 0 and 1, with the two extreme values corresponding to an infinitely thin and perfectly circular shapes, respectively. For ellipses, the roundness reduces to the aspect ratio between the width and length. With this analogy, we can associate the roundness to the landslide aspect ratio, which commonly appears in the landslide literature. By adopting this general definition, calculation ambiguities for scarps with complicated shapes can be avoided. On the other hand, the sphericity r s is defined as the ratio of the surface area of a sphere of a volume V to the total surface area of the landslide mass S 0 + S b (the sum of the free surface S 0 and the failure interface areas S b ), i.e., r s = π 1/3 (6V) 2/3 /(S 0 + S b ). Its value also ranges from 0, an infinitely thin volume, to 1, a perfect sphere. Because the landslide thickness is usually much smaller than the other spanwise dimensions, the sphericity becomes a factor involving the landslide thickness and the slope 3 . These two parameters in the above definitions have long been applied in various landslide studies but have different terminologies, such as the width-length ratio and the depth-length ratio (Taylor et al., 2018).
There appears to be a relation between the roundness and sphericity with the present landslide inventory ( Figure 4B). The data are somewhat evenly distributed in the range of the roundness and sphericity; i.e., no favorable clustering spots of data are found. Instead, the sphericity generally increases with increasing roundness. Among the sites, FID12 is likely a statistical outlier because of the distinctive gap between it and the other data. Inspecting the site in Figure 2, we find that the scarp boundary of FID12 appears to have a peculiar shape, colloquially a dumbbell shape, with landslide mass biasedly distributed at both end lobes but its relative roundishness reduces the sphericity. The site is likely composed of two distinguished landslides instead of an integrated one. We, however, do not perform further manipulations on this site other than simply exclude it from the accuracy assessments in section 6.
The other important landslide geometric parameters, such as the mean and maximum depths are calculated and listed in Table 2 for facilitating the later accuracy assessments with the failure interface predictions. The related discussion will resume FIGURE 2 | Landslides in the Chishan River catchment, Jiaxian District, Kaohsiung City (Central Geological Survey, 2010). Negative values of the color legends indicate landslide scarp depths and the positive values indicate the deposit thickness in meters. The HLIN prefixes mark the Hsiaolin landslide and its associated landslides. The sites outlined in red solid lines are the deep-seated landslides investigated in the present study. The TWD97 coordinate system is used.
in section 6, after the description of the method of the smooth minimal surface, assessment indices and casewise applications.

SMOOTH MINIMAL SURFACE
Unlike the full mechanical approach involving fracturing or plasticity, the failure interface is not a determined result but a prescribed prerequisite. Based on observations, circular (2D) or spherical (3D) shaped failure surface profiles are usually applied for slopes with homogeneous and isotropic materials. To fit the application scenario, we relax the surface to a smooth minimal surface. In our application, the minimal smooth surface is obtained by giving an assumed failure volume V with the constraint that its boundary on the free surface matches that from geological field investigations. This type of surface is chosen because it has a close relation to the spherical failure surfaces in 3D slope slability analysis and can degenerate to the spherical surfaces for some special cases, cf. Appendix 1. Therefore it may also inherit the same property that the method fits better for slopes with isotropic and homogenous materials, following the same reasons of the traditional slope stability analysis.
Determining the smooth minimal surface is an optimization process. Let (x, y, z) be the Cartesian coordinates with z vertically pointing upwards, and let the failure interface z b ≡ z b (x, y) be a smooth differentiable surface, as shown in Figure 5. Mathematically, the area of the surface can be calculated as such that the minimum surface can be acquired by minimizing (1) by varying z b ; i.e., Frontiers in Earth Science | www.frontiersin.org     subjected to the constraints and In the equations, ∇ is the 2D gradient operator (∂/∂x, ∂/∂y); p is the xy−projection area of the failure scarp area; and Ŵ is the boundary of p . The first constraint states that the boundary elevation of the failure interface is equal to the elevation of the free surface, z 0 ≡ z 0 (x, y); the scarp surface is below the free surface; and the second FIGURE 5 | Conceptual sketch, symbols, and coordinate system. The landslide site FID5, defined in Figure 2 is used as an example. Symbols z 0 ≡ z 0 (x, y) and z b ≡ z b (x, y) are the free surface and underground failure interfaces, respectively. p and Ŵ are the projection of the scarp area and outline on the xy plane, respectively. The scarp failure interfaces and depth are largely exaggerated for clarity in this sketch, and the color map conceptually indicates the scarp depth.
constraint states that the failure volume is the prescribed volume V.
The scheme is implemented in the open source finite element platform FreeFEM++, ver. 3.47 (or above) (Hecht, 2012), and the integrated Interior Point OPTimizer (IpOpt, Wächter and Biegler, 2006) optimization scheme in FreeFEM++ is applied to find the numerical solution. In the scheme, the Jacobian and Hessian matrices of the target function (1) are used to accelerate the numerical convergence: and The above two expressions are written with the help of the notation of the Fréchet derivative for facilitating the numerical scheme implementation (Vergez et al., 2016). Similarly, the Jacobian matrix of constraint (4) is Equations (5), (6), and (7) are evaluated at discrete points in each mesh element. Supporting external C++ functions are coded to handle the DEM, scarp boundary outline inputs and outputs to other software for geographic information systems (GIS), such as ArcGIS, qGIS, and GRASS. Initially, the outline is uniformly discretized with a prescribed number into a set of linear edges, and the mesh of the scarp area is generated. Then, the optimization (2) for the minimal surface is executed.
Depending on the precision requirement of the mesh, multiple passes of mesh adaptation can be performed based on the calculated failure interface. As seen in the formulation and calculation principle, the smooth minimal surface can be equivalently replaced with other mathematical surfaces, provided that they are phenomenologically reasonable and their descriptive expressions are available. It is also in principle possible to design transformations that can further manipulate the predicted failure surfaces to incorporate with site specific geological settings or hydrogeological conditions. Nevertheless, because the method is a phenomenological approach, the accuracy of the failure surface prediction has to be verified with actual landslide data. In what follows, we use the current smooth minimal surface to proceed the application example and accuracy assessment. The assessment indices and procedure are generic and applicable to future alternative failure surface prediction or transformation schemes.

ASSESSMENT INDICES
The proposed prediction method for the failure interface is purely mathematical. Intuitively, the failure interface is mimicked by a "soap-bubble" film that encloses the free surface at a specified volume V. Because the film is uniquely defined, it becomes important to assess the fitness, i.e., the closeness between the film and the real failure interface and to establish the error bounds of the method.
Under the current approach, the area, volume and scarp boundary are kept identical for each landslide, and hence, the mean depth of the predicted landslide mass will be equal to the actual one, (V/ p ), in principle. Nevertheless, because of the nature of discrete numerical data and computations, the discretization incurs digital errors, e.g., the DEM is presented at a resolution of grid size and in the minimal surface calculation, forward and backward interpolation is performed between the DEM and mesh system. It is expected that this type of error is bounded by the precision of the mesh, such that we list the normalized grid size¯ as an informative indicator, which is normalized with respect to the equivalent diameter 2R.
To compare the predicted and actual failure interfaces, we start by calculating and finding the crucial quantities in each of them: the mean and maximum scarp depths (h mean and h max , respectively) and the positions of the gravity center r GC and The parameters shaded in light colors indicate normalized non-dimensional assessment indices. EX1 and EX2 are ideal benchmark cases described in Appendix 1. maximum depth point r MD . One can define the measures of discrepancy as the differences of these quantities between the predicted and actual landslide failure interfaces. To eliminate the landslide scale effect, the depthwise quantities are normalized with respect to the actual mean depth; i.e.,Ē mean = (h pred mean − h act mean )/h act mean ,Ē max = (h pred max − h act max )/h act mean , and the spanwise offsets of the gravity centers d GC and maximum depths d MD are then normalized with respect to the equivalent diameter; i.e., The superscripts pred and act obviously represent the model predictions and the actual measurements of these quantities.Ē mean contains only the discretization error as discussed previously.Ē GC andĒ MD somewhat represent the bias of the depth distributions between the predicted and actual failure interfaces.
The ratioĒ max is one of the assessment indices, which merely measures the difference between the maximum depths. The other assessment indices are the statistical properties associated with the prediction discrepancy. Let the prediction discrepancy be the difference between the actual and predicted scarp depths, where (x, y) ∈ p . The statistical properties that can be computed include the distribution function of the prediction discrepancy δ, its mode(δ) (the most frequent discrepancy), standard deviation std(δ), and their non-dimensionalized counterparts,Ē mode andĒ std (normalized with respect to h act mean ). Note that std(δ) is virtually the root mean square error (RMSE) between the predicted and actual depths according to its computational principle . For completeness, the linear regression model, h act =βh pred + α, is also calculated for each landslide site. To facilitate comparisons among landslides of different size scales, the intercept α of the regression model is again nondimensionalized with respect to h act mean ; i.e.,ᾱ = α/h act mean . Together with the slopeβ and the coefficient of determination R 2 , these regression parameters are tabulated in Table 3.
The aforementioned statistical quantities are commonly applied in landslide studies. The limitations of these indices are that they rely on the fitness of a single point (e.g.,Ē max ,Ē MD ), on the averaged properties of the scarp (Ē mean ,Ē GC ), or on the discrepancy distribution of the scattered grid data (Ē mode ,Ē std ,ᾱ, β, R 2 , etc.). To the authors' point of view, one important factor has not yet been properly addressed by these indices, and this factor is the likelihood of the patterns between the predicted and actual scarps. The pattern refers to the landslide depth distribution in the proximity of any given point in the scarp area, and this pattern is highly dependent on the neighboring area. The patterns are omitted in the former statistical indices because  the elevation of each mesh grid is treated as an independent random variable.
To include in the spatial topographic patterns into the current assessment, we adopt the so-called structural similarity index (SSIM), which was developed for and is commonly used in image studies for assessing qualities among various image-processing schemes (Wang et al., 2004). In the computation of this index, the value of each pixel of the image, or here, those of the grids in the scarp area, is replaced by the SSIM calculated in the vicinity window (grid) surrounding the pixel. For each window, the SSIM is defined as (8) where the functions on the right hand side are defined as  Wang et al. (2004). The three expressions in (9) are for comparing the luminance, contrast and structural pattern in image studies but can be analogous to the comparison of the landslide scarp depth, the change in the depth and the pattern of the depth in the present study. Therefore, the SSIM should perform similarly as an assessment index to evaluate the fitness of the prediction. In the implementation, the SSIMs in the vicinity area of each grid are calculated with the help of the standard Gaussian window functions of signal processing and statistics. The window size is 10 DEM grids. The resultant SSIM for each landslide site is then defined as the average SSIM over the scarp area and takes a value between 0 and 1, where 0 indicates complete dissimilarity and 1 indicates perfect identicalness. Tests showed that the resultant SSIM does not sensitively depend on the window size and function. For further statistical theories and derivation details, interested readers are directed to the referenced literature.

APPLICATION EXAMPLES: FID5 AND FID18
In this section, application examples of the minimal surface method for two landslides are presented. There are also two additional benchmark cases for validating the minimal surface optimizations, but because of their simplicity, these benchmark cases are relegated to Appendix 1. The two real cases, one with a good prediction and the other with a fair prediction, are chosen and they are purposely selected to illustrate the comparison of the opposite predictions. The calculation of the statistical assessment indices are also explained in detail, particularly those operational procedures that are not thoroughly described in section 3 and section 4.
We start with the FID5 landslide. The FID5 site is on the slope of the west bank of the Cishan River and has a volume V of 1,167,000 m 3 and an average slope of 35.2 • , inclined toward the E-NE direction. The DEM domain is 600×340 m, in the x and y directions, respectively, with a resolution of 20 m (normalized to¯ = 7.36%). The 3D view of the site has been shown as a conceptual sketch in Figure 5. The post-landslide topography as well as the scarp boundary is plotted in Figure 8A. The projection area, p , of the scarp is then equal to 58, 000 m 2 and the equivalent radius R = 136.0 m. These geometric dimensions lead to a roundness ratio r r of 0.442 and a spherical ratio r s of 0.362.
As the slope and scarp boundary are now defined, we proceed with the preparation of the input data set of the minimal surface optimization scheme. Because FreeFEM++ only provides unstructured triangular meshing for general shapes of computational domains, numerical interpolation algorithms are employed to convey the data between the structured DEM mesh and the unstructured FreeFEM++ computational grids. We do not have any compulsory reasons to adopt higher order accuracy interpolations for the present purpose; thus only first order interpolation algorithms are used.
When the DEM and scarp outline are input, parsed and interpolated, an initial unstructured mesh for the scarp domain p is constructed 4 . The P2 finite elements are used for data arrangement and manipulation. Then, the optimization routine for the minimal surface is executed. To extend the scarp prediction to future mechanical slope stability analyses, which may require a higher accuracy within the FreeFEM++ framework, we perform a second-pass mesh adaptation based on the scarp depth of the initial calculation with a mesh refinement precision error factor of 0.05 4 . The refined mesh for the smooth minimal surface calculation is superposed in Figure 6A. A high mesh density is found in the areas where the topography exhibits large variations, in this case, around the vicinity of the ridges and edges of the slope surface. After the smooth minimal surface is found, it is interpolated back to the structured DEM grids.
The actual and predicted landslide depths are shown in These values lead to dimensionless discrepanciesĒ max = 32.40% andĒ mean = 0.40%, respectively. The prediction discrepancy distribution δ(x, y) is plotted in Figure 6D. The prediction overestimates the scarp depth in the part of the slope with higher elevation, such that δ(x, y) has two major oppositely signed zones aligned adjacently along the downslope direction. This distribution of δ(x, y) leads to the positions of both the maximum depth and gravity center of the landslide mass residing on the upper-slope side compared to the actual landslide mass. The positional offsets of the maximum depth and gravity center are d MD = 102.0 m (Ē MD = 37.53%), and d GC = 11.3 m (Ē GC = 4.17%).
The histogram of the prediction discrepancy δ(x, y) is plotted in Figure 7A. It depicts the frequency distribution of δ(x, y), and the distribution has an approximate symmetric triangular shape except for a minor peak at δ ≈ −9 m. Its mode, mode(δ), is approximately −0.7 m and the standard deviation, std(δ) is 7.3 m, leading to normalized discrepancy ratios ofĒ mode = 3.51% andĒ std = 34.51%. In Figure 7B, the scattered data of h pred and h act are plotted in the regression analysis and the resultant linear regression line is h act = 0.916h pred + 1.84, which provides the normalized intercept parameterᾱ = 0.087. The coefficient of determination, R 2 , is 0.618. The SSIM is 0.864, which is one of the cases with a high score, cf. section 6. In the SSIM computation, the default Gaussian window, of size 11×11 with a standard deviation of 1.5 cell sizes, is used. To compare with other cases, these statistical quantities are tabulated in both Tables 2, 3. FID18 is a landslide site near the Zion village. The landslide has a volume of 567,000 m 3 , an area of 46,800 m 2 , a maximum depth of 31.3 m, and a slope of 29.4 • inclined toward the south. The roundness and sphericity are 0.339 and 0.306, respectively. The DEM has a spatial resolution of 20 m (δ = 8.19%). After performing the smooth minimal surface calculation, the computational mesh, the actual, predicted depths and the prediction discrepancy are plotted in Figure 8. The panels are arranged the same way as in the previous FID5 case. Interestingly, the prediction also overestimated the landslide depth in the upper part of the slope and underestimated the lower part of the slope. Therefore, the prediction discrepancy also shows similar negative-positive-lobed prediction discrepancy zones along the downslope direction as in the previous case.
The accompanied discrepancy histogram and regression analysis are shown in Figure 9. The histogram indicates that the discrepancy has a flatter distribution with a few more irregular minor peaks than that of FID5. This case is identified as a poorer fit because the prediction also has relatively poor linear regression parameters compared to the actual data: the interceptᾱ (0.72, normalized with respect to h act mean ) is large, β (0.28) deviates greatly from 1, and R 2 (0.063) is low. We found consistent indications fromĒ mode (52.47%),Ē std (71.61%) andĒ GC (44.44%), representing larger prediction discrepancies, wider deviations and larger offsets of the gravity center. The SSIM (0.677) also has one of the lowest scores. Nevertheless, despite these comparatively subnormal indices, the prediction is only 2.31% forĒ mean and a satisfactory 20.61% forĒ max , as both the scarp area and volume are constrained in the present failure surface prediction method such that the mean depth has interpolation errors only and the maximum depth is generally proportional to the mean depth (see section 6).

DEPLOYMENT TO THE FULL SET OF LANDSLIDES AND ASSESSMENTS
The method is applied to the 24 landslide sites defined in Figures 2, 3 and Table 1. The resolution of the DEMs varies from 5 to 20 m due to various preparation conditions, such as public release policies and the size of landslide sites. Though with these different resolution settings, the normalized grid size¯ remains within 10% ( Table 2). In the same table, the maximum and mean depths as well as the positional offsets of the maximum depths and gravity centers between the actual and predicted scarps are listed. To make visual comparisons, the depth-related quantities are plotted in Figure 10. The predicted and actual depths excellently match the diagonal lines in Figure 10A. Reorganizing the data, we find that there is also a linear relationship between the maximum and mean depths ( Figure 10B). From the regression model, we have a gross guideline that the maximum scarp depth can be obtained by multiplying the mean scarp depth by a factor of 2.3.
The remaining assessment indices (i.e., the mode, standard deviation, regression parameters of the prediction discrepancy and the SSIM) are tabulated in Table 3. To comprehend at a glance the interrelationship among the assessment indices, we draw the box plots of the normalized indices in Figure 11. These indices are sorted in ascending order by their median values and are divided into two groups, of which the first group contains the slopewise quantities, normalized by 2R and the second involves the depth-related quantities, normalized by h act mean . The discretization is determined by the DEM resolution,¯ , and its values in Figure 11A, are within 10% for the present data-sets. The values of the normalized mean prediction discrepancy,Ē mean in Figure 11B, are small and are all within 10%, as expected. As argued in section 4, this finding is due toĒ mean being constrained by the specified area and volume inputs to the method and the small values arise from the DEM discretization approach and interpolation scheme.
There are two assessment indices associated with the maximum scarp depth: the positional offsetĒ MD and discrepancȳ E max . The two indices both exhibit much larger spreads of their data values than the other indices. The positional offset of the gravity centers of the predicted and actual scarps is interestingly small approximately 6% (median value), dimensionally of 6% · 2R ≈ 0.12R. An important implication of this fact for future incorporation of the mechanical slope stability analysis is that the force balance condition of the predicted scarp mass (cf. free body diagram) may not significantly differ from the actual mass, and consequently, the landslide motion dynamics may bear a close similarity. Research into this proposition is beyond the scope of the present paper and will be reported in follow-up studies.
The indexĒ mode depicts the most frequent prediction discrepancy, and from this definition, it somewhat indicates the skewness of the frequency distribution of δ(x, y), cf. Figures 7, 9. The moderate value range ofĒ mode ( Figure 11B) suggests that the frequency distributions of the δ(x, y) histograms remain reasonably symmetric. As mentioned in section 4, the standard deviation indexĒ std equivalently describes the RMSE between the predicted and actual scarps. The box plot shows that the deviation is approximately 45 ± 10% of the mean scarp depth, or approximately 19 ± 4% of the maximum scarp depth, based on the regression model depicted in Figure 10B. These margins ofĒ std indicate that the present smooth minimal surface is satisfactory for predicting the landslide failure interfaces. Figure 11C presents the value range of the resultant regression parameters and the SSIM.
Finally, we comment on whether there is a relationship between the good predictions from the present method and the landslide scarps. The importance of the answer to this question is that it enables us to estimate the goodness of the prediction without needing to know the actual underground failure interface when the method is applied to hazard mitigation plans for potential landslide sites. For this purpose, the roundness is chosen as the parameter to describe the landslide scarp    (Cleveland et al., 1992). The gray shaded areas indicate 95% confidence intervals. FID12 is an outlier and is excluded from the regression analysis.
FIGURE 13 | Flowchart of the present model for both application scenarios: verification, blocks shaded in light yellows, and hazard mitigation planning, blocks in light greens. The current failure surface prediction model is in light red colors. For failure surface prediction verification, i.e., the main theme of the paper, the computational flow is led by the thick solid brown lines but, for hazard mitigation planning, it is led by the thin dashed green lines.
because this parameter is usually the first obtainable information from topographic surveys. After reviewing the collection of the assessment indices, we selectĒ std , R 2 and the SSIM to measure the fitness of the predictions. These three indices vs. the roundness are plotted in Figure 12. The locally estimated scatterplot smoothing (LOESS) technique (Cleveland et al., 1992), is then applied to determine their relationships with certain degrees of statistical confidence. One exception is made for the noted outlier FID12, cf. section 2, which is excluded from the LOESS analysis but is included in the figure for reference. Comparing the three LOESS results, we can draw a consistent conclusion that the smooth minimal surface performs relatively poorer for a roundness of approximately 0.3. At this value, the standard deviationĒ std (RMSE) has the highest value whereas the coefficient of determination R 2 and structural similarity index SSIM have the lowest values. For values away from a roundness of 0.3, the landslide scarp becomes either slenderer or more circular. Intuitively, these two regimes have lower geometric complexity because they both degrade to simpler 2D-like scarps.
The prediction fitness may thus be associated with the reduction in the geometric complexity. In fact, the same conclusion can be reached with the undiscussedᾱ andβ parameters.
A note to keep in mind is that the conclusion drawn from the above assessments is based on the current simple smooth minimal surface approach for the failure surfaces. The geological settings and hydrogrological conditions are completely omitted. This arrangement is deliberate because it is essential at present to construct a baseline dataset with the simplest setting of the method. The dataset will be the foundation for future comparative studies if alternative failure surface prediction strategies are proposed. For example, the dataset will be used for improvement assessments when we design mathematical transformations to manipulate the failure surface predictions to accommodate site specific geological conditions. On the other hand, whether the current approach fits better for landslides in slopes with uniform rocks also remains to be investigated and the hypothesis can be examined by comparing to the baseline dataset. These related studies will be reported in upcoming papers.

CONCLUSION
Even with protection targets, a large number of potential deepseated landslides in the mountain ranges were identified in a series of renewed hazard mitigation initiatives in Taiwan. When implementing successful, detailed hazard mitigation strategies, tasks such as the determination of the landslide influencing areas and the installation of a monitoring system, need the estimation of the landslide volumes, failure interfaces and other related information. Having observed that these deep-seated landslide sites are represented by polygons in the GIS and that there are regression relations between the landslide scarp areas and volumes, we propose to use smooth minimal surfaces to approximate the landslide failure interfaces. The smooth minimal surface is constructed by minimizing the surface area while keeping the enclosed volume fixed at the value obtained by the scarp area-volume relation. This type of surfaces is chosen because it is closely related to the commonly used spherical-shaped surfaces in slope stability analyses. Consequently, it is expected that the prediction may suit better for slopes with homogeneous and isotropic rocks.
The method, though still in its primitive form, has the potential to be applied to hazard mitigation plans in which higher tolerance in prediction errors is allowed. Therefore, one of the main themes of the paper is to establish the knowledge about the fitness and margins of errors for this method setting. The fitness of the prediction results was verified with 24 landslides that were triggered by excessive rainfall during the 2009 Morakot typhoon. A collection of assessment indices was reviewed and among these indices, the standard deviationĒ std (equivalently, the RMSE), the regression parameters and the SSIM are shown to contain the information about the prediction discrepancy over the entire scarp domain for each landslide site. Using the present landslide dataset, the value range of the indexĒ std was found to be approximately 45 ± 10% of the mean scarp depth. The limited positional offset of the gravity center,Ē GC , indicates that the force balance condition may not significantly differ from the actual landslide mass. Finally, the relation between the prediction fitness and scarp geometry was determined: Better predictions are achieved for either slender or comparatively circular scarps. Overall, the indices reveal that the smooth minimal surface method is able to produce practical, acceptable predictions of deep-seated landslide failure interfaces, despite the omission of stratigraphy, geological structures and hydrogeological conditions.
The method is implemented with ease of use in mind such that operational costs are kept low and automation processes for a large number of landslide sites can be made possible. An additional benefit is that the method can simultaneously generate 3D computational meshes for the landslide scarp mass. These meshes can be applied to many subjects related to hazard mitigation plans. For example, for slope stability management, it is currently ongoing to develop 3D extended schemes for landslide volume and safety factor relations. Integrating with surface deformation data and material constitutive laws, slip displacements on the failure surface can be inversed for deep-seated landslides, in a somewhat similar way to dislocation model in tectonics. The volumetric mesh, or depth distribution, can also be used as the initial mass to perform landslide motion simulation to assess the influenced area. Summarizing the abovementioned potential applications and verification details presented in the present paper, we draw a combined flowchart in Figure 13 for the two application scenarios of the method. The inputs and products of the method are explicitly indicated for each application scenario.

DATA AVAILABILITY STATEMENT
All datasets and computational codes generated for this study are included in the article/Supplementary Material.