Trenggalek Regency is one of the regions in East Java Province, Indonesia. Trenggalek is located at 111° 23' 31.2" to 111° 50' 48.95" E,and 8° 23' 18.46 " to 7° 53' 18.16' S. Trenggalek Regency experiences hydrometeorological drought disasters that occur annually. There were 123 drought-related disasters in Trenggalek Regency throughout 2023, highlighting the waterscarcity problem. The population growth rate of Trenggalek Regency was 0.51% in 2023, in line with the increasing dependence on groundwater as a freshwater source, which reached 135,293 m3. Therefore, a groundwater potential map is urgently needed, as it can help in optimizing the fulfillment of clean water needs in Trenggalek Regency.

Trenggalek Regency has diverse geomorphological conditions. Two-thirds of Trenggalek Regency is mountainous. Trenggalek Regency has elevations ranging from 0 to 1,179 metres above sea level. Trenggalek Regency is a tropical region with an average rainfall of 1,879 to 2,756 mm/year and an average land surface temperature of 24°C to 31°C. The geology of Trenggalek Regency is dominated by the Mandalika formation, which is Neogene in age and consists of interbedded volcanic breccia, lava, tuff, and sandstone.

The selection of conditioning factors can affect the quality of groundwater potential modelling results. The conditioning factors were determined based on a literature review and data availability. This study considered 18 parameters such as TWI, river density (RD), SPI, land cover (LC), normalized difference vegetation index (NDVI), annual rainfall (AR), evapotranspiration (ET), land surface temperature (LST), atmospheric pressure in the surface (SP), elevation, slope, slope aspect, plan curvature (PC), lineament density (LD), geological unit (GU), soil type (ST), complete bouguer anomaly (CBA), and average shear wave velocity to 30m (VS30). The 18 parameters are classified into six factors: hydrology, land cover, climate, topography, geology, and geophysics. Detailed conditioning factor information can be seen in Table 1.

There were 1,840 samples used to train the model, consisting of 740 potential points and 740 non-potential points. The possible points used are spring locations with discharge rates of 10–100 liters/second. Meanwhile, non-potential points were random points located in GU with low porosity (0.2% - 8%). The distribution of training points can be seen in Figure 1. All parameters were resampled in raster format with a spatial resolution of 30m. The details of the data used in this study can be seen in Table 3. The data was processed with spatial data processing software to produce modelling parameters. The parameters can be shown in Figure 2 and Figure 3.

This study uses two types of data: discrete and continuous. To ensure data reliability in the modelling process, pre-processing was carried out. Discrete data are transformed into numerical form using the one-hot encoding (OHE) technique. OHE transforms each class into a binary vector, with 1 if it appears in a particular row and 0 if it does not. The OHE technique can be seen in equation 1 (Samuels, 2024). Meanwhile, for numerical data, the Yeo-Johnson transformer (YJT) and standard scaler are applied. The YJT is a technique that aims to transform the data distribution to near normal (Yeo & Johnson, 2000). On the other hand, the standard scaler or z-score standardization technique is a method to ensure model independence from specific feature scales. The standard scaler transforms data to a mean of 0 and a standard deviation of 1 (Demir & Sahin, 2024). The YJT and standard scaler can be seen in Equations 2 and 3.

(1)

(2)

(3)

Where, Is OHE vector in parameter-I, class j; is the result of the YJT transform; y is original value of the parameter; is original value of the parameter; is a constant estimated by maximum likelihood estimation (MLE); Z is a standardized value; x is the original value before the standardized value; is the mean value; and is the standard deviation ((Demir & Sahin, 2024); (Samuels, 2024); (Yeo & Johnson, 2000)).

This study evaluates the use of ML and DL for groundwater potential mapping in Trenggalek Regency, Indonesia. In addition to using each algorithm individually, this study used the stacking technique to integrate ML and DL models. Each algorithm performed hyperparameter tuning using a tree-structured Parzen estimator (TPE). The two ML algorithms used are RF and GBDT. RF is an ML algorithm that uses bootstrapping and bagging to form multiple DTs. Each DT is constructed with several n-parameters, and n-data are randomly selected. The result of the RF is the voting of the results of each single DT (Breiman, 2001). GBDT is a straightforward combination of DT and gradient descent. GBDT comprises two main processes: boosting and classification. These two processes can improve model performance and minimize overfitting through the regularization process within sequential steps to iterate on the loss function (Chen et al., 2020). Loss function is the difference of prediction and the actual value (Yang et al., 2024).

