Using Geographic Information Systems and the Analytical Hierarchy Process for Delineating Erosion-Induced Land Degradation in the Middle Citarum Sub-Watershed, Indonesia

Land degradation (LD) is an important issue worldwide because it affects food production and people’s welfare. Many factors cause land degradation, but in humid tropical areas, erosion is the main factor. More than 100 countries including Indonesia are affected by LD. Watershed management should be prioritized owing to budget constraints, while on the other side, the area affected by LD is very large compared to the size of the existing land area. The middle Citarum sub-watershed (MCSW) is one of the most degraded drylands in Indonesia, where the environment can be considered a typical humid tropical watershed. The objective of this study was to map degraded lands and prioritize restoration using a combined approach of the universal soil loss equation (USLE), the analytical hierarchy process (AHP) and geographic information systems (GIS) in a multiple-criteria decision analysis (MCDA) environment. The severity of LD was estimated quantitatively by analyzing the parameters of land use and land cover, slope, soil erosion, productivity, and management. The results indicated that the MCSW is dominated by the potentially degraded land classes (38%), followed by the degraded land classes (21%). The prioritization of LD restoration is suggested in the area of very high and high degraded land. The method developed in this research work could be adopted as a tool to guide decision-makers toward sustainable land resource management in humid tropical watersheds affected by LD.


INTRODUCTION
Land degradation (LD) is a dynamic environmental process that can reduce ecosystem functions and disrupts agricultural production. The implications of LD include environmental setbacks and reduced food security as well as diminished sustainable economy (Zhao et al., 2013). LD refers to the reduced or lost biological and economic productivity of agricultural land or forests resulting from human activities (Expert Group on SDG Indicators [IAEG-SDGs], 2016). In the past two decades, LD has affected more than 20% of the vegetated land area and affects more than 1.5 billion people worldwide (United Nations Convention to Combat Desertification [UNCCD], 2017). The global issue of LD has become the target of the Sustainable Development Goals (SDGs). LD can be caused by physical, chemical, and biological factors of the soil (Brevik et al., 2015), as well as human factors (Khaledian et al., 2017). It has been estimated that around 30% of the world's land is degraded (Lal, 2015;Nkonya et al., 2016). Bai et al. (2008) have estimated that the area of the degraded land worldwide is 18.1 m km 2 , 92% of which due to land mismanagement. In the past 45 years, LD has caused world rice yields to decline by an estimated 1.6-2.7% and financial losses to exceed US$ 10.6 billion (Chen et al., 2012). The depletion of natural resources due to LD can occur in dry and subhumid climates alike (Omuto et al., 2014).
Many natural resource processes can cause LD including erosion, deforestation, soil compaction, soil acidity, salting, and desertification (Turner et al., 2016). Among the various forms of LD, water erosion represents the most important form (Hermassi and Amami, 2018). In Indonesia, the main cause of LD is erosion owing to high rainfall intensity. Erosion-induced LD is an environmental problem that occurs worldwide, causing land productivity to decline (García-Ruiz et al., 2013). Indonesia's land area covers 191.19 m ha, which is dominated by 144.47 m ha (75.60%) of dry land; hence, erosioninduced LD can affect a large area. Indeed, dry land is very vulnerable to erosion and degradation. The area of degraded land in Indonesia in 2018 was estimated at 14 million ha (Indonesian Statistics, 2020), and its restoration will be expensive. Thus, it is necessary to prioritize the handling of land affected by LD. As targeted in the SDGs, degraded land needs to be restored. Funding for watershed restoration is limited, and so it is important to define priority areas in order to optimize the results.
Geographic information systems (GISs) can be used throughout the watershed restoration process. Several studies related to the prioritization of watersheds for conservation actions have been carried out, specifically on soil erosion (Singh et al., 2019;Choudhary et al., 2020), watershed morphometry (Mohammed et al., 2018;Tukura et al., 2021), and the sediment yield index (Jang et al., 2013). Furthermore, Ahmad and Pandey (2018) have used the GIS to map LD. Thus, the GIS can be used to map complex and spatial LD (AbdelRahman, 2014).
Multiple-criteria decision analysis (MCDA) integrates GIS models with decision-making processes and is widely used in complex natural resource and environmental management, including watershed management (Halefom et al., 2019). Many MCDA methods are available, one of which uses the analytical hierarchy process (AHP) (Saaty, 2008) to carry out weighting before being integrated with GIS (Hassan et al., 2015;Widiatmaka et al., 2016). MCDA integrating the AHP with GIS can be applied in the study of LD because the model is very flexible, enabling possible decisions to be made based on both quantitative and qualitative data (Weerakoon, 2014). Many researchers have used MCDA approaches for LD mapping. AbdelRahman et al. (2019) considered soil, slope, rainfall, DEM, land use, and spatial land characteristics. Gessesse et al. (2015) studied LD using the land use/land cover (LULC) dynamics model, while in Iraq, Shareef et al. (2020) combined LULC change with the GIS and AHP model in order to consider LULC parameters such as urban area and vegetation.
Many models have been used to estimate soil loss due to water erosion, one of which is the universal soil loss equation (USLE) (Wischmeier and Smith, 1978). This is one of the most relevant methods for assessing soil loss (Alewell et al., 2019). Several researchers (Kefi et al., 2009;López-García et al., 2020) have integrated the USLE and GIS in their studies.
Citarum watershed is a degraded watershed in Indonesia. It plays an important role in the economy and the lives of the inhabitants of the provinces of West Java, for instance, for water supply for irrigation, industry, and energy supply through electric power generation. The Citarum watershed is one of the most eroded watersheds in Indonesia. Therefore, its conservation needs to be prioritized.
The objective of this research was to estimate the LD caused by soil erosion by applying the USLE along with the GIS and the AHP in MCDA approaches, in order that this degraded land area's restoration can be prioritized. The MCDA considered five parameters: LULC, soil erosion, slope, management, and productivity. One part of the Citarum watershed, namely, the middle Citarum sub-watershed (MCWS), was chosen as the study area, owing it to having the characteristics of natural resources in a humid tropical environment. Referring to the extent of a degraded land area in Indonesia, this method for prioritization needs to be developed in typical areas so that it can be extended to other areas.

