Machine Learning for Rainfall-Driven Debris Flow Prediction in Data-Scarce Volcanic Watersheds

Debris flows are highly destructive mass movements in mountainous and hilly terrain, causing severe damage to infrastructure, ecosystems, and human life [1]. These rapid, gravity-driven flows consist of water, soil, rock fragments, and organic matter and often occur suddenly with little warning [2]. Their destructive potential arises from both high kinetic energy and unpredictable occurrence [3]. Frequent debris flow events across Asia, Europe, and South America result in fatalities, infrastructure disruption, and substantial economic losses. Increasing human exposure in hazard-prone mountain valleys, combined with more intense rainfall linked to climate change, has heightened the need for reliable and timely debris flow prediction systems [4].

Traditional debris flow prediction methods commonly rely on empirical rainfall thresholds, statistical analyses, and physically based models [5]. Although rainfall thresholds are useful in some regions, they are highly site-specific, sensitive to local hydrological conditions, and often require long-term data for calibration [6]. The complex interaction of rainfall intensity, antecedent soil moisture, geology, and slope morphology further limits the ability of single parameter approaches to represent debris flow initiation. Physically based models require detailed geotechnical data and high computational effort, restricting their applicability for real-time hazard management [7]. These limitations underscore the need for data-driven approaches capable of learning nonlinear relationships directly from observed data.

In recent years, machine learning (ML) has become an effective tool for hazard prediction and environmental monitoring. By leveraging historical data, these models can identify complex patterns that are difficult to capture with traditional methods and can handle nonlinear, high-dimensional, and noisy inputs [8-10]. In debris flow research, ML enables the integration of multiple predictors, such as rainfall intensity, duration, and timing, into a unified predictive framework. This data-driven approach represents a shift from rigid threshold-based systems toward more flexible models capable of capturing complex hazard dynamics [11-13]. ML-based debris flow prediction models can support disaster risk reduction strategies by strengthening early warning and decision-making frameworks, particularly when aligned with national and regional disaster risk reduction policies. In addition, data-driven hazard prediction frameworks complement adaptation and mitigation strategies for hydro-meteorological disasters, contributing to improved resilience in hazard-prone regions [14].

Despite growing interest in ML for hazard prediction, several challenges remain. Debris flow datasets are typically small, limiting the effectiveness of complex models that require large training samples [15]. In addition, class imbalance, where hazardous events are much rarer than non-events, reduces model sensitivity for minority categories [16]. Practical usability is another key challenge: beyond predictive accuracy, model interpretability and operational integration are essential for disaster management and rapid decision-making [17]. These limitations highlight the need for approaches that balance accuracy, robustness, and usability.

Accordingly, this study aims to develop a reliable and practical framework for predicting debris flow magnitude using rainfall and timing variables. The objectives are threefold: (i) to evaluate the predictive value of rainfall and temporal features, (ii) to establish a robust model suitable for small and imbalanced datasets, and (iii) to translate the selected model into an operational decision-support tool for early warning and risk assessment. Through these objectives, the study contributes both to scientific understanding and to actionable hazard mitigation in debris-flow-prone regions.

The methodological framework adopted in this study is summarized in Figure 1. The process begins with data collection of rainfall and timing variables, followed by data preprocessing to ensure completeness and consistency. An exploratory data analysis (EDA) was then conducted to examine feature distributions, class imbalance, and correlations among predictors. In the model development stage, eight ML classifiers were trained and optimized using cross-validation. Their performance was assessed during the model evaluation stage using Accuracy and macro-averaged F1 score (F1-macro) metrics. To enhance interpretability, model interpretation was carried out through feature importance analysis and simplified decision tree (DT) visualization.

Figure 1. Workflow of the proposed debris flow prediction framework

3.1 Data sources

The Putih River Watershed, situated on the southwestern flank of Mt. Merapi in Magelang Regency, Central Java, represents one of the most dynamic lahar-prone systems in Indonesia [29]. Following the 2010 eruption, which deposited substantial volumes of unconsolidated pyroclastic material across the upper basin, the watershed experienced pronounced alterations in its hydrological and sedimentological regime. This geomorphic disturbance continues to enhance the sensitivity of the channel network to intense monsoonal rainfall, resulting in frequent debris flow events [30]. The watershed extends across volcanic slopes, agricultural land, and densely vegetated areas, as illustrated in Figure 2, forming a complex landscape where natural processes and human activities interact directly.

Figure 2. Location map of the Putih River Watershed in Magelang Regency, Central Java, Indonesia [10]

Figure 3. Field observations conducted in the midstream agricultural zone of the Putih River