CNN and RNN were used as DL algorithms in this study. CNN can solve complex problems and consists of three main layers: convolutional layers (CL), pooling layers (PL), and fully connected layers (FCL) (Zhang et al., 2022). CL is a convolutional filter (also known as a kernel) where the input features expressed by a matrix of dimension N×1×P are convolved with a filter of size n×1×p, where N>n and P>p, resulting in an output feature map. PL is a subsampling technique applied to the feature map generated by CL. After the PL layer, the feature map will be flattened and moved to the FCL layer. In this layer, each neuron is connected to all neurons from the previous hidden layer, so it is called a fully connected (FC) approach (Fang et al., 2020). The overall architecture of the CNN used in this study is illustrated in Figure 4(a). RNNs are designed to process data sequentially by considering hidden states that capture information about previous inputs. RNNs have recurrent connections (Ji et al., 2023), as shown in Figure 4(b). Such connections allow information to loop within the network. The network memorizes the input information of the previous state (t-1) and then applies it to the calculation of the output at state t. The nodes between hidden layers are interconnected, and the input to the hidden layer includes not only the production of the previous input layer but also the output of the hidden layer from the previous state (Mienye et al., 2024).

This study uses two strategies to optimize the prediction model: hyperparameter optimization for each algorithm and stacking. Stacking is a method of successfully integrating different ML or DL models (Sikora, 2017). In the stacking model, heterogeneous prediction models are used, where the output of all the algorithms is given as input to the meta-model to get a more stable model (Khoshkroodi et al., 2024). Logistic regression (LR) is the meta-learner applied in this study, which utilizes the maximum likelihood method to discover the 'best fit' association between input and output (Sun et al., 2018).

TPE is the technique used to tune each algorithm's hyperparameters. TPE is a Bayesian technique used for optimizing hyperparameters that utilizes a kernel density estimator to model the probability distribution of hyperparameters using an objective function (Ishii et al., 2023). TPE imitates the density distribution function to optimize hyperparameters in both the optimal and worst-case settings. The iterative approach is used to select hyperparameters for sampling, evaluate model performance, and update probability density functions. The conditional probability theory p( x|y ) is used, where x is the hyperparameter and y is the loss from using that hyperparameter. The first step in TPE optimisation is selecting the threshold loss (y^*) given the available data, such as based on the median. The formation of probability density functions and is supported by the threshold loss as seen in equation 4 (Rong et al., 2021). Afterwards, other hyperparameter combinations are determined by maximizing the likelihood between and using equation 5 (Islam et al., 2024).

y^*. \end{cases} \end{document} ]]> (4)

(5)

Model evaluation is a validation process for the accuracy and reliability of a prediction model. Evaluation in this study was carried out on two datasets, namely training and testing, to determine whether each model suffered from overfitting. Model evaluation is carried out using four parameters. Model evaluation parameters can be seen in Table 4. Where TP (true positive) is a positive class that is predicted positive; FP (false positive) is a negative class that is predicted positive; FN (false negative) is a positive class that is predicted negative, and TN (true negative) is a negative class that is predicted negative ((Stern, 2021); (Tharwat, 2018)).

Model evaluation was performed using k-fold cross-validation (KFCV). KFCV divides the dataset into n-fold. 1/n fold is used to evaluate the model, and the other n-1 fold is used to train the model. This research uses five folds. This means the training-to-testing ratio is 80:20. The KFCV technique is used to detect overfitting by using the t-test with degrees of freedom (df) n-1, which was four. Model inspection improves our understanding of how to interpret the model. Two additional approaches to model inspection in this study were permutation feature importance (PFI) and partial dependence plots (PDPs).

PFI approximates the relative importance of parameters in a predictive model by permuting each parameter's values. As accuracy varies significantly, the parameter has an important contribution to the model (Altmann et al., 2010). The PDP employs the Shapley additive explanation (SHAP) value, which computes the Shapley value for every sample in the dataset. Shapley values offer a game-theoretic mechanism for reasonably apportioning the prediction output to individual input features. SHAP was initiated by Shapley, (1953) with game theory. It considers every possible combination of parameters in a row of data to calculate the marginal contribution of each parameter. This value is calculated for each parameter value. The SHAP value can be calculated using equation 6. In this equation, ∅ i v (x ) is the SHAP value for feature i with the prediction model v(x); S is the feature subset dataset; F is the entire feature; S ⊆ F / \{i } is a subset of feature F without considering parameter-I; v ( x S ∪{i } ) is the prediction output with subset S considering feature-i; and v (X S ) is the prediction output with subset S without considering i-feature (Lundberg & Lee, 2017).

