Development of Potential Slip Surface Identification Model for Active Deep-seated Landslide Sites: A Case Study in Taiwan

1. Introduction

Deep-seated landslides often result in severe damage due to the greater depth to slip surface and volume of soil and debris mass movement [1]. The massive falling materials cause problems such as burying houses, blocking mountain roads, and obstructing downstream flood-prone areas by clogging river channels [2]. For example, one of Taiwan's most staggering deep-seated landslide disasters in the past decades occurred in August 2009 [3]. A rainfall event with a cumulative rainfall exceeding 2500 mm triggered a large-scale slope collapse disaster, leading to the destruction of the village of Xiaolin and the deaths of 474 people. This sudden and massive slope disaster had a slip surface depth of approximately 80 meters. As a result of such a significant rainfall event, the study of deep-seated landslides has become prevalent in Taiwan.

The failure area and depth of a landslide determine the magnitude of the disaster. When the depth of the slip surface is greater, the failure surface area will also expand with depth. Therefore, the confirmation of the depth of the slip surface plays a crucial role in predicting the scale of the disaster. In terms of mitigating such landslide disasters, having early knowledge of the depth of sliding and the affected surface area allows for the timely implementation of appropriate engineering remediation and disaster preparedness efforts. In order to obtain information concerning the locations of the slip surface, in engineering practice, the most common method for monitoring potential locations of the slip surface is to drill a borehole with the installation of an inclinometer [5]. The inclinometer is an in-situ instrument for displacement monitoring with time that can indicate when and where cumulative slope deformations could happen. Although the inclinometer has high measurement accuracy, this technique requires labor-intensive manual observation and does not have continuous monitoring with time. To overcome this drawback, techniques with real-time monitoring have been developed, including in-placed inclinometer (IPI) [5], ShapeAccelArray (SAA) [6], Time Domain Resistivity (TDR) [7], and Fiber Bragg Grating (FBG)[8]. Compared to the traditional sliding inclinometer, the above four types of instruments and techniques have been successfully introduced to the real-time and distributed monitoring of landslide displacement. However, there may be differences in measurement accuracy, resolution for larger shear deformations, and limitations in measurable deformations among these four methods [9].

Although the aforementioned on-site investigation methods can directly measure the depth of the sliding surface, information about the sliding surface depth is only obtained when the slope undergoes actual changes. Before the deformation trigger of the collapse mass, the scale of the landslide disaster cannot be determined. This situation leads to the inability to formulate prevention and preparation plans for disasters in advance to cope with potential disaster impacts and losses in the future. Therefore, some studies focused on the advance prediction of the disaster scale caused by landslides. The simplest method is the core drilling and visualizing method, which determines the occurrence of the sliding surface based on the lithology of the strata. According to the comparison of past landslide cases with the actual occurrence positions and lithology, it is usually found that the sliding surface is most likely to occur at the soil-rock interface or the location of sheared mud. This perspective is based on the significant differences in the permeability characteristics of the upper and lower formations of the sliding interface. When rainfall infiltrates the slope, groundwater accumulates above low-permeability formations, leading to increased soil weight and reduced internal friction. When the sliding force exceeds the shear strength of the soil, it triggers the sliding of the formation. Therefore, determining the sliding surface location based on the concept that the permeability of rock formations plays a crucial role in groundwater seepage can allow for an early understanding of the potential depth of the sliding surface. However, the drawback of this method is that it is prone to errors due to subjective experiential judgment, and visually determining the permeability characteristics of lithology may also lead to misjudgments.

Additionally, numerical analysis methods are common for predicting potential sliding surface locations in advance. Numerical methods have evolved over time and mainly include the limit equilibrium method (LEM) and numerical modeling methods (finite element, finite difference, and boundary element methods) [10,11,12,13]. However, these methods all have some basic assumptions and require the preparation of hydrogeological models for a collapse site, which consumes more time and human resources. Furthermore, the accuracy of predicting sliding surface locations highly depends on the accuracy of the hydrogeological models built from field investigation data. Therefore, some studies [14,15,16] proposed new methods for predicting sliding surfaces, attempting to overcome the limitations of numerical methods.

The basic concept of the newly proposed method was to integrate core sample record data and core testing data and to compare the data with the actual sliding surface positions, thereby confirming effective reference signal data for determining sliding surface positions. Furuya and Jiang [17] integrated various data, including RQD (rock quality designation), the Equotip hardness values, unit weight, the magnetic susceptibilities, orientation of geologic discontinuities, and digital imaging of core samples, from two deep-seated landslide sites in Japan for determining potential slip surfaces. The investigation results indicated that most data types, except for digital imaging data and magnetic susceptibility, appeared to have a distinct change near the observed slip surfaces for both sites. The Geological Survey and Mining Management Agency attempted to use instruments such as magnetic susceptibility meters and natural gamma ray meters to quantify the physical properties of core samples, aiming to study and evaluate their feasibility and applicability in identifying potential sliding surfaces [18]. Investigation results from a landslide site in southern Taiwan indicated that, regarding natural gamma ray measurements, there were slight variations in values depending on the composition of the rock mass, but overall, there was no significant difference. Additional data were needed to assess potential sliding surface locations if these measurements were to be used. In terms of magnetic susceptibility measurements, it was noted that variations in core magnetic susceptibility were observed due to differences in material composition; however, significant changes were not observed near the sliding depth in this particular case. Therefore, when the lithology of landslide sites was relatively homogeneous, the sensitivity of both instruments mentioned above was limited, posing challenges in evaluating potential sliding surfaces.