Post-eruption sediment supply remains high, sustaining debris flow hazards and channel instability [31]. To characterize current conditions, field surveys were conducted in upstream and midstream reaches of the Putih River, focusing on active sediment pathways in agricultural areas (Figure 3). Measurements included channel geometry, flow depth, bank erosion, and cross-sectional morphology, while sediment samples were collected to characterize grain-size variability.

These field observations were integrated with DEMNAS topography, long-term rainfall records, and historical debris flow data from hydrological stations. A dataset of 33 debris flow events, described by 11 rainfall-related predictors and classified into small, medium, and large magnitudes, was compiled. Together, these datasets provide a robust foundation for ML analysis of rainfall-driven debris flow magnitude in the Putih River Watershed.

3.2 Data description and preprocessing

The dataset comprises 33 debris flow events described by 11 rainfall-related predictors derived from the Agromulyo and Ngepos stations. Variables include station-specific and mean rainfall totals, rainfall start and end times, debris flow onset and termination, and rainfall duration. This limited station coverage reflects typical monitoring constraints in Indonesian volcanic watersheds and represents a realistic low-data setting. Debris flow magnitude is classified into three categories: small, medium, and large.

Before model training, data quality was verified, and no missing values were identified. Rainfall predictors were standardized for distance-based models, while tree-based models used unscaled inputs. The response variable was numerically encoded (small = 0, medium = 1, large = 2), and stratified cross-validation was applied to preserve class proportions and address imbalance.

Exploratory analysis (Figure 4) indicates strong right-skewness in rainfall variables. Agromulyo rainfall ranges from 0.5 to 113.9 mm (mean: 31.0 mm), Ngepos from 0 to 124 mm (mean: 25.0 mm), and average rainfall from 0.25 to 89.45 mm (mean: 28.5 mm). Time-related variables show more uniform distributions, with rainfall typically initiating in the afternoon and debris flow following shortly thereafter. Despite the limited network, these stations provide the most reliable high-resolution rainfall data near debris flow initiation zones, supporting development of lightweight and operational prediction models in data-scarce environments.

Given the limited dataset size (33 events), a repeated stratified 5-fold cross-validation strategy was adopted to reduce overfitting and statistical randomness, with the full evaluation repeated 50 times using different fold partitions. This approach provides a more robust estimate of model generalization under small-sample conditions, with performance summarized using the mean and standard deviation of Accuracy and F1-macro. To address pronounced class imbalance, particularly the limited number of medium debris flow events (n = 5), imbalance-aware learning was incorporated through stratified sampling and cost-sensitive class weighting. All preprocessing and imbalance-handling procedures were implemented strictly within the cross-validation framework to prevent information leakage, and class-specific sensitivity was computed from pooled out-of-fold predictions to ensure stability under extreme data scarcity.

The class distribution of the response variable is shown in Figure 5. Small debris flows are the most frequent (15 events), followed by big events (13), while medium debris flows are underrepresented (5 events). This imbalance poses a risk of model bias toward majority classes and is therefore explicitly addressed during model training.

Boxplots of key predictors (Figure 6) show that rainfall variables are the strongest discriminators of debris flow magnitude. Average rainfall clearly separates small, medium, and big events, with the highest values associated with big debris flows; similar patterns are observed at the Agromulyo and Ngepos stations. Event duration also contributes to discrimination, as medium and big events generally last longer than small ones.

Figure 4. Histograms of predictor variables

Figure 5.  Class distribution of the target variable indicating imbalance among categories

Figure 6. Boxplots of top predictive features across debris flow classes

The correlation heatmap (Figure 7) reveals strong positive correlations among rainfall variables, particularly between station rainfall and average rainfall (r ≈ 0.80). Timing variables are also highly correlated, with start and end times strongly linked to rainfall stop indicators (r > 0.8), indicating consistent temporal behavior across events.

Figure 7. Correlation heatmap of predictors

These patterns indicate that rainfall and timing variables form coherent feature groups, which may introduce multicollinearity in linear models. Tree-based algorithms are less sensitive to such dependencies and can effectively exploit correlated predictors. Accordingly, all preprocessing steps, including scaling, encoding, and class weighting, were applied independently within each cross-validation fold to ensure unbiased evaluation.

The exploratory analysis highlights three key insights: (i) rainfall magnitudes are strongly associated with larger debris flow events, (ii) the dataset is imbalanced, particularly for medium-magnitude events, and (iii) several predictors are highly correlated within rainfall and timing groups. These findings informed subsequent preprocessing and model selection, including stratified sampling and appropriate handling of correlated features.

3.3 Model development for debris flow prediction

The process of developing predictive models for debris flow classification was structured into three stages. In the first stage, eight ML algorithms were evaluated to establish a baseline of performance. These models were carefully selected to represent a wide methodological spectrum, including tree-based classifiers, ensemble learners, kernel-driven methods, and distance-based approaches. A description of each model is provided below.