Location Description
The MCSW has an extension of 71,091.66 ha or 21.94% of the Citarum watershed. The MCSW is located in Purwakarta Regency, Indonesia. Geographically, the MCSW area extends between 107°13′12.4″-107°35′2.7″ east and 6°26′49.9″-6°46′36.7″ south. The MCSW has a tropical climate with a temperature range of 20-30°C and receives 1,000-4,000 mm of rainfall per year. Some soil types that make up the MCSW are Dystrudepts, Eutrudepts, Hapludalfs, and Hapludands. The MCSW is dominated by a hilly morphology with elevations 200-800 m above sea level and slopes ranging from 3 to 40%.

Data Input
This study used a topographic map at a scale of 1:25,000 (Indonesian Geospatial Information Agency [IGIA], 2015), monthly rainfall data for 2010-2019 (Meteorological, Climatological, and Geophysical Agency, 2020), a soil map at a scale of 1:50,000 (Indonesian Center for Agricultural Land Resources Research and Development, 2017), a map of the official regional spatial plan (ORSP) of Purwakarta Regency (Development Planning Agency, 2012), and SPOT 6 imagery of May 26, 2019. Maps and image interpretation results were processed with GIS to obtain thematic maps in order to make soil erosion maps, LD maps, and a priority watershed restoration map. Mostly, the data type used in this study is in a vector format.

Methods
This research was divided into four stages: calculating the weight of each criterion of LD using the AHP model, calculating soil erosion using the USLE model, preparing thematic maps for mapping LD, and making a map of LD and the prioritization of its handling. The research area was divided into three areas with different LULC: protection forests and conservation forests (PFCFs), production forests, and areas protected outside forest areas (PFPO), and cultivation, agriculture/plantation, and other use areas (CAPO). The differentiation between the three areas was based on the fact that the distinctive LULC and area statuses implied differences in the amount of erosion as well as the conservation measures needed. The ORSP map was used to divide the three areas, covering 5.34, 65.58, and 28.97% for PFCF, PFPO, and CAPO, respectively.