In addition to the emerging methods mentioned above for estimating the location of sliding surfaces, some scholars have observed a correlation between pore water pressure and the sliding surface of the stratum. For instance, Vassallo et al. [19] established a three-dimensional numerical model based on long-term rainfall data, pore water pressure, and permeability investigation results, finding that the permeability and pore water pressure of the sliding block are both higher than those of the stable layer beneath the sliding block. Park and Song [20] observed changes in rainfall and pore water pressure in a homogeneous sandy soil slope model, indicating that excessively high pore water pressure and hydrostatic pressure were more likely to cause slope failure. C. Di Maio et al. [21] studied the variation of pore water pressure and permeability in a mudstone landslide site using constant head tests and three-dimensional simulation methods (MODFLOW). The results showed a significant increase in pore water pressure in the sliding zone during rainfall, while the variation of pore water pressure below the sliding surface was stable. On the other hand, there was a significant difference in permeability between the sliding zone and the rock formation below the sliding surface. This variation in permeability can also be found in other landslide-sliding surfaces [22,23,24]. Therefore, these studies imply that if the vertical distribution of permeability characteristics in the strata can be obtained from locations with variations in the coefficient of permeability, this concept appears to have the potential to determine the occurrence location of the sliding surface.

If the use of borehole permeability profile data to predict potential sliding surface locations is a viable approach, the primary task is to address the method of obtaining profile data. Typically, the most direct and accurate method for obtaining such data is through dense-sampling double-packer hydraulic tests. However, this technique is associated with high costs, lengthy testing times, and substantial workforce requirements. To obtain permeability data continuously with depth at a lower cost and more efficiently, some scholars have proposed influential factors related to permeability and established relationships between these factors and permeability coefficients, thereby developing empirical formulas for quickly obtaining permeability profile data. A detailed review of various empirical formulas can be found in the study of Hsu [25]. Most empirical formulas are only applicable to single lithology sites and adopt fewer factors related to permeability characteristics, rendering them unsuitable for complex lithology sites. However, in recent years, Hsu [25] developed a rock mass permeability estimation model (HCB2 model) using composite geological indices applicable to various geological environments (sedimentary, igneous, and metamorphic rocks). Although this model can be applied to different lithology sites, the analysis data used to establish the model were not obtained from deep-seated landslide potential sites. The lithology of landslide sites mainly consists of fractured rock layers and colluvium, and the hydrogeological characteristics of these sites may differ from non-landslide sites, warranting further investigation and comparison.

Therefore, to obtain detailed and continuous hydrogeological information from deep-seated landslide potential sites and meet the needs of disaster prevention management, this study used borehole cores, downhole imaging, and double-packer hydraulic test data from 24 landslide-prone sites in Taiwan to construct a new rock mass permeability estimation model only applicable to deep-seated landslide potential sites. Prediction results from the new model were compared with the HCB2 model. Finally, based on the above results and with displacement data monitored by inclinometer meters, a new model for confirming sliding depth was developed to predict the scale of landslide disasters rapidly.

2. Description of Study Area and Data Preparation

Based on the needs of this study in establishing models for estimating the permeability of landslide rock masses and identifying the location of sliding surfaces, it is necessary to collect relevant data from investigation projects of potential deep-seated landslide areas in the mountainous regions of Taiwan to assist in the construction of the two models. The required data include (1) hydrogeological survey data of the strata where landslides occurred: drilling core samples, borehole imaging data, groundwater level data, and double-packer hydraulic test data, which can serve as the primary data for constructing models to assess the permeability characteristics of landslide rock masses; (2) underground displacement data: time-series monitoring data of inclinometers in landslide areas, serving as investigating the relationship between the simulated distribution of hydraulic conductivity and the actual occurrence locations of sliding surfaces.

Reviewing the current landslide investigation projects in Taiwan, the most comprehensive project is under the jurisdiction of the Geological Survey and Mining Management Agency (GSMMA). From 2007 to 2013 [26,27,28,29,30,31,32], the GSMMA implemented the project entitled "Investigation and Evaluation of Hydrogeological Impacts on Slope Stability in Watersheds across Taiwan." This project mainly focused on investigating and assessing hydrogeological impacts on slope stability in major watersheds across Taiwan. It aimed to strengthen field hydrogeological investigations to understand the extent of hydrogeological impacts on landslide occurrence. Furthermore, it conducted research on the correlation between landslides and related factors. The project targeted deep-seated landslide areas within major northern, central, southern, and eastern Taiwan watersheds. The 24 landslide sites were selected for investigation, experimentation, monitoring, and analysis to understand how infiltration during rainfall alters hydrogeological conditions, leading to a potential trend of recurring rockslides. At each landslide site, three geological boreholes were drilled with depths ranging from 50 to 80 meters, and geological descriptions were prepared based on core samples. Additionally, based on the results of borehole imaging surveys, representative sections of fractured rock masses were selected for double-packer hydraulic tests, totaling 169 test samples with a test section length of 1.5 meters each. Finally, an inclinometer was installed within the deepest borehole at each landslide site to monitor underground displacement, while other boreholes were set up as groundwater observation wells for monitoring groundwater levels. However, some inclinometer boreholes were constructed to monitor groundwater levels simultaneously.

Furthermore, the samples selected for constructing the specific rock mass permeability estimation model did not come from data of a single rock type but encompassed a variety of lithologies. The primary lithologies include sedimentary and metamorphic rocks, while sub-lithologies cover sandstone, shale, alternations of sandstone and shale, sandy shale, siltstone, muddy siltstone, mudstone, argillite, slate, and phyllite. Due to the wide range of lithologies covered by the data samples, the new rock mass permeability estimation model developed in this study possesses the characteristic of adapting to different geological environments.

Finally, the primary data collected for the 24 landslide sites, including inclinometer borehole number, landslide site name, primary lithology, regolith thickness (depth to rock mass), borehole depth, and normal groundwater level, are organized as shown in Table 1. Among them, normal groundwater levels are primarily collected from monitoring wells after the completion of inclinometer borehole installation and a few days of settling without any occurrence of heavy rain. If the inclinometer boreholes do not include groundwater monitoring capabilities, groundwater levels are estimated from data collected at nearby groundwater monitoring stations. Previous studies have suggested a correlation between groundwater levels and slope stability [33,34]. Therefore, this study also used the collected groundwater level data as one of the key factors for developing landslide slip surface identification models.