3.3.1 Decision tree

DTs are hierarchical models that recursively split the input space into regions defined by decision rules. At each internal node, the model selects a feature and threshold that best separates the target classes based on impurity measures such as the Gini index or entropy [32]. For a dataset $D~$with classes $c$, entropy is given by Eq. (1):

$H(D)=-\sum_{c=1}^C p_c \log _2\left(p_c\right)$         (1)

where, ${{p}_{c}}~$is the proportion of class $c$. The DT grows until a stopping criterion is reached, after which pruning may be applied to reduce overfitting. The general structure of a DT used in this study is illustrated in Figure 8, showing the progression from the root node to successive splits and final leaf nodes that represent prediction outcomes.

Figure 8. Structure of a decision tree (DT) model

3.3.2 Random forest (Bagging)

RFs are ensemble methods that combine multiple DTs trained on bootstrapped subsets of the data, with random feature selection at each split. The final prediction is obtained by majority voting across trees. The bagging process reduces variance and improves generalization [33]. If $T~$trees are trained, the ensemble prediction is given by Eq. (2):

$\hat{y}=\text{mode}\left\{ {{h}_{t}}\left( x \right),t=1,2,\ldots ,T \right\}$           (2)

where, ${{h}_{t}}\left( x \right)$ is the prediction of the $t$-th tree. The conceptual structure of an RF voting mechanism is illustrated in Figure 9, showing how multiple trees independently evaluate an input before aggregation.

Figure 9. Example of tree-level decision paths used within a random forest ensemble

3.3.3 Subspace k-nearest neighbors

This method is a variant of k-nearest neighbor (kNN) that uses random feature subspaces for distance calculations, improving robustness and reducing the influence of redundant predictors. For a query point $x$, the class is assigned based on the majority label among its $k$ nearest neighbors under a chosen distance metric (Euclidean) [34]. The prediction is given by Eq. (3):

$d\left( {{x}_{i}},{{x}_{j}} \right)=\sqrt{\underset{m=1}{\overset{M}{\mathop \sum }}\,{{({{x}_{im}}-{{x}_{jm}})}^{2}}}$            (3)

where, $M$ is the dimension of the feature subspace.

3.3.4 K-nearest neighbor, k = 5, standardized

The kNN algorithm is a non-parametric classifier that assigns a class label to a new sample based on the majority class among its $k$closest training instances in the feature space. Distance is typically computed using Euclidean metrics, making the algorithm sensitive to differences in feature scale. To ensure fair distance comparisons, all predictors in this study were standardized to a zero mean and unit variance.

Figure 10. K-nearest neighbor (kNN) classification mechanism before and after applying kNN

A value of $k=5$ was selected to balance bias–variance trade-offs, providing smoother decision boundaries while preventing overfitting associated with very small $k$. An illustration of the kNN classification process is shown in Figure 10, demonstrating how a new data point is assigned to a class based on the nearest labelled neighbours.

3.3.5 Logistic regression (ECOC)

Logistic regression models the posterior probability of class membership through the logistic function [35]. For binary classification, the probability of class $y=1$ is given by Eq. (4):

$P\left( y=1\mid x \right)=\frac{1}{1+{{e}^{-\left( {{\beta }_{0}}+{{\beta }^{T}}x \right)}}}$         (4)

For multi-class classification, this study employs the error-correcting output codes (ECOC) framework, which decomposes the problem into multiple binary logistic regression models. Each classifier corresponds to a column of the ECOC coding matrix, and the final class assignment is determined by matching binary outputs to class codewords. This coding strategy improves robustness through redundancy. An illustrative overview of the logistic regression process is shown in Figure 11.

Figure 11. Logistic regression classification model

Figure 12. Workflow of the random under-sampling boosting (RUSBoost) algorithm

3.3.6 Random under-sampling boosting

Random under-sampling boosting (RUSBoost) is an ensemble method designed for imbalanced datasets. It combines boosting, which iteratively reweights misclassified samples, with random under-sampling (RUS) of the majority class to reduce class imbalance [36]. For an ensemble of $T$ weak learners, the final hypothesis is a weighted combination of the individual classifier outputs, expressed as Eq. (5):

$H\left( x \right)=\text{sign}\left( \underset{t=1}{\overset{T}{\mathop \sum }}\,{{\alpha }_{t}}\,{{h}_{t}}\left( x \right) \right)$            (5)

where, ${{\alpha }_{t}}~$denotes the weight assigned to the $t$-th weak learner ${{h}_{t}}\left( x \right)$, typically a DT with limited depth. A schematic overview of the RUSBoost workflow showing RUS sampling, iterative reweighting, weak learner construction, and final ensemble prediction is presented in Figure 12.

