GIS-Based Soil Erosion Risk Assessment in the Watersheds of Bukidnon, Philippines Using the RUSLE Model

: The sustainability of watersheds for supplying water and for carbon sequestration and other environmental services depends to a large extent on their susceptibility to soil erosion, particularly under changing climate. This study aimed to assess the risk of soil erosion in the watersheds in Bukidnon, Philippines, determine the spatial distribution of soil loss based on recent land cover maps, and predict soil loss under various rainfall scenarios based on recently reported climate change projections. The soil erosion risk assessment and soil loss prediction made use of GIS and the RUSLE model, while the rainfall scenarios were formulated based on PAGASA’s prediction of drier years for Bukidnon in the early-future to late-future. Results showed that a general increase in soil loss was observed in 2015, over the period from 2010 to 2020, although some watershed clusters also showed a declining trend of soil erosion, particularly the Agusan-Cugman and Maridugao watershed clusters. Nearly 60% of Bukidnon has high to very severe soil loss rates. Under extreme rainfall change scenario with 12.61% less annual rainfall, the soil loss changes were only +1.37% and − 2.87% in the category of none-to-slight and very severe, respectively. Results showed that a decrease in rainfall would have little effect on resolving the excessive soil erosion problem in Bukidnon. Results of this study suggest that having more vegetative land cover and employing soil conservation measures may prove to be effective in minimizing the risk of soil erosion in the watersheds. This study provides valuable information to enhance the sustainability of the watersheds. The erosion-prone areas identiﬁed will help decision-makers identify priority areas for soil conservation and environmental protection.


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, 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.

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.

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° (https://disc.gsfc.nasa.gov/datasets/TRMM_3B42_Daily_7/summary (accessed on 10 November 2022)). For this study, the annual accumulated TRMM_3B42_Daily v7 was taken from the Geospatial Interactive Online Visualization ANd aNalysis Infrastructure (GIOVANNI) website (https://giovanni.gsfc.nasa.gov/giovanni (accessed on 10 November 2022)). GIOVANNI is a webbased application developed by the GES DISC that provides access to Earth science remote sensing data (https://gpm.nasa.gov/data/sources/giovanni (accessed on 10 November 2022)). The TRMM 3B42_Daily v7 was used in this study because it has the least overall monthly bias and most closely matches the rainfall distribution observed at weather stations, particularly the dry days and torrential rain days, across the entire Philippines, compared to other gridded rainfall products [23]. During this study, the available TRMM_3B42_Daily v7 data were from 1998 to 2019. The Interferometric Synthetic Aperture Radar (IfSAR) Digital Elevation Mod (DEM) with 5 × 5 m resolution was used to delineate the province's major river wat sheds. This was obtained from the National Mapping and Resource Information Autho ity (NAMRIA) of the Philippines. The DEM was projected into WGS 84/UTM zone 51 The elevation map of Bukidnon is shown in Figure 2a.
The soil map of Bukidnon was extracted from the soil map of the Philippines a was downloaded from the Geoportal Philippines website (https://geoportal.gov.ph (a cessed on 1 December 2022)). For some areas in the soil map with soil texture classified undifferentiated, the Food and Agriculture Organization (FAO) Soil Map of the Phil pines [24] and the Harmonized World Soil Database (HWSD) [25] were used as a referen to determine the soil textural class in those areas. The soil map of Bukidnon is shown Figure 2b.
Three land cover maps for 2010, 2015, and 2020 were used in this study. The la cover maps of Bukidnon were extracted from the 2010, 2015, and 2020 land cover maps the Philippines. The land cover maps and the administrative boundary map were dow loaded via the Geoportal Philippines website (https://geoportal.gov.ph (accessed on 1 D cember 2022)), managed by NAMRIA. The land cover maps are shown in Figure 3.

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 valueadded 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 • (https://disc.gsfc.nasa.gov/datasets/TRMM_3B42_Daily_7/summary (accessed on 10 November 2022)). For this study, the annual accumulated TRMM_3B42_Daily v7 was taken from the Geospatial Interactive Online Visualization ANd aNalysis Infrastructure (GIOVANNI) website (https://giovanni.gsfc.nasa.gov/giovanni (accessed on 10 November 2022)). GIOVANNI is a web-based application developed by the GES DISC that provides access to Earth science remote sensing data (https://gpm.nasa.gov/data/ sources/giovanni (accessed on 10 November 2022)). The TRMM 3B42_Daily v7 was used in this study because it has the least overall monthly bias and most closely matches the rainfall distribution observed at weather stations, particularly the dry days and torrential rain days, across the entire Philippines, compared to other gridded rainfall products [23]. During this study, the available TRMM_3B42_Daily v7 data were from 1998 to 2019.
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 (https://geoportal.gov.ph (accessed on 1 December 2022)). For some areas in the soil map with soil texture classified as undifferentiated, the Food and Agriculture Organization (FAO) Soil Map of the Philippines [24] and the Harmonized World Soil Database (HWSD) [25] were used as a reference to determine the soil textural class in those areas. The soil map of Bukidnon is shown in Figure 2b.
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 (https://geoportal.gov.ph (accessed on 1 December 2022)), managed by NAMRIA. The land cover maps are shown in Figure 3. 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 (https://geoportal.gov.ph (accessed on 1 December 2022)). For some areas in the soil map with soil texture classified as undifferentiated, the Food and Agriculture Organization (FAO) Soil Map of the Philippines [24] and the Harmonized World Soil Database (HWSD) [25] were used as a reference to determine the soil textural class in those areas. The soil map of Bukidnon is shown in Figure 2b.
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 (https://geoportal.gov.ph (accessed on 1 December 2022)), managed by NAMRIA. The land cover maps are shown in Figure 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: 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: 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: where A s 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.

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.

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. Figures 4-6 show the map of the RUSLE factors used in calculating the annual soil erosion rates.

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. Figures 4-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-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 LSfactors 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 (a) (c) (b) 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.  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. 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.  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.

Soil Loss
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.  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. 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.   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.

Soil Loss
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 LSfactor. 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.  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.

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.
Tables 5 and 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.

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.