Figure 1. Locations of the 24 landslide sites.

Figure 1. Locations of the 24 landslide sites.

3. Methodology

This study targeted mountainous areas of Taiwan with deep-seated landslide potential. Utilizing past hydrogeological survey data, a model for estimating rock mass permeability was constructed. The aim was to reveal the distribution of permeability characteristics in the vertical direction of rock formations and attempt to explain the mechanism of landslide occurrence from the perspective of permeability distribution patterns. The research findings were expected to provide rapid feedback on disaster scale information for landslide-prone areas with deep-seated landslide potential, contributing to disaster prevention efforts. The associated approaches and procedures, including (1) selecting influential factors for characterizing the hydraulic properties of fractured rocks in active landslide sites; (2) multicollinearity analysis for the selected influential factors; (3) principal component analysis (PCA) for determining the weight of each factor; (4) model establishment and evaluations; (5) developing a rule for identifying potential slip surface, are provided as follows.

3.1. Selection of Influential Factors

Prior to developing a rock mass permeability estimation model for landslide sites, influential factors that may affect the permeability of disturbed rock mass must be considered. Based on the material composition and structural characteristics of fractured rock masses in landslide areas, this study referred to previous research [35,36,37,38,39,40,41,42,43,44] to select potential factors related to rock mass permeability characteristics. The most commonly considered factors include depth index (DI), rock quality designation (RQD), fracture aperture (FA), and fracture density (FD). Earlier models were developed based solely on individual influencing factors. However, subsequent research has considered multiple factors in model development [39,42,45], and the research results have demonstrated that predictive models using multiple factors outperform those using single factors. However, the models established above were often limited to application in single lithology sites. In recent years, Hsu [25] attempted to develop predictive models applicable to various geological environments, considering up to eight influencing factors. Eventually, six different combinations of factors were used to build predictive models. Although Hsu [25] has improved the accuracy and limitations of predictive models, the models were built using data collected from non-collapsible strata in the mountainous areas of Taiwan. However, hydrogeological characteristics are significantly more complex for sites with high collapsibility potential compared to non-disturbed strata conditions. Currently, there is no systematic construction of similar predictive models. Therefore, this study, based on survey data collected from sites with high collapsibility potential in Taiwan and referring to Hsu's [25] six-factor model, selected six influencing factors, including RQD, DI, FA, FD, gouge content designation (GCD), and lithology permeability index (LPI), to construct a predictive model for the permeability of disturbed strata.

In brief, RQD measures the quality of the rock mass, with lower values indicating more fractured rock, which correlates with higher permeability. DI represents the influence of in-situ stress on rock mass permeability, theoretically leading to decreased permeability with depth. GCD calculates the content of shear mud in a specific section of the rock mass, aiming to correct the drawback of RQD in neglecting shear mud, thereby identifying a critical factor in determining whether fractures have groundwater-conducting properties. LPI represents the overall permeability of the intact rock mass and is related to the rock's pore characteristics, excluding the influence of fracture characteristics on permeability. Finally, FA and FD describe the structural features of the rock mass and are crucial indicators of rock mass permeability, indirectly influencing the groundwater flow. Detailed scoring definitions and criteria in response to permeability for each influencing factor can be found in Hsu [25] and will not be repeated in this paper.

However, regarding the scoring system for the LPI factor across various lithologies, Hsu [25] initially designed it based on core samples obtained from relatively undisturbed rock mass sites. However, the sites of this study mainly consist of disturbed geological strata. According to the core investigation data from the collapse sites, lithologies such as weathered slate, weathered shale, weathered argillite, and argillite frequently appear, which were not accounted for in the original LPI scoring system. Therefore, this study made partial modifications to the LPI scoring system. In Table 2, in addition to the LPI score table developed by Hsu [25], additional rock mass scoring values specific to collapse sites were added (as shown in the last four rows of Table 2). The scoring ranges from 0 to 1, where a higher score indicates a higher permeability. Lastly, due to the complexity of lithologies in collapse areas, the calculation of LPI values in double-packer test sections may involve different lithological proportions, and therefore, the calculation was performed by combining various lithological ratios.

3.2. Multicollinearity Analysis for Selected Influencing Factors

Since this study employs multiple-factor statistical analysis, it is necessary to conduct multicollinearity tests on the selected influencing factors to examine whether there are multicollinearity issues among them. Multicollinearity can be assessed using one commonly used indicator-tolerance and variance inflation factor (VIF). Tolerance refers to the percentage of variance in a predictor variable that is not explained by other predictor variables, with values ranging from 0 to 1. A tolerance value close to 1 indicates the absence of multicollinearity, while values closer to 0 suggest multicollinearity issues. However, there is no formal threshold to determine the presence of multicollinearity among factors [46]. Although no universally accepted standards exist, some studies have explored the relationship between tolerance values and multicollinearity. For instance, Myers [47] suggested that tolerance values below 0.1 indicate severe multicollinearity issues, while Menard [48] proposed that values below 0.2 indicate potential multicollinearity problems.

Additionally, the reciprocal of tolerance, known as the variance inflation factor (VIF), is used to indicate the extent to which the variance of coefficient estimates is influenced by multicollinearity. Like tolerance, there is no formal cutoff value for VIF. However, VIF values exceeding ten indicate multicollinearity issues, although this criterion may require smaller values for weaker regression models. For example, Midi et al. [46] recommended that VIF values exceeding 2.5 in logistic regression analyses warrant attention to multicollinearity issues. Finally, this study conducted multicollinearity tests on the six influencing factors mentioned in Section 3-1 using the aforementioned diagnostic metrics to guide model construction.