(6)

Hyperparameter tuning is one of the crucial steps of constructing an optimal prediction model (Jafari et al., 2023). TPE reduces computational cost by selecting the most appropriate pair using a Bayesian-based method, rather than evaluating all pairs of hyperparameters. The optimum hyperparameter configurations for each algorithm are presented in Table 5. In this study, the number of iterations is limited to 100 trials, with test-set accuracy as the objective value. Although there is no universal judgment on the number of trials that ensures maximum accuracy or convergence, several studies have used a setting of 100 trials and achieved convergent results ((Fan et al., 2022); (Kilic et al., 2024); (Zhao et al., 2025)).

Optimal hyperparameter optimization using TPE can be visualized in a line diagram, as shown in Figure 5. Figure 5 illustrates the relationship between trial n and the objective value. According to the figure, the RF, CNN, and GBDT algorithms have an accuracy of more than 0.9 at the initial state or trial to 0. Meanwhile, the RNN algorithm has an accuracy of 0.49 at the initial state. It happens because, in the initial state, the hyperparameter configuration is randomly selected. The best trial or convergence conditions for the RF, GBDT, CNN, and RNN algorithms are 4, 84, 25, and 86, respectively. The highest objective values in the RF, GBDT, CNN, and RNN algorithms are 0.952, 0.962, 0.962, and 0.945, respectively. Overall, GBDT and CNN have the highest objective values, but CNN finds convergent conditions earlier than GBDT. While it has a very low aim value at the initial state, the RNN algorithm experiences a drastic increase in the third trial.

Groundwater potential maps are generated by applying each prediction model to pre-processed parameter data, including categorical data encoding, transformation, and standardization for numerical data. There are seven prediction models generated in this study, namely RF, GBDT, RNN, 1D CNN, RF-GBDT stacking, 1D CNN-RNN stacking, and RF-GBDT-RNN-1D CNN stacking. The output of the modelling process is the probability of groundwater potential. The range of probability values is 0 to 1. The closer the value is to 1, the higher the groundwater potential. Meanwhile, the closer the value is to 0, the lower the groundwater potential. The groundwater potential probability is classified into five classes: very low (0 -0.2), low (0.2 -0.4), medium (0.4 -0.6), high (0.6 -0.8), and very high (0.8 -1). Equal intervals were used to demonstrate the consistency and interpretability of groundwater potential mapping results expressed as continuous probability values-equal intervals allowed for an objective comparison of groundwater potential class interpretations across prediction models.

Natural and quantile classification techniques are highly dependent on the distribution of each prediction model result, leading to less objective comparative analysis. A few spatial modelling studies have used the equal interval technique to facilitate comparisons among models ((Bi et al., 2025); (Hossen et al., 2025); (Nurwatik et al., 2022); (Prasad et al., 2020)). The groundwater potential map produced in this study can be seen in Figure 6. The map is visualized with graduated colors ranging from light blue to dark blue. Light blue represents very low potential, while dark blue represents very high potential. Overall, the probability distribution almost has the same pattern between prediction models, where high probability is found in the west of the study area. Meanwhile, low probabilities are found in the north, some in the southwest.

Figure 7 presents the percentage area of each potential class and the percentage of spring locations within each class. Each model indicates that the class with the largest area is the low class, with an area of 22.32% to 31.29%. The very low class is the smallest area, except for the RF model. The very low class has an area percentage ranging from 12.57% to 17.66%. For the RF model, medium is the smallest-area potential class, with 16.03%. The percentage of potential points ought to be related in a direct relationship to the groundwater potential class. The lowest potential class has the least percentage of points at 4.05%, and the highest percentage is in the very high class at 38.11%. Most models show an inverse trend between the low and medium classes, with more potential points in the class of low potential. However, in general all models showed that the class of very high potential contained the most significant proportion of potential points.

