Climate Change Impacts on Soil Erosion in a Tropical Watershed: Insights from SINGV Regional Climate Modelling and RUSLE

2.1 Study area

The study was conducted in the Sub-Biyonga watershed, a sub-basin of the Limboto River Basin located in Gorontalo Province, Indonesia. The watershed plays an important role in maintaining the hydrological balance of Lake Limboto, which has increasingly experienced environmental pressures due to land-use change and watershed degradation (Olii & Ichsan, 2020). Moreover, the Lake Limboto catchment area is included in the conservation priorities regulated under Presidential Regulation (Perpres) Number 60 of 2021 concerning the Rescue of National Priority Lakes (2021). Sub-Biyonga covers an area of approximately 72.07 km² (00°36′–00°44′ N, 122°56′–123°40′ E) and is characterized by a transition from lowland plains to hilly terrain. Gentle slopes (<15%) dominate the watershed landscape, while steeper slopes in the upstream areas are more susceptible to surface runoff and soil erosion processes (Jaya et al., 2024).

The region experiences a humid tropical climate, with an average annual rainfall of about 2,063 mm and mean temperatures ranging from 26–27 °C (Jaya et al., 2024). Ongoing land-use conversion from forest to agricultural land has reduced vegetation cover and increased runoff potential within the watershed (Olii & Ichsan, 2020). These environmental and geomorphological characteristics make the Sub-Biyonga watershed a suitable case study for assessing soil erosion dynamics and their potential responses to future climate scenarios. Figure 1 presents the study area, including catchment boundaries, drainage network, and elevation.

Figure 1. Location map of the Biyonga sub-watershed within the Limboto River Basin, Gorontalo Province, Indonesia.

Panels show (a) the location of Sulawesi Island, (b) the Limboto watershed, and (c) the Biyonga sub-watershed including elevation, river network, watershed boundary, and hydrometeorological stations.

2.2 Data sources

This study integrates multiple datasets to analyze soil erosion dynamics in the Biyonga watershed. The datasets include observed rainfall records, regional climate model (RCM) precipitation projections, remote sensing products, and topographic data used to derive the parameters of the Revised Universal Soil Loss Equation (RUSLE). Observed rainfall data were obtained from three rain gauge stations within the watershed: Biyonga, Dulamayo, and Hepuhulawa. These records were used as reference data to evaluate and bias-correct precipitation outputs from climate models. Future precipitation projections were derived from SINGV Regional Climate Model (SINGV-RCM) simulations within the CMIP6 framework (Prasanna et al., 2024), using six climate models: ACCESS-CM2, EC-Earth3, MIROC6, MPI-ESM 1–2-HR, NorESM2-MM, and UKESM1–0-LL. The projections were analyzed under two climate scenarios, SSP245 and SSP585, representing intermediate and high greenhouse gas emission pathways.

Spatial datasets were used to derive the RUSLE factors. Rainfall erosivity (R) was estimated from bias-corrected precipitation data, soil erodibility (K) was obtained from a global soil dataset (Panagos et al., 2014), and the slope length and steepness factor (LS) was derived from the FABDEM digital elevation model (30 m resolution) (Hawker et al., 2022). Land cover information was obtained from the ESRI Global Land Cover dataset, while vegetation conditions used to estimate the cover-management factor (C) were derived from Sentinel-2 NDVI. All spatial processing was performed using the Google Earth Engine platform. All datasets were harmonized to a consistent spatial resolution and clipped to the Biyonga watershed boundary to ensure compatibility in the soil erosion modeling process. A summary of the datasets used in this study is provided in Table 1.

2.3 Bias correction and evaluation

Regional climate model (RCM) precipitation outputs often contain systematic biases relative to observed rainfall due to limitations in model resolution, parameterization of atmospheric processes, and the representation of regional climate variability. To improve the reliability of precipitation projections for hydrological applications, a bias correction procedure was applied to the RCM outputs prior to their use in soil erosion modeling (Cannon et al., 2015).