3.3. Weights of Influential Factors through Principal Component Analysis

Theoretically, the permeability characteristics of rock masses are controlled by various influencing factors, and these factors may have different impacts on the permeability characteristics of fractured rock masses. By assigning appropriate weights to the influential factors and then combining the factors, better correlations with permeability characteristics may be obtained. To extract the weight values of each influential factor, this study adopted the Principal Component Analysis (PCA) method, a multivariate analysis technique.

PCA can be conducted using either an unstandardized covariance matrix or a standardized correlation matrix. When the units of the analyzed variables differ, using the correlation matrix is generally preferable. The eigenvalues and eigenvectors are then calculated from this matrix. The variance of each principal component is equal to its corresponding eigenvalue, with the first principal component having the most considerable variance, followed by the second principal component, and so on. The eigenvalues can also be used to compute the contribution proportion of each principal component (the percentage of total variance accounted for by that principal component), which reflects the explanatory power of the entire set of variables and serves as a basis for selecting principal components. The eigenvectors are the coefficients in the linear combinations of the principal components, and the coefficients can be derived from the loadings. Loadings are the correlation coefficients between the variables and the principal components., calculated using the eigenvalues, eigenvectors, and standard deviations. The loading formula is shown in Equation (1).

$$r_{y_j{x}_i}={\displaystyle \frac{a_{ji}\sqrt{{\lambda }_j}}{s_i}}$$

(1)

where the loading $r_{y_j{x}_i}$ is the correlation coefficient between the j-th principal component and the i-th variable. A larger loading value indicates a greater influence of the principal components on the variables; ${\lambda }_j$ is the eigenvalue of the j-th principal component; $s_i$ is the standard deviation of the i-th variable; and $a_{ji}$ is the eigenvector coefficient.

Therefore, the eigenvector coefficient $a_{ji}$ can be derived from Equation (1) as follows.

$$a_{ji}={\displaystyle \frac{r_{y_j{x}_i}}{\sqrt{{\lambda }_j}}}{s}_i$$

(2)

After obtaining the eigenvector coefficients, the linear combination coefficients of each principal component can be derived. The next step is to combine the linear combination coefficients of each principal component with their variance contribution ratios. Using the variance contribution ratio of each principal component as weights, the eigenvectors and variance contribution ratios of the more explanatory principal components are integrated to calculate the new coefficients of each variable in the new linear combination model.

The selection rule for the more explanatory principal components refers to selecting the first n principal components with the most significant variance contribution ratios among all j principal components, with the cumulative variance contribution rate of these n principal components exceeding 80%. Assuming that the first n principal components represent the original variable indicators, the new coefficients of the variables can be considered as a weighted average of the old coefficients of the variables in the linear combination of these n principal components, using their variance contribution ratios as weights. Therefore, the coefficient matrix b of all variables in the composite score model is given as follows:

$$\mathrm{b}=\left[\begin{array}{c}{b}_1\\ b_2\\ \begin{array}{c}⋮\\ b_i\end{array}\end{array}\right]={\displaystyle \frac{a_1\times c_1+a_2\times c_2+\dots +a_n\times c_n}{c_1+c_2+\dots +c_n}}$$

(3)

where $a_n$ represents the eigenvector matrix of the n-th principal component; $c_n$ represents the contribution ratio of the n-th principal component.

Finally, since the sum of the coefficients (weights) for all variables i must equal one, all variable weights need to be normalized based on the variable coefficients in the composite score model. This results in the final weight score $w_i$ for each variable given as follows.

$$w_i={\displaystyle \frac{b_i}{b_1+b_2+\dots +b_i}}$$

(4)

By using the above analysis process with available hydrogeological investigation data, the number of principal components with significant contributions can be selected and combined with the contribution proportions to calculate the weight value for each influencing factor. The determined weights can be utilized to establish the rock mass permeability index and further construct a rock mass permeability estimation model.

3.4. Model Establishment and Hypothesis Testing

As discussed in Section 3.1, six geological factors may influence the hydraulic properties of rock mass. Each factor impacts the permeability values differently. For example, if the RQD value of the rock core shows poor rock mass quality, theoretically, the permeability should be higher. However, if there is gouge filling in the fractures, the permeability of the rock mass may decrease. Similarly, suppose two core sections have the same RQD value but are at different depths. In that case, the deeper section may have reduced permeability due to increased stress closing the joints or rock matrix porosity. Therefore, a comprehensive permeability index for fractured rock masses may need to integrate various geological indicators. The six geological indicators (RQD, DI, GCD, LPI, FA, and FD) were integrated into this study. Considering the principal component analysis in Section 3.3 to obtain the weight of each geological factor, the factors were linearly combined by multiplying each factor by its corresponding weight value of w_i. Equation (5) defines the hydraulic conductivity potential index (HCPI) for landslide-prone bedrock. The formula for HCPI is given as follows:

$$\mathrm{H}\mathrm{C}\mathrm{P}\mathrm{I}=\left[\right(1-{\displaystyle \frac{RQD}{100}})\times w_1]+(\mathrm{D}\mathrm{I}\times w_2)+\left[\right(1-\mathrm{G}\mathrm{C}\mathrm{D})\times w_3]+(\mathrm{L}\mathrm{P}\mathrm{I}\times w_4)+(\mathrm{F}\mathrm{D}\times w_5)+(\mathrm{F}\mathrm{A}\times w_6)$$

(5)

HCPI can quantify the permeability of the rock mass with the scores of the given six factors, except for FA, ranging between 0 and 1. Higher values indicate a greater contribution of the corresponding parameter to increasing the rock mass permeability. A higher HCPI value represents permeable rock masses in the test section. Subsequently, different mathematical models (such as simple linear, logarithmic, exponential, or power law) can be selected, and regression analysis techniques can be applied to establish an estimation model for the permeability of fractured rock masses in landslide formations. The goodness of fit for each model can be determined by the R² value, the correlation coefficient between permeability and HCPI.