Calculating Weight Using the Analytical Hierarchy Process
The LD map was developed based on the parameters: LULC, slope, soil erosion, productivity, and management. The AHP was applied to weight these parameters, via three steps: creating a hierarchy of problem structures, priority analysis by comparing each parameter in pairs to obtain weights, and verification by calculating the consistency in pairwise comparisons. Five experts in the field of land degradation and watershed management were involved in determining the weights. The consistency of expert opinions was verified by determining the consistency ratio value (Saaty, 2008). The weights were given to describe the extent to which each criterion contributed to LD. Five experts were involved in weight estimation to create a pairwise comparison matrix. A scale of 1-9 was used to identify the importance of the LD parameters: if two criteria were equally important, it was given a weight of 1, while a weight of 9 indicated that one criterion was more important than the others. The sub-criteria scoring assessment was based on the weight of each sub-criterion determined by the experts.
where A is the average annual soil loss per unit area (t ha −1 y −1 ), R is the rainfall erosivity factor (MJ mm ha −1 h −1 y −1 ), K is the soil erodibility factor (t h MJ −1 mm −1 ), LS is the slope length and slope steepness factor, C is the cover and management factor, and P is the support and conservation practice factor.

Rainfall Erosivity Factor
The annual rainfall values of 2010-2019 were used to calculate the R factor based on monthly rainfall using Eq. 2 (Li et al., 2014). The R factor reflects the impact of raindrops on the ground (Wischmeier and Smith, 1978).
where P i is the monthly precipitation.

Soil Erodibility Factor
The K factor reflects the rate of soil loss per rainfall erosivity index. The K factor varies from 0 to 1; a low K value indicates less susceptibility to soil erosion. K values were calculated using Eq. 3 (Hammer, 1978).
where M is (very fine sand % + silt %) x (100-clay %), a is the organic matter content (%), b is the soil structure code, and c is the soil permeability code. Analysis of the soil texture and the soil organic content was carried out in the laboratory of the Department of Soil Science, IPB University, Indonesia.

Slope Length and Slope Steepness
LS is the ratio of soil loss on a standard slope with a length of 72.6 feet and a slope of 9%. The values of length (L) and slope (S) can be gathered from the table in Wischmeier and Smith (1978).

Cover and Management Factor
This factor reflects the effect of planting and tillage on soil erosion rates. It ranges from 1 to 0. A value equal to 1 indicates no land cover and is treated as barren land, whereas a value close to 0 indicates a very strong cover effect, the soil being well-protected from erosion. The value of C in this study was obtained from the C index table (Wischmeier and Smith, 1978) and applied to the land use map.

Conservation Practices Factor
This factor is a reflection of the application of soil conservation practices such as planting in strips and terracing. In this study, the p value was obtained from the land cover and the P table (Wischmeier and Smith, 1978) prevailing in Indonesia.

Preparation Input Map for the Land Degradation Map LULC
The LULC map was obtained from a visual interpretation by onscreen digitizing using ArcGIS 10.6.1 of SPOT 6 satellite imagery. The LULC is classified into 10 classes based on the Indonesian

Productivity
Productivity is an important factor when assessing degraded land in agricultural areas. It was calculated based on comparing the productivity achieved by rice and maize with the average productivity under traditional management. The productivity classifications for LD were divided into very high, high, moderate, low, and very low (Indonesian Ministry of Forestry [IMF], 2013).

Management
Management was used to assess LD based on the completeness of parameters: management aspects, the existence of area boundaries, security and supervision, and whether or not counseling is implemented. The management criteria in determining LD were divided into good, moderate, and bad.

Mapping Land Degradation and Prioritizing Restoration
LD was analyzed using five thematic map inputs. Each criterion was reclassified based on the score of each sub-criterion. An overlay was carried out using ArcGIS software on the reclassified thematic map, multiplied by the respective criterion weights generated from the AHP. Subsequently, the LD map was reclassified into five classes: very highly degraded, highly degraded, degraded, potentially degraded, and non-degraded. The prioritization of watershed management is very important in natural resource management, especially in the context of integrated watershed development and limited watershed restoration funding. Planners and decision makers should therefore pay attention to restoration prioritization. In this study, prioritization restoration was carried out by assigning a rating to the level of LD, divided into three classes: high priority, moderate priority, and low priority.