3.3.7 Support vector machine with RBF kernel (SVM-RBF, ECOC)

SVMs are margin-based classifiers that separate classes by maximizing the margin between support vectors [37]. The RBF kernel allows nonlinear decision boundaries to be given by Eq. (6). In this study, ECOC was used to extend binary SVMs to multi-class prediction.

$K\left( {{x}_{i}},{{x}_{j}} \right)=\text{exp}\left( -\gamma \parallel {{x}_{i}}-{{x}_{j}}{{\parallel }^{2}} \right)$         (6)

3.3.8 Support vector machine with linear kernel (SVM-linear, ECOC)

A linear SVM attempts to find the hyperplane that maximizes the margin between classes [38]. The equation is given by Eq. (7).

$\text{f}\left( \text{x} \right)=\text{sign}\left( {{\text{w}}^{\text{T}}}\text{x}+\text{b} \right)$           (7)

where, $w$ represents the weight vector normal to the hyperplane and $b$ is the bias term. The maximization of the margin between support vectors enhances generalization, particularly when the classes are linearly separable.

For multi-class classification, the ECOC strategy was employed. ECOC decomposes the multi-class problem into multiple binary SVM classifiers, each corresponding to a column of the encoding matrix. The final class prediction is determined by selecting the class whose codeword has the minimum decoding loss relative to the set of binary outputs. An illustration of the kernel mapping concept showing how non-linear class distributions can become linearly separable in a transformed feature space is provided in Figure 13.

Figure 13. Visualization of the kernel trick

ML classifiers were selected to balance predictive performance, robustness to small datasets, and interpretability. Tree-based models capture nonlinear rainfall thresholds, ensemble methods (RF, RUSBoost) assess variance reduction and imbalance handling, and distance- and kernel-based models provide benchmark comparisons under identical data conditions.

All models were optimized using grid-search tuning within cross-validation to reduce overfitting. DTs were tuned for depth (3–10), RF for trees (50–200), kNN for neighbors (k = 3–9), logistic regression for C (0.01–10), SVM for kernels (linear, RBF) with C = 0.1–100 and γ = 0.001–1, and RUSBoost for learners (50–200) and learning rate (0.01–1). Optimal models were selected using cross-validated F1-macro scores.

3.4 Performance evaluation of the developed models

Evaluating ML models requires reliable and widely adopted performance metrics, particularly for imbalanced classification problems [39, 40]. In this study, overall Accuracy (Acc) and the F1-macro were selected as the primary evaluation metrics to provide a systematic and comparable assessment of model performance.

Accuracy measures the proportion of correctly classified samples and offers an intuitive indicator of overall predictive success. However, in imbalanced datasets, Accuracy can be biased toward majority classes and may not adequately reflect performance on minority events [41]. To address this limitation, F1-macro was employed, as it computes precision and recall independently for each class and then averages them, ensuring equal weighting of minority and majority debris flow categories.

Beyond threshold-based metrics, discriminative capability was assessed using the macro-averaged AUC-ROC in a one-vs-rest framework, while probabilistic reliability was evaluated through calibration curves comparing predicted probabilities with observed frequencies. Computational efficiency was also examined by averaging training and prediction times across cross-validation folds. Together, Accuracy and F1-macro provide a robust evaluation framework, capturing both overall correctness and balanced performance across debris flow magnitudes [42], with their mathematical definitions presented in Eqs. (8) and (9):

$\text{Accuracy}=\frac{\text{TP}+\text{TN}}{\text{TP}+\text{TN}+\text{FP}+\text{FN}}$           (8)

$\mathrm{F} 1-\mathrm{macro}=\frac{1}{\mathrm{C}} \sum_{\mathrm{i}=1}^{\mathrm{C}} \frac{2 \cdot \text {Precision}_{\mathrm{i}} \cdot \text {Recall}_{\mathrm{i}}}{\text {Precision}_{\mathrm{i}}+\text {Recall}_{\mathrm{i}}}$        (9)

This study has several limitations. The dataset is small (33 events) and highly imbalanced, particularly for medium debris flows, which may limit generalization and sensitivity for intermediate hazard levels. Feature engineering is deliberately simple and based only on rainfall magnitude and timing from two stations, excluding other influential factors such as terrain properties, soil moisture, antecedent rainfall, and land cover. In addition, model performance may be affected by data drift and climate change, as future rainfall–debris flow relationships may differ from historical patterns. Finally, the model was developed for a single watershed, and its transferability to other regions with different geomorphological and climatic conditions remains untested. Future work should therefore focus on expanding event inventories and integrating multi-source data to improve robustness and generalization.

The authors gratefully acknowledge Universitas Muhammadiyah Yogyakarta and RKI-PRN Equity Research Grant 2025/2026 for providing institutional support and research facilities that enabled this study.