In this study, a hybrid bias correction framework was adopted by combining quantile mapping (QM) for the historical period and the quantile delta method (QDM) for future projections. This approach allows the correction of systematic biases in modeled precipitation while preserving the climate change signal contained in the projected datasets (Gupta et al., 2019).

To preserve the climate change signal in future projections, the quantile delta method (QDM) was applied. Unlike conventional quantile mapping, which may distort projected climate signals, QDM adjusts modeled precipitation while maintaining the relative change between historical and future simulations (Xavier et al., 2022). The bias-corrected future precipitation can be expressed as:

For rainfall conditions with sufficient wet-day samples, a gamma distribution was applied to represent precipitation distributions (Gupta et al., 2019). The probability density function of the gamma distribution is given by:

The performance of the bias-corrected precipitation datasets was evaluated using three statistical metrics: the Pearson correlation coefficient (r), root mean square error (RMSE), and the Kling–Gupta efficiency (KGE). In addition to these statistical metrics, model performance was visually assessed using Taylor diagrams, which summarize the agreement between simulated and observed precipitation in terms of correlation, standard deviation, and centered root-mean-square difference (Izzaddin et al., 2024). The bias-corrected precipitation datasets obtained from this procedure were subsequently used to estimate rainfall erosivity and to drive soil erosion simulations using the RUSLE model.

2.4 RUSLE model

Soil erosion in the Biyonga watershed was estimated using the Revised Universal Soil Loss Equation (RUSLE), which is widely applied for assessing long-term average annual soil loss caused by water erosion (Hagras, 2023). The RUSLE model integrates rainfall characteristics, soil properties, topography, land cover, and conservation practices to estimate soil erosion risk (Renard et al., 1997). The model is expressed as:

2.4.1 Rainfall erosivity (R factor)

The rainfall erosivity factor (R) was derived from annual precipitation data obtained from the bias-corrected regional climate model outputs. Following commonly applied empirical relationships for tropical regions used by Lihawa et al. (2025) in Tomini Bay watershed, the rainfall erosivity factor was estimated using the following equation established by Babu et al. (2004):

2.4.2 Soil erodibility (K factor)

The soil erodibility factor (K) represents the susceptibility of soil particles to detachment and transport by rainfall and runoff. In this study, the K factor was derived using soil physical properties obtained from the SoilGrids global soil database provided by ISRIC (Hengl et al., 2017). Soil properties including clay, sand, silt, soil organic carbon, coarse fragment volume, and bulk density were extracted for the upper soil layer (0–30 cm depth) by averaging the values from three depth intervals (0–5 cm, 5–15 cm, and 15–30 cm).

Soil organic carbon values were converted to soil organic matter using a standard conversion factor (Huang et al., 2009). In addition, soil bulk density was used to estimate packing density, which represents the degree of soil compaction (Jones et al., 2003). These soil properties were subsequently used to determine soil texture classes (Hengl et al., 2017) and soil structure classes (Bagarello et al., 2009) based on empirical classification rules commonly applied in soil erosion studies (Panagos et al., 2014).

To account for soil hydraulic behavior, soil permeability was estimated using saturated hydraulic conductivity $K_{\mathit{sat}}$ obtained from the HiHydroSoil v2.0 dataset (Simons et al., 2020). The permeability class was determined based on hydraulic conductivity thresholds combined with hydrologic soil group information.

The soil textural factor M, which represents the influence of soil particle size distribution on erodibility, was calculated as:

Finally, the soil erodibility factor was calculated using the following empirical formulation:

2.4.3 Topographic factor (LS factor)

The topographic factor (LS) represents the combined effects of slope length and slope steepness on soil erosion. In this study, the LS factor was derived from a 30 m resolution digital elevation model (DEM) obtained from the FABDEM dataset (Hawker et al., 2022). From the DEM, slope gradient, slope angle, and slope aspect were calculated to characterize terrain conditions within the watershed.

Slope length was estimated using the upstream contributing area derived from the MERIT Hydro dataset (Yamazaki et al., 2019), which represents the accumulated flow contributing to each grid cell. This contributing area was used to approximate the effective slope length influencing runoff concentration. The slope length component (L) was calculated using the following equation:

The final LS factor was then calculated as:

This approach allows the representation of terrain-driven erosion processes by incorporating both slope length and slope steepness effects derived from high-resolution topographic data.

2.4.4 Cover management (C factor)

The cover management factor (C) represents the protective effect of vegetation cover against soil erosion. In this study, the C factor was estimated from the Normalized Difference Vegetation Index (NDVI) derived from Sentinel-2 imagery (Almagro et al., 2019). NDVI was calculated using the following equation:

2.4.5 Support practice (P factor)

The support practice factor (P) represents the ratio of soil loss with a given conservation practice to the corresponding soil loss under conditions without conservation measures (Renard et al., 1997). In large watershed-scale studies, detailed information on soil conservation practices is often unavailable and difficult to map. Consequently, many studies assign a uniform P value of 1 to represent the absence of conservation practices (Gashaw et al., 2019; Liu et al., 2020).

In this study, the P factor was estimated based on a combination of land use/land cover and slope conditions to reflect potential conservation effects associated with vegetation cover and terrain characteristics. Cropland areas were assigned different P values depending on slope classes, considering that the effectiveness of conservation practices decreases with increasing slope. This approach follows the concept proposed by Wischmeier and Smith (1978), where P values vary according to slope conditions. Non-agricultural land cover types such as built-up areas and bare land were assigned a value of 1.0, representing the absence of effective conservation practices.

All spatial datasets were processed using the Google Earth Engine (GEE) cloud-based geospatial analysis platform to generate the input layers required for the Revised Universal Soil Loss Equation (RUSLE) model. The GEE environment enables large-scale geospatial processing and efficient integration of multi-source datasets such as satellite imagery, digital elevation models, and global soil databases. Data pre-processing involved filtering, reprojection, spatial harmonization, and clipping to the watershed boundary to ensure consistency in spatial resolution and coordinate systems. The derived factor layers were then combined using the multiplicative RUSLE framework to estimate spatially distributed annual soil loss across the watershed.

2.5 Validation RUSLE model

Model validation is an important step in soil erosion studies to ensure that the estimated results are reliable and applicable. Ideally, validation should be performed by comparing model outputs with field-based measurements obtained from erosion plots or monitoring stations. However, such ground-based erosion data were not available for the Biyonga Watershed. Consequently, an alternative validation strategy was adopted.

In this study, the model results were evaluated by comparing the estimated mean annual soil loss with values reported in previous peer-reviewed studies conducted in tropical regions with similar agro-ecological and topographic characteristics. This comparative validation approach has been widely applied in soil erosion modelling studies facing similar data limitations (Lihawa et al., 2025; Olii et al., 2023), and it provides a reasonable basis for assessing model performance when direct field observations are unavailable.

2.6 Future scenario analysis

Future soil erosion dynamics in the Biyonga watershed were assessed by integrating bias-corrected regional climate model (RCM) precipitation projections with the RUSLE framework. Climate projections were derived from the SINGV-RCM ensemble (Prasanna et al., 2024), driven by six CMIP6 global climate models (ACCESS-CM2, EC-Earth3, MIROC6, MPI-ESM 1–2-HR, NorESM2-MM, and UKESM1–0-LL) under two emission scenarios, SSP245 and SSP585.

Future precipitation datasets were generated for two projection periods, 2040 and 2060, representing mid-century climate conditions. The bias-corrected precipitation was then used to estimate the rainfall erosivity factor (R) using an empirical relationship between precipitation and rainfall erosivity. The soil erodibility (K), slope length and steepness (LS), cover management (C), and support practice (P) were assumed to remain constant due to the absence of reliable projections for soil properties, terrain characteristics, and land management practices. Consequently, changes in projected soil erosion are primarily attributed to variations in rainfall erosivity.

Annual soil erosion was estimated using the RUSLE model for each climate scenario. To quantify climate change impacts, projected soil erosion rates were compared with baseline conditions using the relative change:

Figure 2. Conceptual framework and analytical workflow integrating regional climate modeling, rainfall erosivity estimation, and RUSLE-based soil erosion assessment.