Weight of Criteria
The pairwise comparison and weight values of the criteria obtained from the AHP are presented in Table 1. According to the judgments of the five experts, the highest weight values were given to the LULC in the PFCF and PFPO areas, with weights of 0.519 and 0.502, respectively. These high weight values can be explained by the fact that in these two forest areas, changes in the LULC easily cause land to become degraded. Such findings are in line with the results of Narendra et al.'s (2019) research on the island of Lombok, Indonesia, which has similar environmental conditions. In CAPO areas, the highest weights were given to management and productivity, which can be explained by the fact that if mismanagement occurs, the result is decreased land quality and thereby degraded land. Management and productivity are thus important factors in agricultural areas, as also insisted by several researchers (Abera et al., 2020;Han et al., 2020). Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 710570 5 AHP is a qualitative approach, in which the primary weakness is subjective and unable to incorporate uncertainty (Hong et al., 2019). The AHP is highly dependent on the selection of appropriate experts to determine the priority scale in pairwise comparisons because it will cause biased opinions among experts (Chen et al., 2021). In this study, anticipation has been done by involving five competent experts in the field of land degradation. This was proved by the consistency of the experts involved in determining the weight, which indicates that this study's results were very good. Indeed, the consistency ratio was below 10%, as displayed in Table 1.

Soil Erosion Predicted
The soil erosion predicted using the USLE model in the MCSW, ranging from 0.87 to 495.30 t ha −1 y −1 , is presented in Figure 1E and Table 2. The results indicated that soil erosion has a wider range than the soil erosion in Indonesia presented by Adimihardja (2008), which ranges from 35 to 220 t ha −1 y −1 . However, the results showed the same trend as recognized by Taslim et al. (2019) in 15 watersheds in Indonesia, where the erosion rate was found to range from 0 to 564 t ha −1 y −1 . The large range found in our study suggests high diversity of natural resources, especially LULC. This range of soil erosion was in the humid tropics, as reported by Labrière et al. (2015).
The results indicated that the MCSW was dominated by the light erosion class (15.20-65 t ha −1 y −1 ), which accounted for 25% and was found in PFPO. The very heavy soil erosion class (>495.30 t ha −1 y −1 ) was found in 7.30% of the study area and was mainly in the PFPO region. PFPO represents the largest area, occupying 65.58% of the study area. The PFPO region is mainly covered by forests, plantations, non-cultivated vegetation, and shrubs. Based on the USLE factor, the PFPO area is characterized by an LS factor of 1.2, a C factor of 0.05-0.6, and a K factor ranging from 0.18 to 0.21. This area is located on an 8-25% slope. The soil textures were found to be soft and medium. These physical conditions, coupled with human factors, cause this area to face all classes of soil erosion, from very severe to slight, albeit dominated by the former. As has been shown by Vahabi and Nikkami (2008), different soil characteristics influence soil erosion.
The CAPO area was found to be dominated by rice fields and settlements, located on a 0-8% slope. The CAPO area is mostly composed of metamorphic rocks with the main soil types Eutrudepts and Dystrudepts. The soil texture ranges from soft to moderately soft, while soil drainage ranges from good to poor, and soil depth from deep to very deep. This area has a C value of 0.05-0.75, a K factor ranging from 0.06 to 0.32, and an LS factor of 0.25. Consequently, erosion in the CAPO area tends to be quite slight. However, in this environment, which is strongly influenced by human activity, erosion can vary (Dotterweich, 2013).
The PFCF area is located on steep slopes (15-40%) with a coarse to medium-coarse soil texture, good to fast soil drainage, and shallow to deep soil depth. This area is dominated by the protected forest and conservation forest. While the USLE factor in this area is the K value ranging from 0.14 to 0.32, the LS factor was found to be 1-4.25, and the C factor 0.03-0.75. The main soil erosion class in this area was found to be slight. However, the very severe erosion class was found in locations bordering the PFPO area due to land clearing activities, which disturb the forest area and cause soil erosion.

Thematic Inputs Map and Degraded Land Map
The impacts of LULC, erosion, slope, productivity, and management on LD are presented in Figure 1 and Table 2. The LULC in the MCSW was found to be dominated by plantations (24.23%) and paddy fields (24.62%). The results of the reclassification based on the LULC indicated that MCSW were dominated by the medium-degraded land classes (45%) ( Figure 1A; Table 2). The impact of slopes on LD is presented in Figure 1B and Table 2. The result of the reclassification of slope Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 710570 7 classes indicated that the MCSW is dominated by potentially degraded classes (36.14%) and slightly degraded class (28.15%) land classes. The impact of productivity on LD is presented in Figure 1C and Table 2. The productivity factor was one of the criteria used to assess LD in the CAPO areas. The reclassification of the productivity factor revealed that 59.68% of the MCSW has low productivity, meaning that it is highly degraded. The impact of management on LD is presented in Figure 1D and Table 2. The results of the management factor reclassification showed that 47.20% (CAPO areas) has moderate management factors, followed by 31.92% (PFPO areas) with poor management, rendering these places prone to LD.

