1. Introduction
Soil erosion by water is a naturally occurring process associated with the hydrologic cycle. At small spatial scales, soil particle detachment by rainfall impact may predominate, but at larger spatial scales other water erosion processes are likely to dominate. Erosion of soil from catchment areas and sediment deposition in waterways can reduce storage capacity in the reservoir and degrade downstream water quality. Increasing soil erosion can also decrease soil fertility and crop yield [1], severely threatening local, national, and global food production systems and environmental sustainability [2,3]. Soil erosion severely restricts agricultural land use by reducing the productive potential of soils. It also contributes to water pollution by introducing suspended matter and nutrients into bodies of water. In addition, land that has been eroded becomes susceptible to other environmental impacts. Soil erosion is primarily caused by unsustainable agricultural practices, forest clearance, overgrazing, mining, and construction activities [4].
Soil erosion is considered one of the worst environmental issues in the Philippines [2,5]. Many parts of the country are highly susceptible to soil erosion because of their steep topographic conditions, which are compounded by the occurrence of severe rainfall events, land cover degradation, poor farming practices, and other soil-related factors [2]. Reduction in reservoir storage capacity is significantly related to increased soil erosion in the catchment areas. This is true in the case of the Pulangi IV reservoir in Bukidnon, Philippines. It was reported by Deutsch et al. [6] that the Pulangi IV reservoir was silting up at a rate of about one meter per year at the dam and that sedimentation had reduced the reservoir capacity by approximately 50%. Moreover, this contributed to the premature deterioration of hydropower turbines and frequent power outages [6]. Bukidnon is a watershed–landlocked province and the catchment area of various rivers in Mindanao.
To some extent, soil erosion can be mitigated. Better land use planning can help reduce the long-term threat of soil erosion [7]. There is a need for a straightforward and practical approach to estimating and mapping soil erosion risk that uses readily available data to improve water and soil conservation initiatives [3]. Estimating the risk of soil loss and its spatial distribution is critical for a successful assessment of soil erosion. Erosion-induced soil loss can be calculated using a prediction model such as the Revised Universal Soil Loss Equation (RUSLE) [8]. RUSLE is an empirical model for estimating soil erosion and is practical for predicting soil loss on hillslopes [2,9]. The model has been widely adopted due to its simple and straightforward computational input requirements compared to other conceptual and process-based models [10]. Numerous studies on soil erosion risk assessment and soil loss prediction around the world have been conducted using RUSLE in conjunction with Geographic Information Systems (GIS) and remote sensing techniques [1,2,3,4,8,10,11,12,13,14]. GIS has great potential in soil erosion inventory for soil erosion modeling and erosion hazard assessment [1] because it facilitates the manipulation, integration, analysis, and display of large amounts of spatial data, and can provide spatial distribution information on erosion [15,16].
According to the climate projection of the Philippine Atmospheric, Geophysical, and Astronomical Services Administration (PAGASA) [17], Mindanao will generally have a decreasing rainfall trend in 2050 (2036–2065). In 2050, seasonal rainfall in Bukidnon is projected to decrease significantly under high- and medium-range emission scenarios, most notably during the December–January–February (DJF) and March–April–May (MAM) seasons. Overall, the projected annual rainfall will reduce as well. The equivalent annual rainfall change based on the projected seasonal rainfall change under high-range and medium-range emission scenarios were estimated at −1.74% and −8.32%, respectively. A more recent report about the Philippine’s climate extremes revealed that the projected climate will be drier across the country, with more severe conditions expected in Visayas and Mindanao. Bukidnon, in the early-future (2020–2039), mid-future (2046–2065), and late-future (2080–2099), is still projected to have drier years [18]. Therefore, it is critical to incorporate climate projection scenarios when assessing the future risk of soil erosion.
While previous studies have attempted to characterize the soil erosion characteristics in some watersheds in Bukidnon [5,19,20], no study on GIS-based soil erosion assessments in all Bukidnon watersheds exists in the published literature, particularly in a predictive mode based on recent land use data and on recently reported climate change scenarios by PAGASA. Adornado and Yoshida [21] have previously used the RUSLE model to assess soil erosion in Bukidnon. However, the land cover in Bukidnon has changed significantly over time, and the climatic conditions are also expected to change. Additionally, technological advancements have occurred in recent years, land cover maps are updated regularly, the DEM data has much higher resolution, and satellite precipitation data are already available. Thus, this study aims to generate a more updated GIS-based soil erosion risk assessment in the watersheds in Bukidnon, assess the spatial distribution of soil loss based on the land cover maps in 2010, 2015, and 2020, and predict soil loss under various rainfall scenarios based on recently reported climate change projections.
2. Materials and Methods
2.1. Study Area
The province of Bukidnon is a landlocked province in the middle of Mindanao Island, southern Philippines, and is geographically located between 7°18.42′–8°38.22′ N latitude and 124°15.6′–125°31.74′ E longitude, as shown in Figure 1. The topography of the province is predominantly hilly and mountainous, especially in the eastern portion, and the other two mountain ranges in the west, have an average elevation of 915 m, and a range of 22 to 2867 m above sea level (see Figure 2a). It has rolling uplands, deep canyons, valleys alternating with the low plains, and terrain characterized by deep ravines and dense forest mountains in several mountain ranges. Due to its relatively high elevation, Bukidnon remains relatively cool and moist throughout the year. Bukidnon has a developing agricultural-based economy and is primarily a producer of rice, corn, sugar, coffee, rubber, flower, fruits, vegetables, poultry, and livestock.
2.2. Datasets
For rainfall, the TRMM_3B42_Daily v7 was used in this study. These data are a daily accumulated precipitation product generated from the research-quality three-hourly Tropical Rainfall Measuring Mission (TRMM) Multi-Satellite Precipitation Analysis (TMPA) [22]. Accordingly, it is produced at the National Aeronautics and Space Administration (NASA) Goddard Earth Sciences Data and Information Services Center (GES DISC), as a value-added product. A simple summation of valid retrievals in a grid cell was applied for the day data. These rainfall data in millimeters are available in raster format with a resolution of 0.25° × 0.25° (
The Interferometric Synthetic Aperture Radar (IfSAR) Digital Elevation Model (DEM) with 5 × 5 m resolution was used to delineate the province’s major river watersheds. This was obtained from the National Mapping and Resource Information Authority (NAMRIA) of the Philippines. The DEM was projected into WGS 84/UTM zone 51 N. The elevation map of Bukidnon is shown in Figure 2a.
The soil map of Bukidnon was extracted from the soil map of the Philippines and was downloaded from the Geoportal Philippines website (
Three land cover maps for 2010, 2015, and 2020 were used in this study. The land cover maps of Bukidnon were extracted from the 2010, 2015, and 2020 land cover maps of the Philippines. The land cover maps and the administrative boundary map were downloaded via the Geoportal Philippines website (
2.3. RUSLE
The Revised Universal Soil Loss Equation, or the RUSLE model, is a well-known and widely used empirical model for estimating soil erosion [9]. It was developed using the five following factors to estimate the average annual soil loss (A): rainfall erosivity factor (R), soil erodibility factor (K), slope length and steepness factor (LS), cover management factor (C), and conservation practice factor (P). RUSLE is expressed as:
(1)
The rainfall erosivity factor (R) quantifies the kinetic energy impact of rainfall and predicts the rate and amount of runoff associated with that precipitation [9]. Originally, to determine the R-factor, rainfall intensity data is required [10,26]. Due to the unavailability of rainfall intensity records for Bukidnon, the equation for R-factor (Equation (2)) used by Salvacion [2] in Marinduque, Blanco and Nadaoka [27] in Laguna Lake Watershed, and Adornado et al. [28] in Quezon, Philippines was adopted. Accordingly, this simplified equation produced a result that was acceptable for tropical and humid subtropical climatic zones [9,28]. The equation is as follows:
(2)
where R is the rainfall erosivity factor (MJ mm/ha/h/y), and P is the average annual precipitation (mm). Three sets of R-factor maps were generated with three different sets of average annual precipitation maps. The three sets of the average annual precipitation maps represented the average annual precipitation of 1998–2009, 2003–2014, and 2008–2019. Each period consists of an equal number of years, with a seven-year overlap between subsequent periods. This was done because this study used three land cover maps at different periods, as previously mentioned. Furthermore, the calculation of the average annual precipitation using no less than ten years of precipitation data was based on the methodology of the previous studies [28,29,30,31].The soil erodibility factor (K) indicates the soil’s resistance to erosion caused by raindrop impact, as well as by the runoff generated from the rainfall. Soil erodibility is determined by geological and soil characteristics such as structure, texture, parent material, porosity, and organic matter content. Regardless of the soil concentration of sand and clay, the silt content is directly related to soil erodibility [9,26]. This study adapted the K-factor values from David [32] for each soil textural class, which was also adapted by Salvacion [2].
The LS-factor in RUSLE is the slope length (L) and steepness (S) factors combined to reflect the effect of regional topography on the rate of soil erosion. Cumulative runoff increases in both amount and rate as the slope lengthens. As the land slope increases, the runoff velocity increases proportionately, resulting in massive erosion [9,28]. The LS-factor was generated using Equation (4). Equation (4) is a method developed by Moore and Burch [33] and Moore and Wilson [34] and was used by Hrabalíková and Janeček [35] in which it was proven as one of the better options to calculate the LS-factor. According to Andreoli [11], this expression is appropriate for areas with complex topography, including plateaus, terraced ledges, and mountains, because it considers the convergence and divergence of the flows. Using the DEM, flow accumulation and slope were generated in Quantum GIS (QGIS) through the r.flow Geographic Resources Analysis Support System (GRASS) algorithm, and slope algorithm of the Geospatial Data Abstraction Library (GDAL), respectively, in the QGIS toolbox. The equations are as follows:
(3)
(4)
where As is the specific catchment area, FA is the flow accumulation in each grid cell and its value corresponds to the number of flowlines that traverse that grid cell, cell size is the resolution of the grid (for this study 5 m), the value of m is 0.4, n is equal to 1.3, and β is the slope angle [9,11,35].The C-factor values represent how cropping and management practices affect the erosion rate. It is inextricably linked to land use types and is a factor in reducing soil erosion vulnerability. It is defined as the ratio of soil loss from land under specific conditions to the equivalent loss from clean-tilled, continuous fallow. Essentially, vegetative cover prevents raindrops from colliding with the soil surface and dissipates the kinetic energy of rainfall before it reaches the soil surface, slowing down runoff, thus facilitating infiltration; hence, the amount and type of vegetation cover has a significant effect on soil loss. C-factor is directly related to the vegetation type, stage of growth, and percentage of cover [1,9]. In this study, the C-factor values from David [32] and Delgado and Canters [36] were used as a reference in assigning the C-factor for each land cover class of the three land cover maps.
The conservation or support practice factor (P) indicates the effects of implementations that decrease the rate and amount of runoff, thereby reducing the amount and rate of soil erosion. It indicates the proportion of soil loss caused by a particular support practice compared to soil loss caused by upward and downward slope, contour farming, and tillage. Primary support practices include strip cropping, contour farming, terracing, cross-slope cultivation, and grassed waterways. P-factor values are calculated as the ratio of soil loss caused by a specific support practice to soil loss caused by row farming in both upward and downward slope conditions [9]. The value of the P-factor ranges from 0 to 1, where a value close to 0 indicates good conservation practice while values approaching 1 indicate poor or no erosion control practice [37]. Since no records document the extent and adoption of conservation practices in the province, though some may have adopted them, a value of 1 was assigned for the P-factor for the entire province.
2.4. Annual Rainfall Change Scenario
Initially, a baseline average annual rainfall scenario was established before formulating annual rainfall change scenarios. Since the accumulated annual TRMM_3B42_Daily data used in this study was from 1998 to 2019, the same period was used in generating the baseline scenario condition. Thus, the average annual rainfall for the period 1998–2019 was used to generate the R-factor of the baseline scenario.
As previously mentioned, Bukidnon is expected to experience drier years in the future. In Bukidnon, the projected rainfall on the early-future (2020–2039), mid-future (2046–2065), and late-future (2080–2099), and under both the moderate emission (RCP4.5) and high emission (RCP8.5) scenarios, range from −3.70 to −12.61% from the baseline value [18]. On this basis, three annual rainfall scenarios were generated (Table 1) to assess the risk of soil erosion under the projected rainfall conditions. Rather than using the six scenarios proposed by DOST-PAGASA et al. [18] to represent each future category and emission scenario, only three scenarios were considered because all projected annual rainfall values are consistently lower than the baseline value. To generate the rainfall amount of each scenario, the amount of rainfall equivalent to the percent change in rainfall, as shown in Table 1, was subtracted from the baseline average annual rainfall scenario. The results in each scenario were used to obtain the corresponding R-factor.
3. Results
3.1. RUSLE Factor Distribution
GIS was used in this study to generate the RUSLE factors, calculate the average annual soil erosion rates, and produce maps showing the distribution of these factors, the soil erosion rates, and the soil erosion risk map in Bukidnon. Cell-by-cell calculations of the mean annual soil erosion rates [1] were conducted; thus, the RUSLE factors were prepared in raster format using QGIS. Figure 4, Figure 5 and Figure 6 show the map of the RUSLE factors used in calculating the annual soil erosion rates.
As previously mentioned, the TRMM_3B42_Daily accumulated annual precipitation data in raster format downloaded from GIOVANNI were used as the dataset to generate the average annual rainfall in the province. Since three land cover maps were used in the study, three maps for the average annual precipitation were also created. Each represents the average for 1998–2009, 2003–2014, and 2008–2019, respectively. Using the average annual rainfall from 1998–2019 as a reference, it was observed that most of the northern and western parts of the province in 1998–2009 were wetter while the eastern areas were drier. In addition, the average annual rainfall in 2003–2014 was generally higher in all areas of the province compared to other time periods. In 2008–2019, however, the province was predominantly wetter, with a few drier areas on the western side. The average annual precipitation maps were used to calculate the R-factor maps using Equation (2), thus resulting in three R-factor maps, as shown in Figure 4. Generally, values of the R-factor are much higher in the northeastern and eastern parts of the province, but lower values of the R-factor can be found in the northwestern part of the province, as shown in Figure 4.
Using Equations (3) and (4), the LS-factor of the province was obtained. The spatial distribution of the LS-factor is shown in Figure 5a. Approximately, 48.3% of the province’s LS-factors are below 5, followed by 22.98% within 5–10, 24.16% within 10–25, 4.12% within 25–50, and only approximately 0.54% above 50. There are certain areas in the province with LS-factors over 500, but they only represent approximately 0.00008% of the province’s total area. Higher LS-factor values can be found in the mountain ranges, in deep canyons, and in nearby river networks, whereas lower values can be found in the plains and perhaps the cropland areas in the province. Extremely high values of the LS-factor are expected, given that the range elevation in the province is relatively high, having a standard deviation of 445 m. Additionally, the slope in the province ranges from 0 to 77° with 16.08° as the average and 12.27° as standard deviation, and approximately 20.33% of the province is mountainous and extremely steep, with slopes exceeding 26.6° or 50%, resulting in some extremely high LS-factor. Most of the LS-factors over 50 are in areas with a slope above 50%. LS-factors over 50 were also observed by Salvacion [2] in Marinduque, Philippines. Using the same expression of the LS-factor used in this study, Andreoli [11] was able to obtain higher LS-factors over 350 and others over 2000 [38,39].
Figure 5b shows the soil erodibility factor (K) of the province, which ranges between 0.19 and 0.6. Lower K-factor values are more prevalent in the eastern side of the province, while higher values are generally located in the northern side and a few other areas. Most of the province has a lower K-factor, ranging between 0.19 and 0.3.
Three C-factor maps were generated based on the land cover maps of 2010, 2015, and 2020. As shown in Figure 6, the C-factor values range from 0 to 1. During C-factor classification based on the land cover classes, the 0 value was assigned for water bodies while 1 was set for barren and built-up areas. The rest of the land cover classes were reclassified based on data available from the existing literature [2,32,36,40]. As land cover changes over time, so does the C-factor distribution in the province.
3.2. Soil Loss
Figure 7 shows the spatial distribution of soil erosion rate in Bukidnon during 2010, 2015, and 2020. Because the province is a watershed area of the neighboring provinces, the Bukidnon watershed areas were clustered into seven watershed clusters based on the classification considered by Rola et al. [41], as shown previously in Figure 1. These watershed clusters include Tagoloan in the north; Agusan-Cugman, and Cagayan in the northwest; Upper Pulangi in the central and eastern side; Maridugao in the southwest; Lower Pulangi in the south; and Davao-Salug in the southeastern side of Bukidnon.
On average, the RUSLE model-predicted soil erosion rates in Bukidnon range between 312 to 363 t/ha/y based on the periods under consideration. The predicted values range from 0 to 114,275 t/ha/y. The predicted soil erosion rates may appear to be exceedingly high. However, a study conducted in Marinduque, Philippines, showed that the RUSLE predicted soil erosion rates is around 120 t/ha/y on average, and can vary from 0 to as high as 20,767 t/ha/y [2]. These high erosion rates were observed in an island province that has a drier climate and less complex topography than Bukidnon. Additionally, the plot experiment on soil erosion in a few selected areas in Bukidnon revealed that a month of 5 mm rainfall can cause up to 229 t/ha of soil loss on certain plots with a slope of 15%. Although soil deposition can also occur, reaching 400 t/ha in some plots [42], this could also mean that somewhere near the plots an accumulated soil loss equal to that amount is also possible. Thus, the RUSLE-predicted soil erosion rates obtained in this study are comparable with the empirical evidence obtained from previous studies in the Philippines. As illustrated in Figure 7, areas prone to soil erosion are primarily found along the buffer zones of the major mountain ranges in Bukidnon. These areas are frequently found in deep river canyons and rolling areas that were previously used for crop production. The predicted soil erosion rates in these areas are extremely high, exceeding 300 t/ha/y on average. It can be observed in Figure 7 that the areas with lower soil erosion rates are the plains and mountain ranges that have good vegetative cover, the protected areas of the province.
Figure 8 illustrates the mean value of the average annual soil erosion rate for each watershed cluster. The thin lines in the graph represent the magnitude of standard deviations on top of the mean soil erosion rates. The predicted soil loss in 2020 is lower than in 2015 in most areas of the watershed clusters, except for Tagoloan. This decline may be attributable to a wetter climate in 2003–2014, compared to 2008–2019, that was used in deriving the R-factor. The exception in Tagoloan may be attributed to the decrease in vegetative cover in 2020, especially in the forest areas, despite the climate in 2008–2019 being drier than in 2003–2014. The soil erosion rates in the Maridugao watershed cluster are slightly decreasing. In contrast, the Upper and Lower Pulangi, and the Davao-Salug watershed clusters have an increasing soil erosion rate from 2010 to 2020, of which the highest can be observed during 2015. The increase in soil erosion in the Upper Pulangi watershed cluster may have contributed more to the accumulation of sediments in the Pulangi reservoir.
As shown in Table 2 and Figure 8, Maridugao and Tagoloan have the highest mean soil erosion rates among the watershed clusters. From 2010 to 2020, the Tagoloan watershed had the highest maximum soil erosion rates. This can be attributed to a combination of higher LS-factor and R-factor values in the region. In fact, the cluster with the highest LS-factor can be found in Tagoloan watershed cluster. In contrast, the Davao-Saug cluster in 2010, the Agusan-Cugman cluster in 2015, and the Cagayan watershed cluster in 2020 all have lower maximum values of soil loss rates. Table 2 contains statistical values for soil erosion rates predicted by RUSLE by watershed cluster.
Several areas in Bukidnon had extremely high values of the predicted soil erosion rates, as shown in Figure 7. This could be the result of setting the P-factor value equal to one for the entire province and could also be attributed to the resolution of the DEM that was used to generate the LS-factor. The resolution may affect the computation of the LS-factor. Hrabalíková and Janeček [35] predicted soil loss rates closer to those observed while using a similar LS-factor expression to that used in this study, and a DEM with 1 m resolution. Hence, to avoid overestimating the predicted soil loss, a much higher DEM resolution would be preferable, as would documentation of conservation practices in the area.
Based on the predicted values of soil erosion rates shown in Table 3, the areas with very severe soil erosion in Bukidnon decreased slightly in 2020 compared to 2015, by about 1.77%, equivalent to 16,564 ha. Though it can be observed that there was a slight increase in the severe areas from 2015 to 2020, the moderate and high soil erosion areas had decreased by 2.03, and 3.66%, respectively. This decrease has led to an increase in the extent of areas under none to slight soil loss categories by about 7.1%, equivalent to 66,445 ha. When the results in Table 3 were compared with those reported by Adornado and Yoshida [21], the extent of severe and very severe areas had increased by about twofold. The extent of severe areas went from 6.66% to 13.45% on average, while very severe areas went from 15.19% to as high as 34.16% on average. The twofold increase on these areas has been prevalent since 2010. This may be due to the decrease in the areas classified as none to slight, high, and very high soil loss categories. The differences of the results may also be attributed to the nature of the model inputs used in the earlier study. In that study, a DEM was derived from 100 m interval contour lines to obtain an LS-factor. The generated DEM they used had less details about the topography of the province and this could have affected the calculation of the LS-factor. Additionally, their C-factor was generated from an analog land cover map and the ASTER image available during that time [21]. Nevertheless, both results indicate the potential severity of soil erosion in Bukidnon as the rates of the very high to very severe soil erosion areas have not dropped significantly from 2015 to 2020, based on the predicted soil erosion rates, and this poses a serious problem on the sustainability of land and water resources in this province.
3.3. Soil Loss under Annual Rainfall Change Scenarios
As expected, the predicted rate of soil loss (Table 4) decreased under the very severe category while it increased under none to slight, high, very high, and severe soil loss categories, in all the rainfall scenarios, as all generated rainfall scenarios had negative percent changes. The areas categorized as none to slight increased with a range of 0.4 to 1.37% against the baseline scenario, equivalent to 3743 to 12,821 ha. The extent of the very severe areas fell within the range of 0.79 to 2.87%, which means around 7393 to 26,858 ha of land will no longer experience very severe soil loss if annual rainfall decreases by 3.17 to 12.61% in future decades. However, this reduction in very severe areas will also result in an expansion of areas classified as having high, very high, and severe soil loss. In addition, approximately no more than 0.9% (8423 ha) of the province of Bukidnon may experience a decrease in the moderate soil loss rates areas. The extent of very severe soil loss areas may reduce significantly but remains small compared to the total size of the province.
Table 5 and Table 6 depict the extent of areas under the baseline scenario and the third rainfall change scenario (R3) at various soil loss categories for each watershed cluster. The extent of areas for each soil loss category are expressed in percent relative to the total area of each watershed cluster. The results indicate that the extent of the none to slight soil loss category for most of the watershed clusters will increase relative to baseline values under the R3 scenario. By contrast, the extent of the very severe areas will decrease relative to baseline values under the R3 scenario. The extent of the none to slight soil loss categories will increase between 0.48 and 1.98%, whereas the extent of the very severe areas will decrease between 2.49 and 3.89% of the area of the watershed clusters.
As shown in Table 6, the Upper Pulangi, Davao-Salug, and Cagayan watershed clusters will consistently and significantly have larger areas with none to slight soil loss. On the other hand, Lower Pulangi, Maridugao, and Agusan-Cugman watershed clusters will have a smaller proportion of areas with none to slight soil loss but will have a greater proportion of areas with very severe soil loss in future decades. Relative to the watershed size, an extremely large extent of area with very severe soil loss will be in the Lower Pulangi watershed cluster. The distribution of soil loss, particularly in Maridugao and Lower Pulangi watershed clusters, appears to be negatively skewed. This means that there is a greater proportion of areas that have experienced very high to very severe soil loss, as shown in Table 5, and this will still happen, as shown in Table 6, despite an overall decrease in the very severe areas and an increase in the none to slight soil loss areas, as shown in Table 4. Significant variation can be found when comparing the soil erosion rates of the baseline scenario (Table 5) and the R3 scenario (Table 6) in all soil loss categories, except for the none to slight and very severe categories. The extent of areas under the moderate soil loss category in Agusan-Cugman and Lower Pulangi will potentially increase while the rest of the watershed clusters will decrease. All the watershed clusters will experience increases in the high, very high and severe soil loss rate categories.
The predicted soil loss using the RUSLE model at various time periods and under the rainfall change scenarios provided insights on the soil loss behavior in Bukidnon watershed areas. GIS was crucial in identifying the distribution of soil loss both in provincial and watershed cluster scales. The predicted soil loss under various rainfall change scenario may serve as baseline information for determining the potential soil loss in future decades under future rainfall conditions. Results showed that a reduction of approximately 12.61% of the annual rainfall could reduce the extent of very severe areas up to 2.87% and could increase the none to slight soil loss areas at 1.37% of the total area of Bukidnon. This implies that the reduction in rainfall alone will have little effect on the severity of erosion in the province. This suggests that increasing land cover and adapting soil conservation measures may be a more effective way to mitigate the severity of soil erosion than just anticipating for the province to experience less rainfall. Moreover, it is necessary to map and monitor areas adopting and implementing soil conservation measures and practices in order to compare the soil erosion in the areas with and without conservation measures. Not only will this help reduce soil loss, but it will also allow for a more accurate prediction of soil erosion. Nevertheless, the maps presented in this study may help planners in identifying priority areas for soil conservation measures, particularly those with excessive soil loss. Neighboring provinces downstream of the major rivers in Bukidnon should consider coordinating their efforts to implement conservation measures with those in the province of Bukidnon, as they are not exempted from the effects of excessive soil erosion in Bukidnon.
4. Discussion
The RUSLE model has proven to be a practical method for assessing soil erosion risk at the watershed scale and the impact of the rainfall change to soil erosion-prone areas in the province of Bukidnon. Excessive soil loss occurs on steep hillslopes with less vegetation that are devoid of support and conservation practices. As observed during the generation of the LS-factor, the LS-factor is sensitive to the resolution of DEM. The higher the resolution of DEM, the wider the range of values of the LS-factor. This is because the expression of the LS-factor resulted in much higher maximum values at a higher resolution [39], which may add to the uncertainty of the predicted soil erosion rates. Michalopoulou et al. [39] suggested a method to avoid overestimation of the LS-factor; however, it has not been evaluated and validated whether this strategy prevents overestimation or can lead to underestimation of the LS-factor. In addition, different land cover classification of land cover maps at different time periods made it difficult to compare the predicted soil loss rates between each period and from earlier studies, as changes in land cover classification entail a different C-factor value. Nevertheless, the results indicated that predicted soil loss under drier rainfall change scenarios was less significant compared to the size of the province, implying that the severity of soil erosion in the province may not necessarily be reduced just by experiencing less rainfall alone. Increased land cover and the adoption of conservation measures such as conservation agriculture may reduce the risk of extreme soil erosion more than anticipating future rainfall reductions.
The new updated information generated in this study could serve as the basis for the formulation of policies geared towards soil conservation and environmental protection in the province. The approach developed and employed may also be extended to other erosion-prone provinces and regions in the country. The application and integration of RUSLE and GIS to identify the areas most vulnerable to soil erosion will allow policymakers and decision-makers to identify priority areas to focus in implementing soil erosion control measures in the future. It is highly recommended to use DEM with higher resolution, whenever available, and promote the implementation and documentation of soil conservation and support practices in the province. It is also recommended that future studies consider alternative expressions for generating the LS-factor that would limit its overestimation and underestimation, and to use field measurements for evaluation and validation purposes.
Conceptualization, I.G.D. and V.B.E.; methodology, I.G.D.; formal analysis, I.G.D. and V.B.E.; writing—original draft preparation, I.G.D. and V.B.E.; writing—review and editing, I.G.D. and V.B.E.; visualization, I.G.D.; supervision, V.B.E.; funding acquisition, I.G.D. and V.B.E. All authors have read and agreed to the published version of the manuscript.
The study was approved by the University of the Philippines Los Baños, Graduate School.
Not applicable.
Not applicable.
This study was conducted under the auspices of DOST-SEI and ERDT program. Acknowledgement is due to all agencies that provided the needed data used in this study.
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.
Footnotes
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Enlarge this image.Figure 4. The R-factor map of Bukidnon derived from the average annual rainfall for the period (a) 1998–2009, (b) 2003–2014, and (c) 2008–2019.
Enlarge this image.Figure 7. The annual soil loss map of Bukidnon in (a) 2010, (b) 2015, and (c) 2020.
Rainfall change scenarios.
| Scenario | Rainfall Change (%) | Description |
|---|---|---|
| R1 | −3.70 | Dry: Early future (2020–2039) |
| R2 | −8.30 | Very dry: Late future (2080–2099) |
| R3 | −12.61 | Extremely dry: Late future (2080–2099) |
Statistics of the predicted RUSLE erosion rates per watershed cluster.
| Watershed |
Soil Erosion Rates (t/ha/y) | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| Mean | Minimum | Maximum | Standard Deviation | |||||||
| 2010 | 2015 | 2020 | 2010 | 2015 | 2020 | 2010 | 2015 | 2020 | ||
| Agusan-Cugman | 363 | 338 | 300 | 0 | 42,192 | 29,127 | 51,003 | 564 | 479 | 443 |
| Cagayan | 227 | 241 | 220 | 0 | 46,693 | 47,834 | 35,126 | 421 | 403 | 394 |
| Tagoloan | 466 | 458 | 472 | 0 | 84,508 | 91,421 | 114,275 | 789 | 742 | 795 |
| Maridugao | 475 | 482 | 441 | 0 | 42,155 | 37,769 | 41,749 | 645 | 646 | 605 |
| Lower Pulangi | 389 | 449 | 414 | 0 | 58,496 | 61,037 | 49,218 | 556 | 580 | 551 |
| Upper Pulangi | 209 | 281 | 265 | 0 | 46,396 | 54,213 | 47,484 | 409 | 514 | 498 |
| Davao-Salug | 196 | 442 | 401 | 0 | 25,556 | 55,941 | 44,924 | 361 | 697 | 654 |
The extent of erosion in Bukidnon 1.
| Category | Soil Loss Rates |
Area Affected by Erosion (%) | ||
|---|---|---|---|---|
| 2010 | 2015 | 2020 | ||
| None to slight | 0–5 | 18.79 | 18.40 | 25.50 |
| Moderate | 5–15 | 9.55 | 9.68 | 7.65 |
| High | 15–50 | 11.14 | 10.12 | 6.46 |
| Very high | 50–150 | 15.34 | 12.20 | 12.32 |
| Severe | 150–300 | 14.87 | 12.62 | 12.87 |
| Very severe | >300 | 30.31 | 36.97 | 35.20 |
Note: 1 Bukidnon total area = 935,846 ha, based on the calculation using GIS.
The extent of soil erosion under rainfall change scenarios in Bukidnon 1.
| Category | Area Affected by Erosion (%) | |||
|---|---|---|---|---|
| Baseline (B) | R1 | R2 | R3 | |
| None to slight | 25.64 | 26.04 | 26.54 | 27.01 |
| Moderate | 7.55 | 7.29 | 6.96 | 6.65 |
| High | 6.50 | 6.65 | 6.86 | 7.07 |
| Very high | 12.39 | 12.67 | 13.04 | 13.41 |
| Severe | 12.97 | 13.20 | 13.49 | 13.77 |
| Very severe | 34.96 | 34.16 | 33.12 | 32.09 |
Note: 1 Total area considered = 935,846 ha, based on the calculation using GIS.
Percent (%) of the extent of soil erosion under baseline rainfall scenario in Bukidnon 1.
| Watershed |
Area (ha) | Extent of Soil Erosion (%) | |||||
|---|---|---|---|---|---|---|---|
| None to Slight | Moderate | High | Very High | Severe | Very |
||
| Agusan-Cugman | 22,202 | 9.43 | 7.54 | 17.40 | 19.06 | 13.98 | 32.61 |
| Cagayan | 124,202 | 31.38 | 9.13 | 6.80 | 14.56 | 14.06 | 24.08 |
| Tagoloan | 167,959 | 21.95 | 6.10 | 7.36 | 11.94 | 10.83 | 41.86 |
| Maridugao | 65,233 | 14.54 | 3.55 | 5.87 | 14.46 | 16.43 | 45.14 |
| Lower Pulangi | 164,739 | 10.10 | 4.09 | 7.93 | 15.88 | 17.34 | 44.69 |
| Upper Pulangi | 337,527 | 34.39 | 10.65 | 5.34 | 10.36 | 11.29 | 28.03 |
| Davao-Salug | 53,986 | 36.87 | 4.53 | 2.32 | 5.65 | 9.71 | 40.97 |
Note: 1 Total area considered = 935,846 ha, based on the calculation using GIS.
Percent (%) of the extent of soil erosion under R3 rainfall change scenario in Bukidnon.
| Watershed |
Area (ha) | Extent of Soil Erosion (%) | |||||
|---|---|---|---|---|---|---|---|
| None to Slight | Moderate | High | Very High | Severe | Very |
||
| Agusan-Cugman | 22,202 | 10.34 | 7.80 | 18.40 | 19.31 | 14.40 | 29.77 |
| Cagayan | 124,202 | 32.99 | 7.94 | 7.55 | 15.78 | 14.34 | 21.39 |
| Tagoloan | 167,959 | 23.19 | 5.32 | 8.06 | 12.49 | 11.65 | 39.33 |
| Maridugao | 65,233 | 15.04 | 3.42 | 6.59 | 15.98 | 17.40 | 41.55 |
| Lower Pulangi | 164,739 | 10.58 | 4.15 | 8.70 | 17.27 | 18.52 | 40.80 |
| Upper Pulangi | 337,527 | 36.38 | 9.14 | 5.69 | 11.33 | 11.99 | 25.54 |
| Davao-Salug | 53,986 | 38.29 | 3.31 | 2.49 | 6.56 | 11.15 | 38.26 |
References
1. El Jazouli, A.; Barakat, A.; Ghafiri, A.; El Moutaki, S.; Ettaqy, A.; Khellouk, R. Soil erosion modeled with USLE, GIS, and remote sensing: A case study of Ikkour watershed in Middle Atlas (Morocco). Geosci. Lett.; 2017; 4, 25. [DOI: https://dx.doi.org/10.1186/s40562-017-0091-6]
2. Salvacion, A.R. Delineating soil erosion risk in Marinduque, Philippines using RUSLE. GeoJournal; 2022; 87, pp. 423-435. [DOI: https://dx.doi.org/10.1007/s10708-020-10264-7]
3. Hu, Y.; Tian, G.; Mayer, A.L.; He, R. Risk assessment of soil erosion by application of remote sensing and GIS in Yanshan Reservoir catchment, China. Nat. Hazards; 2015; 79, pp. 277-289. [DOI: https://dx.doi.org/10.1007/s11069-015-1841-4]
4. Tania, M. Assessment of soil erosion risk in Fizes River Catchment using USLE model and GIS. ProEnvironment; 2013; 6, pp. 595-599.
5. Ella, V.B. Simulating soil erosion and sediment yield in small upland watersheds using the WEPP model. Land Use Change in Tropical Watersheds: Evidence, Causes and Remedies; Coxhead, I.; Shively, G. CAB International: Wallingford, UK, 2005; pp. 9-125.
6. Deutsch, W.G.; Busby, A.L.; Orprecio, J.L.; Bago-Labis, J.P.; Cequiña, E.Y. Community-based hydrological and water quality assessments in Mindanao, Philippines. Forests, Water and People in the Humid Tropics; Bonell, M.; Bruijnzeel, L.A. Cambridge University Press: Cambridge, UK, 2005; pp. 134-150.
7. Tehrany, M.S.; Shabani, F.; Javier, D.N.; Kumar, L. Soil erosion susceptibility mapping for current and 2100 climate conditions using evidential belief function and frequency ratio. Geomat. Nat. Hazards Risk; 2017; 8, pp. 1695-1714. [DOI: https://dx.doi.org/10.1080/19475705.2017.1384406]
8. Fathizad, H.; Karimi, H.; Alibakhshi, S.M. The estimation of erosion and sediment by using the RUSLE model and RS and GIS techniques (Case study: Arid and semi-arid regions of Doviraj, Ilam province, Iran). Int. J. Agric. Crop Sci.; 2014; 7, pp. 304-314.
9. Ghosal, K.; Das Bhattacharya, S. A review of RUSLE Model. J. Indian Soc. Remote Sens.; 2020; 48, pp. 689-707. [DOI: https://dx.doi.org/10.1007/s12524-019-01097-0]
10. Gashaw, T.; Tulu, T.; Argaw, M. Erosion risk assessment for prioritization of conservation measures in Geleda watershed, Blue Nile basin, Ethiopia. Environ. Sys. Res.; 2018; 6, 1. [DOI: https://dx.doi.org/10.1186/s40068-016-0078-x]
11. Andreoli, R. Modeling erosion risk using the RUSLE equation. QGIS and Applications in Water and Risks; Baghdadi, N.; Mallet, C.; Zribi, M. Wiley-Blackwell: Hoboken, NJ, USA, 2018; pp. 245-282.
12. Bekele, B.; Gemi, Y. Soil erosion risk and sediment yield assessment with universal soil loss equation and GIS: In Dijo watershed, Rift valley Basin of Ethiopia. Model. Earth Syst. Environ.; 2021; 7, pp. 273-281. [DOI: https://dx.doi.org/10.1007/s40808-020-01017-z]
13. Belasri, A.; Lakhouili, A. Estimation of soil erosion risk using the universal soil loss equation (USLE) and geo-information technology in Oued El Makhazine Watershed, Morocco. J. Geogr. Inf. Syst.; 2016; 8, pp. 98-107. [DOI: https://dx.doi.org/10.4236/jgis.2016.81010]
14. Mhangara, P.; Kakembo, V.; Lim, K.J. Soil erosion risk assessment of the Keiskamma catchment, South Africa using GIS and remote sensing. Environ. Earth Sci.; 2012; 65, pp. 2087-2102. [DOI: https://dx.doi.org/10.1007/s12665-011-1190-x]
15. Duarte, L.; Teodoro, A.C.; Gonçalves, J.A.; Soares, D.; Cunha, M. Assessing soil erosion risk using RUSLE through a GIS open source desktop and web application. Environ. Monit. Assess.; 2016; 188, 351. [DOI: https://dx.doi.org/10.1007/s10661-016-5349-5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27184749]
16. Tayebi, M.; Tayebi, M.H.; Sameni, A. Soil erosion risk assessment using GIS and CORINE model: A case study from western Shiraz, Iran. Arch. Agron. Soil Sci.; 2017; 63, pp. 1163-1175. [DOI: https://dx.doi.org/10.1080/03650340.2016.1265106]
17. PAGASA. Climate Change in the Philippines; PAGASA: Quezon City, Philippines, 2011.
18. DOST-PAGASAManila ObservatoryAteneo de Manila University. Philippine Climate Extremes Report 2020: Observed and Projected Climate Extremes in the Philippines to Support Informed Decisions on Climate Change Adaptation and Risk Management; Philippine Atmospheric, Geophysical and Astronomical Services Administration: Quezon City, Philippines, 2021; 145p.
19. Paningbatan, E.P. Identifying soil erosion hotspots in the Manupali river watershed. Land Use Change in Tropical Watersheds: Evidence, Causes and Remedies; Coxhead, I.; Shively, G. CAB International: Wallingford, UK, 2005; pp. 126-132.
20. Alibuyog, N.R.; Ella, V.B.; Reyes, M.R.; Srinivasan, R.; Heatwole, C.; Dillaha, T. Predicting the effects of land use change on runoff and sediment yield in Manupali river subwatersheds using the SWAT model. Int. Agric. Eng. J.; 2009; 18, pp. 15-25.
21. Adornado, H.; Yoshida, M. Assessing the adverse impacts of climate change: A case study in the Philippines. J. Dev. Sus. Agric.; 2010; 5, pp. 141-146.
22. Goddard Earth Sciences Data and Information Services Center. TRMM (TMPA-RT) Near Real-Time Precipitation L3 1 day 0.25 degree × 0.25 degree V7, Edited by Andrey Savtchenko, Greenbelt, MD, Goddard Earth Sciences Data and Information Services Center (GES DISC) 2016. Available online: https://disc.gsfc.nasa.gov/datasets/TRMM_3B42RT_Daily_7/summary (accessed on 20 April 2021).
23. Peralta, J.C.A.C.; Narisma, G.T.T.; Cruz, F.A.T. Validation of high-resolution gridded rainfall datasets for climate applications in the Philippines. J. Hydrometeorol.; 2020; 21, pp. 1571-1587. [DOI: https://dx.doi.org/10.1175/JHM-D-19-0276.1]
24. BSWM. Updating the Harmonized World Soil Database (HWSD): Correlation of Philippine Soils into FAO’s World Reference Base for Soil Resources (WRB); BSWM: Quezon City, Philippines, 2013.
25. FAOIIASAISRICISSCASJRC. Harmonized World Soil Database (Version 1.2); FAO: Rome, Italy, IIASA: Laxenburg, Austria, 2012.
26. Wischmeier, W.H.; Smith, D.D. Predicting Rainfall Erosion Losses: A Guide to Conservation Planning; Agriculture Handbook No. 537; U.S. Department of Agriculture: Washington, DC, USA, 1978.
27. Blanco, A.C.; Nadaoka, K. A comparative assessment and estimation of relative soil erosion rates and patterns in Laguna Lake watershed using three models: Towards development of an erosion index system for integrated watershed-lake management. Proceedings of the Symposium on Infrastructure Development and the Environment; Quezon City, Philippines, 7–8 December 2006.
28. Adornado, H.A.; Yoshida, M.; Apolinares, H.A. Erosion vulnerability assessment in REINA, Quezon Province, Philippines with raster-based tool built within GIS environment. Agric. Inf. Res.; 2009; 18, pp. 24-31. [DOI: https://dx.doi.org/10.3173/air.18.24]
29. Chatterjee, S.; Krishna, A.P.; Sharma, A.P. Geospatial assessment of soil erosion vulnerability at watershed level in some sections of the Upper Subarnarekha River basin, Jharkhand, India. Environ. Earth Sci.; 2014; 71, pp. 357-374. [DOI: https://dx.doi.org/10.1007/s12665-013-2439-3]
30. Kim, S.M.; Choi, Y.; Suh, J.; Oh, S.; Park, H.D.; Yoon, S.H. Estimation of soil erosion and sediment yield from mine tailing dumps using gis: A case study at the Samgwang mine, Korea. Geosystem Eng.; 2012; 15, pp. 2-9. [DOI: https://dx.doi.org/10.1080/12269328.2012.674426]
31. Lee, G.S.; Lee, K.H. Scaling effect for estimating soil loss in the RUSLE model using remotely sensed geospatial data in Korea. Hydrol. Earth Syst. Sci. Discuss; 2006; 3, pp. 135-157.
32. David, W.P. Soil and Water Conservation Planning: Policy Issues and Recommendations. J. Philipp. Dev.; 1988; 15, pp. 47-84.
33. Moore, I.D.; Burch, G.J. Physical basis of the length-slope factor in the universal soil loss equation. Soil Sci. Soc. Am. J.; 1986; 50, pp. 1284-1288. [DOI: https://dx.doi.org/10.2136/sssaj1986.03615995005000050042x]
34. Moore, I.D.; Wilson, J.P. Length-slope factors for the revised universal soil loss equation: Simplified method of estimation. J. Soil Water Conserv.; 1992; 47, pp. 423-428.
35. Hrabalíková, M.; Janeček, M. Comparison of different approaches to LS factor calculations based on a measured soil loss under simulated rainfall. Soil Water Res.; 2017; 12, pp. 69-77. [DOI: https://dx.doi.org/10.17221/222/2015-SWR]
36. Delgado, M.E.M.; Canters, F. Modeling the impacts of agroforestry systems on the spatial patterns of soil erosion risk in three catchments of Claveria, the Philippines. Agroforest. Syst.; 2012; 85, pp. 411-423. [DOI: https://dx.doi.org/10.1007/s10457-011-9442-z]
37. Rellini, I.; Scopesi, C.; Olivari, S.; Firpo, M.; Maerker, M. Assessment of soil erosion risk in a typical Mediterranean environment using a high resolution RUSLE approach (Portofino promontory, NW-Italy). J. Maps; 2019; 15, pp. 356-362. [DOI: https://dx.doi.org/10.1080/17445647.2019.1599452]
38. Gong, W.; Liu, T.; Duan, X.; Sun, Y.; Zhang, Y.; Tong, X.; Qiu, Z. Estimating the soil erosion response to land-use land-cover. Water; 2022; 14, 742. [DOI: https://dx.doi.org/10.3390/w14050742]
39. Michalopoulou, M.; Depountis, N.; Nikolakopoulos, K.; Boumpoulis, V. The significance of digital elevation models in the calculation of LS factor and soil erosion. Land; 2022; 11, 1592. [DOI: https://dx.doi.org/10.3390/land11091592]
40. FAOIFAD. Strategic Environmental Assessment, an Assessment of the Impact of Cassava Production and Processing on the Environment and Biodiversity; FAO: Rome, Italy, 2001.
41. Rola, A.C.; Suminguit, V.J.; Sumbalan, A.T. Realities of the Watershed Management Approach: The Manupali Watershed Experience; PIDS Discussion Paper Series, No. 2004-23; Philippine Institute for Development Studies (PIDS): Makati City, Philippines, 2005.
42. Marin, R.A.; Jamis, C.V. Soil erosion status of the three sub-watersheds in Bukidnon Province, Philippines. AES Bioflux; 2016; 8, pp. 194-204.
© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Details
1 Department of Agricultural and Biosystems Engineering, College of Engineering, Central Mindanao University, Musuan, Maramag 8710, Philippines; Land and Water Resources Engineering Division, Institute of Agricultural and Biosystems Engineering, College of Engineering and Agro-Industrial Technology, University of the Philippines Los Baños, Los Baños 4031, Philippines
2 Land and Water Resources Engineering Division, Institute of Agricultural and Biosystems Engineering, College of Engineering and Agro-Industrial Technology, University of the Philippines Los Baños, Los Baños 4031, Philippines