3.1 Bias correction of RCM precipitation

The performance of bias-corrected regional climate model (RCM) precipitation was evaluated using the Pearson correlation coefficient (r), root mean square error (RMSE), and Kling–Gupta efficiency (KGE) at three rain gauge stations: Biyonga, Dulamayo, and Hepuhulawa ( Table 2). Overall, the corrected models showed moderate skill in reproducing monthly rainfall variability, with correlation values generally below 0.4.

At the Biyonga station, ACCESS-CM2 exhibited the best performance (r = 0.36; RMSE = 127.07 mm; KGE = 0.356). At Dulamayo, the MME median produced the highest correlation (r = 0.329), while UKESM1–0-LL achieved the highest KGE (0.303). At Hepuhulawa, ACCESS-CM2 again showed the highest correlation and KGE (r = 0.232; KGE = 0.225), whereas UKESM1–0-LL yielded the lowest RMSE (90.87 mm). In contrast, MIROC6 and EC-Earth3 generally showed weaker performance across stations.

The comparative performance of the models is illustrated in the Taylor diagram ( Figure 3), which summarizes the agreement between model simulations and observations in terms of correlation, variability, and centered RMSE. ACCESS-CM2 and UKESM1–0-LL lie closest to the observational reference point, while the MME median provides relatively balanced performance across stations. These results indicate that bias correction improves the statistical consistency of RCM precipitation, although some limitations remain in reproducing rainfall variability. The corrected precipitation datasets were therefore used for subsequent climate impact analyses.

Figure 3. Taylor diagrams showing the performance of individual regional climate models and multi-model ensemble (MME) statistics in reproducing monthly precipitation at the Biyonga, Dulamayo, and Hepuhulawa stations.

The diagrams summarize model agreement with observations in terms of correlation coefficient, standard deviation, and centered RMSE.

3.2 Spatial distribution of baseline soil erosion

The spatial distribution of baseline soil erosion in the Biyonga watershed was estimated using the RUSLE model ( Figure 4). The results indicate pronounced spatial variability across the watershed, which reflects the heterogeneous nature of its biophysical characteristics. This variability is primarily influenced by differences in topography, rainfall erosivity, soil properties, and land cover conditions, all of which interact to control erosion intensity.

Figure 4. Spatial distribution of the RUSLE factors and estimated historical soil erosion in the Biyonga watershed.

Panels show (a) rainfall erosivity factor (R), (b) soil erodibility factor (K), (c) slope length and steepness factor (LS), (d) cover management factor (C), (e) conservation practice factor (P), and (f ) estimated historical soil loss (t ha⁻¹ yr⁻¹).

Estimated soil erosion rates range from 0.08 to 410.66 t ha⁻¹ yr⁻¹, with a mean of 40.20 t ha⁻¹ yr⁻¹ and a median of 34.34 t ha⁻¹ yr⁻¹ ( Table 3). These values indicate that the watershed generally experiences moderate erosion, although localized areas exhibit substantially higher erosion intensity. The relatively high standard deviation (31.22 t ha⁻¹ yr⁻¹) further highlights the heterogeneous nature of erosion processes within the watershed. Spatially, higher erosion rates are concentrated in areas with steeper slopes and higher LS factor values, which reach up to 45.49. In contrast, areas characterized by gentler slopes and denser vegetation cover tend to exhibit lower erosion rates.

Among the RUSLE factors, rainfall erosivity exhibits relatively high values across the watershed, with a mean of 853.22 MJ mm ha⁻¹ h⁻¹ yr⁻¹, indicating the strong influence of rainfall in driving erosion processes. In contrast, the soil erodibility factor shows moderate susceptibility (K = 0.0219), suggesting that soil properties contribute to erosion but are not the primary controlling factor. Meanwhile, the mean C factor (0.375) and P factor (0.646) reflect varying levels of vegetation cover and conservation practices, which locally modulate erosion intensity across the watershed. Taken together, these results indicate that while land cover and management practices play an important role, the overall erosion pattern is predominantly governed by the interaction between topography and rainfall erosivity. This dominance explains the spatial concentration of higher erosion rates in areas with steep slopes and intense rainfall exposure. Therefore, these baseline conditions provide a critical reference for assessing future changes in soil erosion under different climate scenarios, enabling a clearer interpretation of how shifts in rainfall patterns may alter erosion dynamics over time.