Land Degradation Prioritization
In this article, certain criteria have been specified to study the degree of degradation and then create a map of prioritizing for watershed restoration. The distribution of LD in the MCSW is presented in Figure 2 and Table 3. The results show that the area of degraded land in the MCSW is 29,911.95 ha (48.70%), comprising highly degraded (6.96%), degraded (24.64%), and mediumdegraded classes (17.10%). The area of the non-degraded class is 31,506 ha (51.30%), divided into potentially degraded (36.19%) and non-degraded (15.11%). The degradation class map demonstrates that the area with very high LD is distributed around the forest PFPO areas. This class is mainly characterized by very severe and severe erosion classes. Non-degraded areas are concentrated in the CAPO area, distributed in the southern part of the study area and covering 15.85% of the area. The non-degraded class is dominated by paddy fields and settlements where the soil erosion can be described as either slightly eroded or moderately eroded. The dominant LD class in this study area is "potentially degraded" referring to the land that has not yet become degraded, but that has the potential to become degraded due to the changes in the LULC and inappropriate land use management. Indeed, changes in the LULC can cause severe soil erosion and LD (Firdaus et al., 2014). Inappropriate land use management is the main cause of LD, as has been discussed by Mohawesh et al. (2015). Potentially degraded land spreads into the CAPO and PFPO areas dominated by very slight and slight soil erosion.
The prioritization of watershed management was carried out to identify degraded land with high erosion activity, so that appropriate conservation measures can be taken to minimize soil erosion. The watershed was classified into three priority zones: high priority (highly degraded class and degraded class), medium priority (slightly degraded class), and low priority (potentially degraded class). The prioritization classes can be summarized as follows: 1) high priority, referring to the areas in the MCSW with high soil erosion (28.24%) and hence needing to be prioritized for soil conservation; 2) medium priority (17.94%), characterized by a moderate slope, fine texture, deep solum, good drainage, and metamorphic rocks; and 3) low priority (37.96%), referring to areas facing very light erosion as well as those that are potentially vulnerable, that is, in their existing condition, there is no significant erosion, but if they are not protected, soil erosion will occur. Therefore, soil conservation measures are required to protect against erosion. In our study, prioritization was based on the spatial unit of the LD class. This approach differed from many other studies in watershed prioritization, which have generally been based on prioritizing sub-watershed units (Choudhary et al., 2020). The two approaches are increasingly becoming complementary in watershed management.

CONCLUSION
Most of the study area is characterized by a low rate of soil erosion (38%), while a small part (i.e., 4,275.8 ha, or 7.30%) manifests very severe soil erosion. The resulting map indicates the degree of LD. The results demonstrated that the watershed comprises degraded land (very high, high, and moderate degraded classes), occupying 29,911.95 ha (49%), and non-degraded land (potentially degraded and non-degraded land classes), occupying 31,506 ha (51%) of the study area. The highly degraded class was found to cover an area of 7% of the total area, predominantly in the PFPO areas. The non-degraded class was found to cover an area of 15.85%, particularly in the CAPO areas. The analysis of LD classes and prioritization for handling revealed that 16,538.02 ha (28.24%) of the MCSW areas located on slopes >15% and where soil erosion is very severe or severe need to be treated with soil conservation measures. The results of this study indicate that areas with forest landcover face less LD than cultivated areas, and so forest clearing in cultivation areas needs to be accompanied by land conservation efforts. The study has found that the most significant parameter causing soil erosion-induced LD was LULC as it demonstrated the highest weights in PFCF and PFPO areas. In agricultural areas (CAPO), the highest weights were assigned to the management and productivity parameters. In terms of methodology, the results show the effectiveness of integrating the USLE model with the GIS and AHP in an MCDA environment to map LD classes and the spatial prioritization of watershed restoration. The combination of methodologies used in this study can be used as tools for prioritization management in humid tropical watersheds.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
WA and IN were the main contributors of this article, and both contributed to the conceptualization, methodology, analysis, validation, writing the draft, and article preparation. WD, JS, SL, MD, TT, YS, and DS were the contributor members. WD provided support in conceptualization, methodology and analysis, and writing the manuscript. JS, SL, MD, TT, YS, and DS contributed to the field survey and reviewed the manuscript.