The predictive model that has been established was evaluated through the KFCV method and utilized the ACC, MCC, CK, and ROC-AUC assessment metrics. The evaluation results are represented through the average value and standard deviation. Standard deviation can be used to quantify the robustness of each metric across the five cross-validation iterations. The results are shown in Table 6. A low standard deviation for each metric indicates stable performance across all folding schemes. Overall, all models achieved high prediction accuracy (ACC > 0.9; MCC and CK > 0.85). The GBDT algorithm demonstrated better consistency and stability than other single algorithms, as evidenced by its lower performance variance. GBDT, which iterates on bias values, can produce more stable prediction models with the best performance. Although both are built on DT, GBDT is superior to RF because it minimizes prediction variance while placing greater emphasis on bias. However, the performance of GBDT and RF is not significantly different, with only a 0.01 difference on several parameters, and GBDT is relatively more stable than RF. The bagging technique, which randomly selects features, makes RF less stable than the boosting technique used in GBDT.

Stacking learning (SL-ML) models that integrate RF and GBDT achieve the best performance and generalization on the test data, thereby confirming the safe and beneficial combination of the two algorithms. The SL-ML technique outperformed because it complemented the weaknesses of each weak learner. The SL technique reduced the bias and variance of the prediction values produced by each weak learner. SL-ML had advantages over SL-DL and the combination of SL-ML-DL due to the nature of the data used in the groundwater potential prediction model. The DL algorithm outperforms ML when there is a highly complex relationship between dependent and independent variables. This study confirms that the relationship between dependent and independent variables in mapping groundwater potential in the study area is not too complex. DT-based algorithms such as RF and GBDT are better at capturing the relationship between dependent and independent variables in this case study.

In the DL group, RNN outperformed CNN. Although RNNs are designed to solve sequential problems such as time-series prediction, research confirms that RNNs also excel at nonsequential data. RNN builds internal hierarchical relationships among independent variables, enabling it to capture better relationships with dependent variables in the dataset used in this study. This technique outperforms CNN-based convolutional feature extraction in this study.

Model performance can be visualized using the ROC-AUC curve, which is shown in Figure 8. Several models show an AUC value of 1.00 (perfect) on training data, but decrease to 0.98–0.99 on testing data. These models are RF, GBDT, and SL-ML. An ROC-AUC value of 1.00 is standard in ensemble models based on decision trees such as RF, GBDT, and SL-ML. The same phenomenon has been observed in several studies using DT-based ensemble models ((Jin et al., 2024); (Ouali et al., 2023); (Prasad et al., 2025)). However, it should be noted that the decreases are not significant, with the resulting gap ranging from 0.02 to 0.01. This phenomenon shows that prediction models can learn the relationship between dependent and independent variables well enough to generalize new data (testing) effectively. The lowest AUC in the training dataset is for RNN, at 0.97. The highest AUC value in the testing dataset is the stack ML model (0.99). The other models have an AUC value of 0.98. In general, for every model included: DL, DL-ML, GBDT, RNN, CNN, and Stack ML, the AUC values from prediction models using the training and testing data are all more than 0.95, which indicates that all the models have excellent performance.

The t-test is performed to assess whether there is a statistically significant difference between the two evaluations on the training and testing data. This was done to determine whether the model was overfitting. The H0 set was that the two datasets were not statistically different at a 95% confidence level. The t-test result refers to the P − value (PV) which is the probability value to support H0. The H0 is acceptable if the PV is higher than 0.05. The test results for each model were represented in Table 6 with the PV component. Based on the PV value, it can be inferred that none of the models overfit because the PV value is greater than 0.05. In general, all the prediction models performed well and did not exhibit overfitting.

PFI analysis identifies the most influential parameters in each prediction model. The results of the PFI analysis were conducted on all four basic learning models, as can be seen in Figure 9 Error! Reference source not found.. The PFI values are presented in a bar chart, with the x-axis ranging from 0 to 1 and the y-axis showing the variable names, sorted in order of highest PFI values. According to these results, GU is the most influential parameter in all models with PFI values ranging from 0.242 to 0.349. Other parameters that significantly influence the four prediction models are CBA, RD, TWI, and slope. There are nine parameters with a PFI value of 0 in the RNN model: NDVI, ST, LC, AR, LST, slope aspect, VS30, LD, and PC. Meanwhile, the other three models each have one parameter with a PFI value of 0.