3.3 Validation RUSLE model

Direct validation of soil erosion models ideally requires field measurements from erosion plots or sediment monitoring stations. However, such data are rarely available in large tropical watersheds due to logistical and monitoring constraints. Similar limitations have been reported in watershed-scale erosion studies in Gorontalo (Olii et al., 2023). Therefore, a comparative validation approach was applied by comparing the estimated mean soil erosion with results reported in previous studies conducted within the same regional watershed system draining into Tomini Bay (Lihawa et al., 2025). The Biyonga watershed analysed in this study represents one of the sub-watersheds contributing to the Tomini Bay drainage system.

A study conducted in the Tomini Bay area reported average soil erosion rates ranging from 56.06 to 57.68 t ha⁻¹ yr⁻¹ using the RUSLE model. The baseline erosion estimated in the present study (40.88 t ha⁻¹ yr⁻¹) falls within the typical range reported for tropical watersheds in Sulawesi. Differences between the two estimates may arise from variations in model inputs and parameterization. The Tomini Bay study derived rainfall erosivity from CHIRPS precipitation data, estimated the soil erodibility factor (K) based on soil texture classification, and used land cover data representing conditions in 2020. In contrast, the present study incorporated updated land cover data for 2024 and different parameterization of RUSLE factors. Despite these differences, the magnitude of estimated soil erosion remains comparable with regional studies, indicating that the model outputs provide a reasonable representation of erosion conditions in the Biyonga watershed.

3.4 Projected soil erosion under climate scenarios

Projected soil erosion provides insight into how future climate variability may alter watershed-scale erosion dynamics. Future projections indicate relatively stable soil erosion patterns in the Sub-Biyonga watershed across climate scenarios. Mean soil erosion is projected to range from 46.33 to 49.03 t ha⁻¹ yr⁻¹, with the lowest value under SSP585 (2040) and the highest under SSP585 (2060). Median values remain similar (39.63–41.35 t ha⁻¹ yr⁻¹), indicating limited changes in the overall distribution of erosion intensity. Although minimum erosion values remain very low (0.09–0.10 t ha⁻¹ yr⁻¹), maximum values exceed 500 t ha⁻¹ yr⁻¹, reflecting localized areas with extremely high erosion potential ( Table 4).

Across all scenarios, Class II erosion remains dominant, covering approximately 54–56% of the watershed ( Table 5). Class III erosion represents the second largest share (26–28%), followed by Class I erosion (16–17%). In contrast, Class IV erosion occupies less than 1.2% of the watershed, while Class V erosion is absent in all projections. Overall, these results indicate that future climate scenarios may slightly increase erosion magnitude but do not substantially alter the overall spatial pattern, with moderate erosion continuing to dominate the watershed ( Figure 5).

Figure 5. Spatial distribution of projected soil erosion classes (t ha⁻¹ yr⁻¹) under SSP245 and SSP585 climate scenarios for 2040 and 2060 in the Biyonga watershed.

Soil erosion was classified into five severity classes based on RUSLE thresholds: very slight (0–15), slight (15–60), moderate (60–180), severe (180–480), and very severe (>480 t ha⁻¹ yr⁻¹).

3.5 Spatial hotspots of erosion increase

The analysis of soil erosion change (Δ) reveals a clear tendency toward increasing erosion across the Sub-Biyonga watershed under future climate scenarios. Mean soil erosion is projected to increase by 6.12–8.83 t ha⁻¹ yr⁻¹, with the largest increase occurring under SSP585 in 2060 ( Table 6). Median changes range from 5.26 to 6.37 t ha⁻¹ yr⁻¹, indicating that a large portion of the watershed experiences moderate increases in erosion rates. Although the minimum change remains negligible (0.01 t ha⁻¹ yr⁻¹), maximum increases reach 124.50 t ha⁻¹ yr⁻¹, highlighting the emergence of localized erosion hotspots.