Additionally, two statistical analyses were conducted to ensure a reliable regression result for the best model. The first statistical analysis inspected the regression model’s assumptions for residuals (independence using the Durbin-Watson statistic, constant variance using Spearman’s rank correlation, and normality using the Kolmogorov–Smirnov test). In contrast, the second statistical analysis investigated the statistical significance of the regression analysis model using the F test statistic and p values. Based on the regression analysis outcomes, the best model for predicting the permeability of fractured rock masses was decided. The regression and statistical analyses mentioned above were carried out using the commercial software SigmaPlot.

3.5. Developing a Rule for Identifying Potential Slip Surface

To summarize and clarify whether the distribution pattern of permeability in a landslide area is related to the occurring location of the sliding surface, this study used the landslide-specific rock mass permeability estimation model established in Section 3.4 to create continuous hydraulic conductivity profiles for the 24 landslide sites, showing the variation with borehole depth meter by meter. Next, the observed inclinometer data from the 24 landslide sites, which had clear records of sliding, were superimposed onto the hydraulic conductivity bar charts (as illustrated in Figure 2). This work helped observe whether there were significant changes or other trends in hydraulic conductivity values near the monitored sliding depths. This process explored the feasibility of identifying sliding surface locations using information on the hydraulic conductivity profile. Subsequently, a more comprehensive method for estimating the potential sliding depth was established by combining borehole core data, groundwater level data, and the relationship between these data and the sliding surface locations. This method was then validated by applying it to sites with known sliding depths to check if it could correctly identify the sliding surface locations. Finally, the identification rule developed in this study can be applied when inclinometers have not yet detected significant displacement. This application can assist disaster prevention decision-making agencies in estimating the potential scale of deep-seated landslides in advance.

4. Results and Discussion

4.1. Descriptive Statistics for Six Influential Factors

Establishing the permeability estimation model requires first integrating the hydraulic conductivity of existing double-packer test sections and calculating the scores of the influential factors for each test section for subsequent principal component analysis and regression modeling. This study compiled 169 double-packer test section data from 48 boreholes in the 24 landslide areas. It sequentially calculated the RQD, DI, GCD, LPI, FA, and FD values for each double-packer section according to the method described in Section 3.2. The descriptive statistical results of the calculated factor scores are shown in Table 3. The overall average RQD value for all lithologies is approximately 0.479, with sedimentary rocks averaging 0.656 and metamorphic rocks averaging 0.291. According to Bieniawski's [49] rock core quality determination method, the overall rock core quality is poor, with fragmented metamorphic rock formations. Although the rock core quality of sedimentary rock formations is better than that of metamorphic rocks, it only reaches a fair level according to Bieniawski's [49] criteria. These results indicate that the overall rock core quality in landslide areas is generally poor, as expected.

Further observation of the statistical analysis results regarding FD and FA show that the average values for both parameters in metamorphic rock formations are greater than those in sedimentary rock formations, indicating a more developed fracture network in metamorphic rocks, which is consistent with the RQD results. Moreover, the descriptive statistics concerning GCD indicate that the gouge content in sedimentary and metamorphic rocks is comparable. Since DI and LPI represent the relative depth and lithology of the sampling section, their mean values and standard deviations have less physical significance and are not discussed here.

Finally, examining the six geological characteristic indicators selected for this study reveals that higher RQD and GCD values indicate better core quality and higher gouge content in the rock core section. Both higher values lead to lower hydraulic conductivity K, meaning RQD and GCD are negatively correlated with K. Conversely, the remaining four indicators are positively correlated with K. To ensure that all indicator factors have a positive correlation with K, this study transformed the RQD and GCD values to (1-RQD) and (1-GCD). This transformation changes the correlation to positive without affecting the order and spacing of the values. After the transformation, all factors show a positive correlation with K, thereby enhancing the convenience of interpreting the results of subsequent multicollinearity analysis and principal component analysis.

4.2. Multicollinearity Analysis Results of the Selected Factors

When establishing regression models using possible geological factors, checking for multicollinearity among them is vital to prevent redundancy in the independent variables. This could artificially enhance the explanatory and predictive power of certain variables, leading to inaccurate regression model predictions. Therefore, this study used the Variance Inflation Factor (VIF) to assess multicollinearity. The criteria for VIF are as follows: when VIF is greater than ten, it indicates a strong multicollinearity problem. The results of the multicollinearity analysis are shown in Table 4. Upon examining the VIF column in Table 4, it is evident that the VIF values for all factors are less than ten, indicating no strong multicollinearity among the indicator factors. Therefore, the regression model can directly include the six geological factors—RQD, DI, GCD, LPI, FA, and FD.

4.3. Results of PCA

This section presents the analysis results of calculating the weight values for each factor using PCA. First, the varimax method was used for axis rotation, and more influential factors were extracted based on principal components with eigenvalues greater than 1, aided by a scree plot to determine the number of components to extract (as shown in Figure 3). The analysis was performed with a cumulative contribution ratio of over 80%.

Table 5 shows the calculation results of each component in the principal component analysis, displaying the eigenvalues, contribution ratios, and cumulative contribution ratios of the principal components. Larger eigenvalues indicate greater importance. As shown in Table 5, the eigenvalues of the first three principal components are greater than 1, with a cumulative contribution ratio of 0.8152 for these three components. The cumulative contribution ratio explains 81.52% of the data variance.

Thus, the first three principal components with their variance contribution ratios were utilized to determine the weights of all factors. First, the factor loadings were derived according to Equation (1), as shown in the upper part of Table 6. Then, the eigenvector coefficients for the three principal components were calculated using Equations (2). The resulting linear combination equations for the three principal components are as follows:

$$\mathrm{F}1=0.589\times \left(1-\mathrm{R}\mathrm{Q}\mathrm{D}\right)+0.058\times \mathrm{D}\mathrm{I}+0.116\times \left(1-\mathrm{G}\mathrm{C}\mathrm{D}\right)-0.068\times \mathrm{L}\mathrm{P}\mathrm{I}+0.505\times \mathrm{F}\mathrm{A}+0.582\times \mathrm{F}\mathrm{D}$$

(6)

$$\mathrm{F}2=-0.14\times \left(1-\mathrm{R}\mathrm{Q}\mathrm{D}\right)-0.474\times \mathrm{D}\mathrm{I}+0.602\times \left(1-\mathrm{G}\mathrm{C}\mathrm{D}\right)+0.622\mathrm{L}\mathrm{P}\mathrm{I}+0.251\times \mathrm{F}\mathrm{A}+0.015\times \mathrm{F}\mathrm{D}$$

(7)

$$\mathrm{F}3=-0.103\times \left(1-\mathrm{R}\mathrm{Q}\mathrm{D}\right)+0.685\times \mathrm{D}\mathrm{I}-0.223\times \left(1-\mathrm{G}\mathrm{C}\mathrm{D}\right)+0.623\times \mathrm{L}\mathrm{P}\mathrm{I}+0.285\times \mathrm{F}\mathrm{A}-0.047\times \mathrm{F}\mathrm{D}$$

(8)

Third, the weighted average of the linear combination coefficients of the three principal components was calculated to obtain the coefficient of the composite score model using Equation (3), as shown in Table 7. The composite model Equation is given as follows:

$$\mathrm{Y}=0.247\times \left(1-\mathrm{R}\mathrm{Q}\mathrm{D}\right)+0.049\times \mathrm{D}\mathrm{I}+0.174\times \left(1-\mathrm{G}\mathrm{C}\mathrm{D}\right)+0.264\times \mathrm{L}\mathrm{P}\mathrm{I}+0.39\times \mathrm{F}\mathrm{A}+0.296\times \mathrm{F}\mathrm{D}$$

(9)

Finally, based on the requirement that the weights of all factors need to be normalized in the composite score model, the weight values for each factor were obtained, as shown in Table 8.

The results of the weight analysis indicate that FA and FD have higher weight values, demonstrating that these two factors are the most significant in landslide areas, implying that the extent of the fracture network has a greater impact on assessing the hydraulic conductivity of fractured rock masses in landslides. Additionally, except for DI, which has a relatively low weight, the other factors are more balanced. The lower influence of DI in the principal component analysis suggests that, given the maximum borehole depth of 80 meters, the effect of tectonic stress is less significant within the limited depth range.

4.4. Establishment and Evaluation of the Predictive Model

Based on the results of the weight analysis on six geological factors through PCA in Section 4.3, each factor was multiplied by its respective weight and then linearly added to obtain the integrated permeability indicator HCPI. Next, by applying common mathematical models (including simple linear, logarithmic, exponential, or power law models) for regression analysis, the model was constructed using double-packer hydraulic test data and the HCPI data calculated from core samples. The coefficient of determination (R²) was used to identify the optimal regression mathematical model, thereby deriving the permeability estimation model for landslide-prone rock masses.

The regression analysis results, as shown in Figure 4, indicate that the power-law formula best fits the relationship between HCPI and K, with an R-squared value of 0.8953, demonstrating good predictive capability.

Furthermore, the optimal regression model underwent an F-test to examine the overall significance of the model and whether the error terms satisfy the assumptions of normality (Kolmogorov-Smirnov statistic), independence (Durbin-Watson statistic), and homoscedasticity (Spearman rank correlation coefficient).

The validation results of the best estimation model are shown in Table 9. The F-value and P-value from the regression analysis indicate that the independent variables of the hydraulic conductivity estimation model are significant predictors of the dependent variable, demonstrating predictive capability. Moreover, the model satisfies the three assumptions for the validation of error terms: the Durbin-Watson statistic of 1.5 (less than 2) indicates no autocorrelation in the residuals; the K-S normality test value is greater than 0.05, indicating the data follows a normal distribution; and the homoscedasticity value is 0.58, suggesting a good correlation between factors and permeability coefficient.

These five validation tests confirm the effectiveness of the permeability estimation model for landslide-prone rock masses. Therefore, in future applications, the developed landslide rock mass permeability estimation model can effectively and cost-efficiently obtain detailed and continuous hydraulic conductivity profiles. This work can be achieved at any active landslide site by obtaining six geological indicator parameters from core samples and double-packer hydraulic test data.

4.5. Comparisons of Previous and New Developed Models

In addition to developing a permeability estimation model for landslide-prone rock layers (named HCPI model), another research objective of this study was to test the differences in hydrogeological characteristics between disturbed and undisturbed fractured rock formations. The disturbed model is the HCPI model, while a brief introduction to the undisturbed model is provided here. Hsu [25] collected hydrogeological investigation data from boreholes in the mid and upper river basin of the Choshui river in Taiwan and developed an integrated permeability estimation model (named HCB2 model) using geological characteristic indicators and hydraulic conductivity data from double-packer test sections. This study successfully proposed a model applicable to the complex and fractured geological environments in the mountainous regions of Taiwan. While Hsu’s model [25] is applicable to sites of different lithologies with a wide application range, the data collected in their study did not specifically consider sites with the potential for stratum sliding. The hydrogeological characteristics of such potentially active landslide sites may differ from typical mountainous sites in Taiwan. Therefore, this study intented to compare hydrogeological properties between HCB2 and HCPI models. Before conducting this comparison, this study expanded the HCB2 model by incorporating more data points. In addition to the original 104 sample points, 352 data points from undisturbed strata in other major mountain basins of Taiwan were collected. The HCB2 model was then reconstructed to make the estimation of hydraulic conductivity more representative (named NHCB2).