Parameters with PFI 0 in the CNN, RF, and GBDT models are ET, SP, and AR, respectively. Visually, other parameters have PFI 0 in the RF and RNN models. However, these parameters have values greater than 0 but less than 0.001. A value of 0 in PFI may be due to PFI's limitations in explaining agnostic models. PFI only calculates the decrease in accuracy when the value of a variable in the test data is randomized to assess its importance. A value of 0 in PFI does not represent that the variable has no theoretical effect on the prediction model. However, this condition occurs when randomizing a particular variable does not decrease the accuracy of the initial model. PFI is highly sensitive to variables that are correlated with each other (Kaneko, 2022). If two variables are highly correlated, randomizing one variable will not significantly affect the initial accuracy, because the other variable already captures the values required by the model. For instance, LST and NDVI have a correlation value of 0.4 according to the analysis results. This yields LST or NDVI PFI values of 0 or near 0. A feature selection process is needed to ensure that the variables are not redundant with one another when creating a prediction model.

To identify the response characteristics of each parameter to the prediction model, this study uses PDP. PDP describes the relationship between a parameter's value and its SHAP value. PDP analysis is only performed on the best basic learning model, namely GBDT. The analysis results can be seen in Figure 10 Error! Reference source not found.. Continuous parameters are represented by a line graph (dark blue line) with an uncertainty polygon (light blue color). The uncertainty polygon represents a confidence interval of 95%. Meanwhile, for discrete parameters, it is defined as the average SHAP value for each class, shown in a bar chart.

The non-linear characteristics of each variable in relation to groundwater potential can be interpreted through the PDP. Topographical features such as elevation and slope have a negative relationship with groundwater potential. This indicates that the higher the surface elevation, the lower the groundwater accumulation. The higher the surface, the faster the runoff process in the area, resulting in less infiltration. It is supported by the condition that steeper slopes can increase lateral flow and reduce the infiltration process. The south-east to north slope direction shows a positive SHAP value due to the shorter solar exposure compared to the north to east direction. The duration of solar exposure can affect the level of evapotranspiration, resulting in a decrease in groundwater quantity in the area. Concave terrain has a positive SHAP value due to its morphological ability to accumulate water properly.

Hydrological conditions such as SPI, TWI, and river density likely have a negative relationship with groundwater potential. SPI and TWI are hydrological parameters derived from topographical factors. Drainage connectivity is represented by river density. The higher the drainage connectivity and TWI, the greater the potential for surface flow accumulation, thereby increasing the probability of lateral infiltration phenomena. SPI represents flow strength; theoretically, the higher the flow strength, the lower the lateral infiltration. However, the PDP results suggest a different interpretation. The higher the SPI value, the higher the groundwater potential. This could occur when supported by permeable geological and soil conditions.

Climatological factors play a significant role in the capability of groundwater recharge and storage. AR with a range of 2,200–2,400 mm/year has a high SHAP value, whereas above this range, the SHAP value decreases significantly. This occurs because high rainfall typically occurs in areas with high topography and is often associated with steep slopes. This condition results in low infiltration rates despite high recharge potential due to high runoff processes. ET shows a negative correlation with SHAP values, meaning that higher ET is associated with greater potential for groundwater, surface water, and vegetation water loss. LST negatively correlates with SHAP, but SHAP values are highest when LST values are too. It is possible that in high-temperature areas, supported by high groundwater runoff (basin areas), with permeable geological conditions. SP negatively correlates with SHAP values, because higher atmospheric pressure can reduce porosity by narrowing the inter-rock gaps beneath.

Land cover and vegetation conditions influence the characteristics of groundwater dynamics. The NDVI parameter is positively associated with groundwater potential. This supports the statement that areas with high vegetation can increase the infiltration process. The type of land cover with the highest SHAP value is residential. Residential areas may reduce infiltration capacity. However, this can be supported by geological conditions with high permeability and recharge potential.

Subsurface conditions greatly influence groundwater transport systems. In the GU parameters, it has the highest range of SHAP values compared to other parameters. The negative SHAP value in GU is obtained by the GU class with low porosity, namely 0% - 8%. Meanwhile, positive SHAP values are obtained by the GU class with a high porosity of 25% to 55%. The soil type parameter in all classes has a negative SHAP value. The lowest negative value is Latosol, which has a low permeability value of less than 0.13 cm/hour. LD has a negative association with the SHAP value because sloping areas (where lienaments are rare) reduce the runoff rate, resulting in a positive effect on groundwater recharge.