Spatially, more than half of the watershed is projected to experience increasing erosion across all scenarios ( Figure 6). The proportion of areas with increasing erosion ranges from 52.71% to 61.27%, while the remaining 38.73–47.29% remains relatively stable. No areas show a decrease in erosion. The largest expansion of erosion-prone areas occurs under SSP585 in 2060, where approximately 61% of the watershed is projected to experience intensified soil loss.

Figure 6. Spatial patterns of projected changes in soil erosion (Δ soil loss) under future climate scenarios in the Biyonga watershed.

Positive values indicate increasing soil erosion relative to the baseline period, while negative values represent areas with decreasing erosion.

These results indicate that future climate variability may substantially amplify erosion risk across the watershed. The concentration of erosion increases in specific areas highlights the development of spatial erosion hotspots, particularly under high-emission scenarios. These emerging hotspots represent priority areas for targeted soil conservation under future climate change.

4.1 Climate-Driven changes in rainfall erosivity

Rainfall erosivity represents one of the most important drivers of soil erosion because it reflects the capacity of rainfall to detach and transport soil particles through raindrop impact and runoff generation. In the Biyonga watershed, the relatively high mean rainfall erosivity value (R = 853.22 MJ mm ha⁻¹ h⁻¹ yr⁻¹) indicates that precipitation plays a dominant role in controlling erosion processes. Such values are typical of humid tropical environments where intense rainfall events and convective storms frequently generate strong surface runoff and sediment transport (Almeida et al., 2021; Yan et al., 2024). Previous studies have shown that soil erosion is highly sensitive to rainfall intensity, with short-duration, high-intensity rainfall events often producing disproportionately large sediment yields compared with lower-intensity rainfall events.

Future climate projections suggest that rainfall dynamics may further intensify erosion processes in the watershed. Although the overall spatial distribution of erosion remains broadly consistent across climate scenarios, the projected increase in mean soil erosion rates ranging from 6.12 to 8.83 t ha⁻¹ yr⁻¹ relative to baseline conditions indicates that future climatic conditions may enhance rainfall erosivity. Moderate increase in erosion amidst noticeable changes in climate imply that erosion is controlled by non-climatic factors such as topography and land cover. This could be the reason for the existence of stable spatial patterns of erosion across scenarios. The largest increases occur under the SSP585 scenario toward the mid-century period, suggesting that higher greenhouse gas emission pathways may amplify the erosive impact of rainfall. Similar findings have been reported in global and regional studies showing that climate change can intensify rainfall erosivity and increase soil erosion risk due to stronger precipitation extremes and altered rainfall regimes (Borrelli et al., 2020; Eekhout & de Vente, 2022).

In tropical environments, erosion processes are particularly sensitive to extreme rainfall events. A relatively small number of high-intensity storms can contribute disproportionately to annual rainfall erosivity values, thereby increasing erosion risks even when total annual rainfall changes only moderately. In the Biyonga watershed, this may apply to localized erosion hotspots rather than across the entire landscape. Several studies have demonstrated that intensified precipitation associated with climate change can significantly increase runoff generation and sediment yield in many regions worldwide (Adeyeri et al., 2024; Yin et al., 2025). The results of this study support this mechanism, indicating that projected climatic changes may strengthen erosive forces in the watershed even without substantial increases in total precipitation.

The projected increase in soil erosion in the Biyonga watershed is also consistent with findings from other tropical and monsoon-dominated regions. Regional climate modeling studies in Southeast Asia indicate that precipitation extremes are likely to intensify under future climate scenarios, potentially increasing soil erosion risks in humid tropical watersheds (Hariadi et al., 2024; Nguyen et al., 2022). These findings align with global assessments indicating that climate change may substantially increase rainfall-driven soil erosion due to intensification of the hydrological cycle and more frequent extreme precipitation events.