Figure 5 shows the comparison of hydrogeological properties between NHCB2 and HCPI models. As shown in this figure, both models cover common metamorphic and sedimentary lithologies, indicating a wide range of lithologies. However, the HCPI model does not include samples of igneous rocks, suggesting potential for future enhancement by collecting similar lithologies to expand the model's applicability. In addition, both models have determination coefficients (R²) above 0.7, indicating a certain degree of explanatory power. The HCPI model performs better than the NHCB2 model, possibly due to the smaller sample size. The analysis samples for the two models differ by approximately three times, and the determination coefficient decreases with fewer samples. However, a larger sample size signifies a more statistically significant model. Future efforts can continue to collect landslide samples to build a more representative model while noting changes in the determination coefficient values with increasing sample size.

Finally, Figure 5 reveals differences in the estimated curve trends between the two integrated indices (HCPI and NHCB2) and K, with more pronounced differences when K is smaller. The NHCB2 model estimates higher values of permeability coefficient K, suggesting higher overall hydraulic conductivity in the major mountainous watersheds of Taiwan. In contrast, the hydraulic conductivity values estimated by the HCPI model are generally lower, implying a geological environment in landslide areas characterized by significant disturbance and more clay shear in geological strata, resulting in lower trends of hydraulic conductivity. Therefore, based on the comparison results, while both models can estimate the hydraulic conductivity of fractured rock masses using geological characteristic factors, the choice of model should consider whether the strata have been disturbed to obtain the required hydraulic characteristics of the site correctly.

4.6. The Relationship between Hydraulic Conductivity Distribution Patterns and the Locations Of Sliding Surfaces

Using the HCPI model, this study produced continuous distribution charts of hydraulic conductivity with depth per meter for each borehole equipped with an inclinometer in the 24 landslide sites. The results showed that most landslide sites exhibited certain invertals where the hydraulic conductivity values significantly increased or decreased with depth. Figure 6 illustrates the hydraulic conductivity profile of the Wanda landslide site, indicating significant variations in hydraulic conductivity values at depths above and below 34 meters (as highlighted in red boxes).

According to previous studies on deep-seated landslides [50,51,52], the upper and lower layers above the sliding surface often exhibit a trend where the upper layer is more permeable (K > 10^-6 m/s,[53]) than the lower layer. The difference in hydraulic conductivity between the upper and lower layers typically falls within one to three orders of magnitude [50,51,52].

To verify whether the landslides collected in this study also exhibit similar trends reported in previous studies, cross-referencing can be done using existing inclinometer data showing clear deformation and hydraulic conductivity profiles. Reviewing the inclinometer monitoring data from the 24 landslide sites collected in this study, it was found that 14 landslide sites already have clear records of sliding in inclinometers. These sites include Caoling, Wanda, Yixing, Songmao, Xinjiayang, Longhua, Hongcaiping, Wushe, Huayuan, Baolong, Lantai, Kanjiao, Taipingshan, and Dayuling landslide sites. Therefore, the inclinometer data from the 14 sites were overlaid with the hydraulic conductivity distribution charts estimated by the HCPI model. Upon comparison, it was observed that for 12 inclinometer boreholes, there were clear interfaces of significant changes in hydraulic conductivity values corresponding to sliding depths. Figure 7 exemplifies the comparison results for the Longhua landslide site. The area within the red box, showing significant differences in hydraulic conductivity values above and below the depth of 44 meters, coincides precisely with the sliding depth recorded by the inclinometer.

Moreover, the average hydraulic conductivity values in the upper segments were greater than 10^-6 m/s, indicating potential groundwater accumulation at the interface, leading to increased pore water pressure and decreased slope stability. However, for two inclinometer boreholes, the relationship between the predicted sliding depth and hydraulic conductivity variation estimated by the HCPI model was less clear, namely at the Yixing and Huayuan landslide sites. Further examination of the inclinometer data at these two sites revealed S-shaped displacement curves, as shown in Figure 8, which differed from the typical significant displacement at a certain depth. This suggests possible compression of the inclinometer due to thick layers of backfilled soil, thin layers of fractured rock, shear zones, or highly weathered rock layers, causing subsidence and negative frictional forces on the inclinometer, resulting in irregular S-shaped buckling phenomena. Additionally, irregularities in the S-shaped curve patterns may also arise from inadequate backfilling during inclinometer installation, as indicated by inconsistencies in multiple measurements of S-shaped curve patterns.

In summary, comparing the continuous distribution diagrams of hydraulic conductivity per meter with depth from boreholes at 14 landslide sites and inclinometer displacement data shows that the actual deformation at 12 landslide sites closely matches the estimated sliding results of the model. However, for two landslide sites, direct comparison was not feasible. Overall, the comparative results are generally consistent, verifying that the developed permeability estimation model can effectively predict potential sliding surfaces.

4.7. Establishment of Identification Rules for Sliding Surface Locations

To more accurately estimate the potential sliding surface locations in landslide areas, this study integrates hydraulic conductivity profile information with literature-recommended weak sections of core samples related to the formation of sliding surfaces [34,52,54,55] and groundwater level data [33,34]. This newly proposed approach led to the establishment of an evaluation rule for identifying potential sliding surfaces. The evaluation method and steps are shown in Figure 9. First, it is necessary to determine whether the hydraulic conductivity profile indicates the presence of a rock layer with an average hydraulic conductivity value greater than 10^-6 m/s and a thickness of more than 2 meters. If this condition is met, the process can proceed to the following evaluation step. The evaluation items can be divided into three categories, with each item's evaluation criteria described as follows:

• The magnitude of hydraulic conductivity variation:

If the upper half of a rock mass interval is more permeable while the lower half is less permeable, groundwater may accumulate in the more permeable upper half of the rock medium. The greater the difference in permeability between the upper and lower sections, the higher the likelihood that water will concentrate in the more permeable upper section, thereby adding more groundwater weight to the upper rock layer. This can reduce the shear strength of the rock layer and lead to rock mass sliding. Therefore, this concept can be used as a basis for weighting the likelihood of sliding in the strata. If the permeability difference between the upper and lower rock strata is between one to two orders of magnitude, the water accumulation capacity is low, thus a weighting score of one is given. If the difference is between two to three orders of magnitude, the water accumulation capacity is higher, thus a weighting score of two is given. When the difference is more than three orders of magnitude, the water accumulation capacity is the highest, thus a weighting score of three is given.

2.

Presence of weak rock mass materials:

The occurrence of deep-seated sliding is often related to the degree of weakness in the properties of the rock material [34,52,54,55]. This study summarized the material properties of sliding surfaces from the 24 landslide sites and found that the materials at the sliding surface typically include shear zones with gouge, highly weathered slate or schist, and water-sensitive shale layers. Therefore, if the core sample contains 'shear zones with gouge,' 'highly weathered slate or schist,' or 'shale layers,' each is given a weighting score of one point.

3.

Position of the groundwater table:

Groundwater can induce weak layers in the strata to become sliding surfaces [33,34]. Therefore, if the potential sliding surface is located below the normal groundwater table, it is assigned a weighting score of one point.

Based on the aforementioned scoring rules, this study further calculated the sliding potential score for each sliding position from 12 landslide sites with clear sliding records. This results in data on sliding surface locations, hydraulic conductivity values for the upper and lower rock mass sections, and scoring results for the 12 borehole drilling records (as shown in Table 10). As shown in Table 10, the sliding potential scores of the 12 existing sliding sections all exceed three points, and all contain weak layers in the core samples. Therefore, this study concludes that if the sliding potential score of a rock layer section exceeds 3 points, it should be considered a potential sliding surface section. According to this conclusion, if future sites without installed inclinometers or sites where inclinometers have not yet detected signs of sliding are encountered, the evaluation rules proposed in this study can be applied. If the total score of a certain rock layer section exceeds 3 points after calculation, that section is likely to induce a sliding surface in the future, thereby predicting the occurrence location of the sliding surface in advance.

Figure 9. Evaluation method and procedure for identification of potential sliding surfaces.

Figure 9. Evaluation method and procedure for identification of potential sliding surfaces.

4.8. Application of the Evaluation Approach - Prediction of Potential Sliding Surfaces

Based on the evaluation method developed in this study for quickly assessing the locations of potential sliding surfaces, this method can then be used to estimate the potential sliding surface depths for the remaining ten landslides where no sliding has yet been detected during the monitoring period (due to the short monitoring period). Through the calculation process shown in Figure 9, a total of 12 potential sliding sections were identified in 9 landslide sites. The depths of each potential sliding surface, hydraulic conductivity values for the upper and lower sections, and scoring conditions of each potential sliding surface in each landslide site are summarized in Table 11. The results show that the sliding surface estimation method designed in this study can identify relatively less stable sliding surfaces in nine of the ten landslide sites, and these sliding depths are more likely to develop into sliding surfaces in the future.

Based on the proposed identification rules, Antong is the only landslide site where no potential sliding sections were identified. According to the rock core photos and core records of the landslide site, it was found that Antong lacks easily sliding-prone strata materials, such as shear mud zones and highly weathered strata. Therefore, when applying the identification method established in this study, the scores for strata materials could not be obtained, resulting in sliding potential scores for the sliding sections being below 3 points. This leads to the determination, based on the collected data, that there is no potential sliding risk. Consequently, the geological conditions of the landslide site are relatively more stable than the other nine landslide sites.

5. Conclusions

This study proposes a preliminary and efficient assessment method to predict the potential slip surface depth in landslide areas using geological characteristic data and the developed rock mass permeability estimation models, thereby reducing investigation time and economic costs. It also seeks to facilitate early preparation and formulation of remediation strategies in disaster prevention efforts, ultimately reducing the occurrence probability of slope disasters.The principal findings can be concluded as follows.

• By collecting hydrogeological data from the 24 active landslide sites across Taiwan, a permeability estimation model (HCPI model) was established using six geological characteristic factors (RQD, DI, GCD, LPI, FA, and FD) related to rock permeability with the multicollinearity analysis, principal component analysis, and regression analysis. The coefficient of determination (R²) of the HPCI model is 0.895, meeting the practical application accuracy requirements. In addition, the HCPI model suits various lithological geological conditions and can be used to estimate the permeability variation per meter with depth along a borehole.
• The HCPI model was compared with the non-disturbed rock mass permeability estimation model "NHCB2". Under the same geological characteristic conditions, the hydraulic conductivity estimated by the HCPI model is smaller than that of the NHCB2 model. This outcome can be attributed to the highly disturbed geological conditions of the landslide sites, and the overall geological environment contains more gouge content, resulting in smaller hydraulic conductivity. Thus, choosing a geologically compatible estimation model is crucial while estimating hydraulic conductivity in the geological environment of mountainous areas.
• The concept of using abrupt changes in hydraulic conductivity data to estimate the depth of potential slip surfaces in landslide areas was feasible and confirmed by 12 landslide cases. According to the investigated landslide sites, the sliding surface occurred within a certain depth range of the strata, where the hydraulic conductivity of the upper layer was greater than that of the lower layer by one to three orders of magnitude. Additionally, the hydraulic conductivity of the upper layer is typically greater than 10^-6 m/s.
• An identification rule predicting sliding surface locations has been established through three criteria: the magnitude of hydraulic conductivity variation, the presence of weak rock mass materials, and the position of groundwater table. This identification rule has been applied to ten landslide sites to seek high-potential sliding zones, which can be used as valuable information for slope disaster prevention management.
• Finally, the HCPI model was constructed based on sedimentary and metamorphic rock data samples, but no samples from igneous rocks. Future work can include igneous rock data and continue collecting more landslide data for analysis, enhancing the model's reliability and practical application value.