Both CBA and VS30 geophysical properties are positively associated with groundwater potential. This condition indicates that the higher the rock density, the higher the groundwater potential in the study area. This may occur because volcanic influences dominate the study area. The Mandalika geological formation dominates the study area. This formation has rock types with medium to high density but is prone to fracturing. The presence of these fractures allows surface water to be transported into groundwater.

This study implemented AI-based algorithms for groundwater potential mapping. The utilization of AI-based algorithms can improve the accuracy and objectivity of groundwater potential mapping. RF and GBDT algorithms represent ML in this study. Meanwhile, RNN and CNN represent the DL algorithm. The integration of ML, DL, and ML-DL algorithms is carried out to optimize the prediction model with the SL approach with LR as meta-learning. Utilized 18 parameters, all prediction models performed well with ACC>0.90; CK>0.85; and MCC>0.85. An interesting finding on the performance of the prediction models was that by applying SL to integrate the algorithms, they can improve 2.04% to 4.45% in ML stacking, 0.44% to 2.32% in DL stacking, and 1.18% to 5.47% in ML-DL stacking compared to the performance of every single algorithm. The application of stacking techniques has been demonstrated to decrease the bias and variance of prediction outcomes significantly (Lu et al., 2023).

Another important finding is that the CBA parameter, which is rarely used in modelling groundwater potential, shows a significant influence on all models in this study. Values of CBA variation represent lateral density characteristics of lithological units (Zaenudin et al., 2020). The most influential variable in each model is GU. Rock characteristics such as porosity affect the absorption rate of surface water into groundwater. The GU variable, as the variable with the highest relative contribution, aligns with research by (Bai et al., 2022), which had a relative contribution value of 0.95. Research by (Nugroho et al., 2024), GU is in the second-highest position with a relative contribution rate of 0.16. However, other studies produced the opposite condition, where the GU parameter had less effect on the prediction model, for example, in a study by (Aslam et al., 2025), where GU was in the last position in the relative contribution level of the prediction model, with a value of 0.15. The less influential parameters in this study are mainly climatological factors; for instance, in the GBDT model, AR has the lowest PFI value, SP is the lowest PI in the RF model, and ET is the lowest PFI in the CNN model. This condition aligned with research by (Hasan et al., 2025) where climatological factors such as rainfall had the lowest relative contribution, at 0.89%. The rainfall variable in the study by (Prasad et al., 2020) was also lacking in significant influence in predicting groundwater potential with a contribution value of 10%. However, studies by (Aslam et al., 2025) and (Bai et al., 2022) had different results, where rainfall had a significant influence on the prediction model. This indicates that each case study has its own characteristics in indicating the potential presence of groundwater.

Compared to previous similar studies using a hybrid algorithm integration approach and AI for groundwater potential mapping ((Arabameri et al., 2020); (Chen et al., 2020); (Chen et al., 2022); (Parasar et al., 2025); (Pham et al., 2019)), this study has advantages in performance optimization through hyperparameter tuning and stacking techniques. Both techniques can optimize the performance of the prediction model, resulting in good and stable evaluation values during cross-validation. In addition, this study also uses modelling parameters that are rarely used in similar studies, namely CBA and VS30. Both parameters have been proven to contribute significantly to modelling using AI-based algorithms through PFI values. Modelling parameters used in this study not only utilize surface phenomena but also subsurface phenomena through parameters on geological and geophysical factors. Hence, the model is more reliable and applicable for decision-making purposes related to freshwater resource management.

This study can provide implications for depicting areas lacking access to groundwater potential in Trenggalek Regency. All models indicate that the study area is dominated by low potential. Some parts of the area, however, have high potential. This information can be leveraged for decision-making in providing equitable access to freshwater resources in Trenggalek Regency during such emergency conditions as drought. The primary contribution of this study is to explore the ML, DL, and SL integration algorithms for groundwater potential. Moreover, this study bridges the black-box problem of the prediction model by applying PFI and PDP. Therefore, it is important to make the model more interpretable. The limitation of this study is not conducting feature selection, which may improve the performance of the prediction model and increase its effectiveness due to reduced computational cost. Further exploration of AI-based algorithms and other parameters is recommended, as the best prediction model in this study may not necessarily perform the following when applied to other study areas due to differences in geological and geomorphological characteristics.