From a modeling perspective, the erosion patterns observed in this study are consistent with previous applications of the Revised Universal Soil Loss Equation (RUSLE) combined with climate projections. Studies integrating climate model outputs with RUSLE have demonstrated that increases in rainfall intensity often translate into higher rainfall erosivity values and subsequently greater soil loss under future climate scenarios (Chuenchum et al., 2020; Weng et al., 2023). In many cases, rainfall variability exerts a stronger influence on erosion rates than land-use changes, highlighting the critical role of precipitation dynamics in controlling soil erosion under climate change conditions (Dash & Maity, 2023).

In the Biyonga watershed, the influence of rainfall erosivity is further amplified by topographic controls. Areas with high LS factor values, reaching 45.49, experience stronger runoff acceleration and higher erosive energy during intense rainfall events. The interaction between rainfall intensity and terrain characteristics is widely recognized as a key driver of erosion susceptibility in tropical watersheds, particularly in landscapes where steep slopes coincide with limited vegetation cover (Bag et al., 2022; Tessema et al., 2020). This interaction likely explains the emergence of localized erosion hotspots under future climate scenarios. Additionally, the relatively modest rate of change suggests that topography and landcover highly influence erosion dynamics in the Biyonga watershed, and this may explain why projected changes remain moderate across the watershed. Climate change may amplify erosion rather than cause major alterations to existing erosion patterns.

4.2 Implications for watershed management

The projected erosion patterns provide important insights for watershed management strategies in the Biyonga watershed. Although moderate erosion (Class II) remains the dominant category across the watershed, covering approximately 54–56% of the total area, the projected increases in erosion rates and the emergence of localized hotspots indicate that certain areas may become increasingly vulnerable under future climate conditions. These hotspots, particularly those identified under the SSP585 (2060) scenario ought to be prioritized for targeted interventions since they are more sensitive to increased rainfall intensity. These findings suggest that watershed management should move beyond reactive mitigation approaches and instead adopt proactive, climate-informed soil conservation strategies.

Topographic conditions represent one of the primary controls on erosion variability within the watershed. Areas characterized by steep slopes and high LS factor values are particularly susceptible to accelerated soil loss because slope gradients increase runoff velocity and erosive energy. Numerous studies have demonstrated that terrain morphology strongly influences soil erosion susceptibility, particularly in tropical catchments where rainfall intensity is high (Li et al., 2024; Wen et al., 2023). Consequently, conservation interventions should prioritize high-risk zones where steep slopes coincide with exposed or sparsely vegetated land surfaces. In the Biyonga watershed, such a slope represents the upstream areas with high LS values as indicated in Figures 3 and 5, and therefore should be considered first for conservation planning.

Vegetation cover also plays a critical role in moderating erosion processes. The moderate average value of the C factor (0.375) suggests that vegetation provides partial protection across the watershed, although substantial spatial variability remains. Areas with dense vegetation generally exhibit lower erosion rates due to increased rainfall interception, root reinforcement, and improved soil structure. In contrast, exposed or sparsely vegetated areas are more vulnerable to soil loss. Previous studies have shown that vegetation cover can significantly reduce runoff and sediment transport by stabilizing soil aggregates and enhancing infiltration capacity (Desta et al., 2021; Wei et al., 2024).

The conservation practice factor (P) further indicates that soil conservation measures are present but not uniformly implemented across the watershed. Expanding the spatial coverage and effectiveness of conservation practices could therefore play a crucial role in mitigating future erosion risks. Practices such as contour farming, terracing, agroforestry systems, and vegetative buffer strips have been widely shown to reduce soil loss and improve land productivity by enhancing soil moisture retention and reducing runoff velocity (Jiang et al., 2024; Mihretie et al., 2022).

The relevance of these findings extends beyond the Biyonga watershed to other tropical catchments in Southeast Asia, where rapid land-use change and climate variability frequently interact to increase erosion susceptibility. Agricultural expansion and deforestation across many tropical landscapes have increased the exposure of soil surfaces to rainfall impact, thereby amplifying erosion risks during extreme precipitation events. Integrated watershed management approaches that combine climate adaptation strategies with sustainable land-use planning are therefore essential for reducing soil degradation and enhancing watershed resilience under changing climate conditions. Since erosion increases were more prevalent under climate scenario SSP585, designing conservation measures ought to be guided with the high intensity rainfall scenarios rather than rely on baseline practices.

4.3 Limitations and recommendations

Despite the insights provided by the projected erosion scenarios, several sources of uncertainty should be considered when interpreting the results. One major source of uncertainty arises from the climate projections used in this study. Global climate models often have limited ability to represent local rainfall processes, particularly in tropical regions where convective and orographic precipitation dominate (Nguyen et al., 2022). In addition, the effectiveness of bias correction is constrained by the limited availability of long-term observational rainfall data, and identifying a single model that best represents regional climate characteristics remains challenging. It is noticeable that the correlation coefficient between modelled and observable rainfall was low (r < 0.4), implying that using climate models may not reproduce local rainfall, hence introducing uncertainty. This uncertainty in estimating rainfall erosivity might in the end impact on the accuracy of the soil erosion results. To reduce these uncertainties, this study applied statistical downscaling and bias correction and employed a multi-model ensemble (MME) approach to improve the robustness of climate projections. Nevertheless, future studies could further enhance projection reliability by incorporating longer observational datasets and evaluating additional climate models to better capture regional rainfall variability.

Uncertainty also arises from the soil erosion modeling approach itself. The Revised Universal Soil Loss Equation (RUSLE) is widely used to estimate long-term average soil erosion and has been widely applied in climate change impact assessments (Chuenchum et al., 2020). However, the model simplifies several complex erosion processes and does not explicitly represent gully erosion, sediment deposition, or event-based hydrological dynamics. Consequently, the model may underestimate erosion in areas where channelized erosion or extreme rainfall events dominate sediment transport processes. Future studies could address this limitation by integrating RUSLE with more process-based hydrological or sediment transport models, or by incorporating event-scale rainfall and runoff data to better capture erosion dynamics during extreme climate events.

In addition, several input parameters used in this study, such as land cover and conservation practices, were assumed to remain constant under future climate scenarios. In reality, land-use change driven by urban expansion, agricultural intensification, or conservation interventions may significantly alter erosion dynamics over time. Therefore, future research should incorporate dynamic land-use change scenarios using land change models or spatial prediction approaches, which would allow a more comprehensive assessment of the combined impacts of climate change and land-use transitions on soil erosion.

Finally, the spatial resolution of input datasets may introduce additional uncertainty in erosion estimates. The use of gridded climate datasets, satellite-derived vegetation indices, and digital elevation models may limit the representation of small-scale terrain features and local land management practices that influence erosion processes. Future studies could reduce this uncertainty by integrating higher-resolution datasets, such as LiDAR-derived terrain models, detailed soil surveys, or field-based erosion measurements, which would improve model calibration and enhance the reliability of spatial erosion assessments. Despite these limitations, the modeling framework provides a valuable first-order assessment of climate-driven soil erosion risk and offers an important baseline for future watershed management and adaptation planning in Sulawesi Island.

4.4 Novelty and contribution of the study

This study offers one of the first applications of bias-corrected high resolution SINGV regional climate model projections with spatially distributed RUSLE modeling to assess future erosion dynamics at the watershed scale. While many previous studies have evaluated soil erosion under climate change scenarios, relatively few have combined bias-corrected regional climate projections with spatial erosion modeling in tropical Southeast Asian environments.

Second, this study highlights the spatial heterogeneity of climate-driven erosion responses. The results demonstrate that future erosion changes are not uniformly distributed across the watershed but instead become concentrated in localized hotspot areas where rainfall intensity interacts with terrain characteristics. This finding emphasizes the importance of spatially explicit erosion assessments rather than relying solely on watershed-scale averages.

Finally, the study demonstrates the value of integrating climate projections with geospatial erosion modeling to support adaptive watershed management strategies. By identifying areas that may become increasingly vulnerable to erosion under future climate conditions, the results provide a scientific basis for targeted soil conservation and climate adaptation planning in tropical watersheds.