Spatial estimation of soil erosion risk using RUSLE model in District Swat, Eastern Hindu Kush, Pakistan

Worldwide, the hazard of soil erosion decreases soil efficiency limiting agricultural production (Farhan et al. 2013; Devatha et al. 2015). Soil erosion is leading to siltation in dams and aggradation of rivers and canals posing a high risk of floods (Khan & Hasan 2016). Alongside, the sediments deposition on the watercourses disturbs reservoirs, increases their maintenance costs, and most importantly reduces their life (Dutta 2016). It has a direct impact on the ecosystem, food production, and economy in mountainous regions (Castillo et al. 2007). Many research studies have been undertaken to recognize these conditions; their results control soil loss and environmental restoration (Rickson 2014). Hence, it is necessary to reduce the risk of soil degradation, for agricultural production, ecological well-being, as well as sustainability (Jie et al. 2002).

Soil erosion is a geomorphic process that causes the loss of topsoil and decreases fertility (Ochoa-Cueva et al. 2015). About 10 million hectares of cultivable land are eroded each year (Comino et al. 2017; Tamiru & Wagari 2021) making vast areas unfertile worldwide (Prasannakumar et al. 2012). The explosive population growth, human-induced interventions over the fragile slope, unplanned land use land cover (LULC) changes, immature geology, and humid to sub-humid climatic conditions are the key triggering factors of soil erosion, which has accelerated land degradation (Abdulkadir et al. 2019; Chakrabortty et al. 2020). The soil loss from croplands is faster than its formation posing future food insecurity (Pimentel & Burgess 2013). It is estimated that in the last five decades nearly 2 billion hectares of land in the world have been degraded, which has decreased about 11.9–13.4% crop production (Montgomery 2007). The food requirements are increasing daily due to the growing population while food production is decreasing due to soil erosion (Pimentel & Burgess 2013). In this regard, spatial estimation of soil erosion is used to reduce soil erosion risk and development of crop management systems (Renard et al. 1997).

In Pakistan, more than 76% of the land is exposed to wind and water erosion (Baig et al. 2013). About 11 million hectares of suitable agricultural land out of 23 million hectares is exposed to water erosion. Each year the Indus River System in the upper catchment areas of Hindu Kush, Himalayas, and Karakoram erode around 40 million tons of soil (Khan et al. 2012), resulting in siltation in the Mangla and Terbela dams which has reduced the life and capacity of the mentioned multi-purpose dams. It also lowers electricity production and water availability (Irshad et al. 2007).

Mountainous regions, especially the fragile slopes, are exposed and vulnerable to landslides, soil loss, and the action of climate. The Hindu Kush region is one of the mountainous regions prone to recurrent hazards of soil erosion (Rahman & Shaw 2015). Population growth, anthropogenic activities over the fragile slope, land use changes, immature geology, and humid to sub-humid climatic conditions are the key triggering factors of land degradation (Rahman & Khan 2011). The application of the Revised Universal Soil Loss Equation (RUSLE) model is a recognized technique for soil erosion assessment throughout the world (Wischmeier & Smith 1978; Renard et al. 1997) due to its usefulness, precision, and authenticity (Pandey et al. 2009; Bonilla et al. 2010; Kinnell 2010). The estimation of soil loss is vital for cropland and water management, comprising sediment transportation and storage (Kouli et al. 2009). Therefore, the purpose of this paper is the estimation of soil erosion hazards using the RUSLE model in the Hindu Kush Region of Pakistan. Several models exist for the estimation of soil erosion that significantly vary in their input data. However, RUSLE is one of the commonly used models to estimate soil erosion risk (Pham et al. 2018). The integration of GIS and remote sensing with the RUSLE further enhances its credibility by geo-visualization of risk zones on a cell-by-cell (Raster) basis (Thapa 2020). The soil erosion risk map can assist the decision-making authorities in designing zones specific soil erosion risk reduction strategies and plans. The integration of GIS and RUSLE involves the use of GIS technology to spatially model and analyze the RUSLE, which is a useful tool for estimating soil erosion rates and can be used to simulate erosion scenarios and identify areas of high erosion risk. It can help in decision-making for controlling and preventing soil erosion (Pal 2016; Thakuriah 2023).

Understanding the causes and effects of soil erosion is essential for developing effective solutions to mitigate its impact in the study area. Therefore, this study is carried out in the Eastern Hindu Kush region; it aims to (1) estimate actual and potential soil erosion and (2) identify the high erosion areas using the RUSLE model for District Swat.

District Swat, a mountainous region, is known for its natural beauty and diverse ecosystems. The topography and climate of the region make it a unique location for studying soil erosion. The study area has a complex network of streams and rivers that originate in the mountains and flow through the district, contributing to soil erosion. The district is also home to diverse agricultural practices, which have an impact on soil erosion rates. Soil erosion has significant implications for the long-term sustainability of the region, including impacts on agricultural productivity and environmental degradation. Despite the importance of soil, there has been limited research on soil erosion. This study contributes to filling this gap in knowledge and provides valuable insights into the factors contributing to soil erosion in the region.

The monthly rainfall data was obtained from the Pakistan Metrology Department, Peshawar for 17 years (2003–2016). The soil texture map was obtained from Soil Survey of Pakistan and Soil survey of Pakistan and Digital Soil Map of the World database. The SPOT-5 (Satellite Pour l'Observation de la Terre) satellite image for the year 2016 was acquired from Space and Upper Atmospheric Research Commission, Pakistan with a 5-meter spatial resolution, which was used to extract LULC. The Shuttle Radar Topographic Mission (SRTM) Digital Elevation Model (DEM) with 30-meter spatial resolution was downloaded from the United States Geological Survey (USGS) open source. This DEM was to calculate the elevation as well as RUSLE LS factor. The dataset used in this study for the RUSLE model is shown in Table 1.

It is the resistance of soil to rainfall and runoff. It is a representation of the erodibility of a particular soil, which is the quantity of soil erosion per unit index (tons ha h ha⁻¹ MJ⁻¹ mm⁻¹) from a standard plot 22.1 m long and 9% slope (Goldman et al. 1986). Different soils resist erosion based on their biological, physical, and chemical attributes. Soil with sand or clay shows greater resistance as compared to silty soil. In addition, organic matter also resists erosion to a higher degree (Renard et al. 1997). Wischmeier & Smith (1978) proposed a soil nomograph used to find the K factor value based on soil texture, organic matter percentage, and soil structure. Schwab et al. (1981) summarized these results in a soil erodibility table (Table 2). These values are used to establish the K factor for the District Swat. The soil map obtained from the Soil Survey of Pakistan has shown that the study area includes six major soil texture classes (Figure 4(b)) based on the spatial distribution, namely, silt loam (46.37%), bare rock (27.46%), sandy loam (15.9), fine sandy loam (2.67%), loam (1.03%), and silt clay (0.77%). Glaciers and water occupy the remaining area. The soil map obtained from the Digital Soil Map of the World database (DSMW) shows that the study area has 0.97% organic matter in the soil. These maps were in the form of a shape file so the study area was clipped for cropping data to the area of interest under clip option in Analysis tool and then it was rasterized under the conversion tool in Arc toolbox.

The cover management factor refers to the effect of soil management as well as soil interrupting activities on erosion (Renard et al. 1997). Vegetation cover shields the earth's surface from the direct effect of precipitation. It slows down the velocity of runoff and allows excess water to infiltrate (Mahapatra et al. 2018). C factor values will be high for the barren soil and will be low for the soil having more plant cover (Behera 2015). C values range between 0 and 1. When the land is barren and has no plant coverage, these values will be high, i.e., 1. Uninterrupted forests and thick grasses have low C values.

In the present day, Satellite Pour l'Observation de la Terre (SPOT) 5-satellite image for the year 2016 with a 5-meter spatial resolution obtained from Space and Upper Atmospheric Research Commission, Pakistan was used to prepare the C factor. In the first step, image analysis was carried out in such a way that the study area was extracted from the image. The coordinate system was then changed from WGS (World Geodetic System) 1984 to the Projected Coordinate system UTM (Universal Transverse Mercator) Zone 42. Then, maximum likelihood supervised classification was carried out taking 100 training samples/signature files from each land cover class. The image was then classified into 10 land cover classes, i.e., agriculture, barren land, forest, built-up area, water bodies, glaciers and snow, pastureland, orchards, roads, and rangeland (Figure 4(d)). After image classification, it is necessary to find the accuracy of the classified image to verify the information obtained from this image. In this study, the classified image is verified using visual interpretation technique through freely available high-resolution Google earth images. Then, the classified image was assigned C values from Table 3 extracted from different literature based on the experts' knowledge.

It is the soil erosion ratio with specific support practice to subsequent soil erosion with uphill and downslope farming and tillage practices. These practices affect erosion by altering the slope gradient, flow pattern, and direction (Renard & Foster 1983). P factor values vary from 0 to 1. When terrace farming and contouring practices are employed, the P factor value will be low; on the other hand, it will be high when cultivation is directly made on the hill slopes (Wischmeier & Smith 1978). For the study area, the classified SPOT 5 satellite image and DEM were used to prepare the P factor. The agricultural land was assigned slope values from DEM using the standard table (Table 4) proposed by Wischmeier & Smith (1978), while all the non-agricultural areas were assigned P value = 1. The agricultural and non-agricultural land use were then combined in the union option under Analysis tool and it was rasterized under conversion tool in the Arc toolbox.

This approach is a common alternative when direct measurements of soil erosion are not feasible, and it relies on identifying patterns and associations between soil erosion and other observable phenomena, such as landslides. While this method has some limitations, it can still provide valuable insights into soil erosion patterns and help inform decision-making processes related to erosion control and management.

In this study, the soil erosion rate was estimated by applying the RUSLE model in a GIS environment. It is the first effort to undertake spatial estimation of soil loss and delineation of the study area into the different erosion risk zones. The RUSLE approach estimates the rate of soil loss on the exposed slope for any temporal extent by utilizing five factors, namely, R factor, LS factor, K factor, P factor, and C factor. The combination of steep slopes, intense rainfall, and erodible soil has made the study region more susceptible to erosion. RUSLE is one of the most suitable models for the estimation of soil loss (Gashaw et al. 2018). In many studies, the same model has been implemented having similar environmental and geographical conditions and obtained acceptable results, few are Kumar & Kushwaha (2013), Panagos et al. (2015), Bastola et al. (2019), and Thapa (2020).

Various types of soil have various degrees of susceptibility to erosion. Soil erodibility, or the K factor, is the inherent vulnerability of soil to erosion caused by rainfall and runoff. Numerous physical and chemical characteristics of soil can determine soil erodibility. The RUSLE model considers the physical attributes of soil that play the most crucial role in determining soil erodibility are primary particle size distribution, soil structure, organic matter contents, and soil permeability. In the study, the K factor was prepared from the soil map obtained from the Soil Survey of Pakistan and DSMW. The computed K factor values showed that K values range from 0 to 0.45 tons ha h ha⁻¹MJ⁻¹mm⁻¹ with a mean of 0.26 tons ha h ha⁻¹MJ⁻¹mm⁻¹ (Figure 7(b)). The computed K factor was in vector format, which was rasterized in ArcGIS software and located under the conversion tool. Nearly, 47% of the study area is having silty soil, which is the most sensitive soil to erosion. Higher K values are found in the middle and northern parts of the district. Due to poor infiltration capacity, water erosion is highly expected especially in the northern part of the study region. The relative susceptibility of the soil groups is silt loam > fine sandy loam > loam > sandy loam > silty clay.

There is a strong relationship between the LS factor and soil erosion. Erosion is more in the areas where the slope is steep and slope length is more. The study area is characterized by a higher and steep slope. Nearly 70% of the basin has a steeper slope than 40%. In this study, the LS factor is derived from SRTM DEM, as it is widely used in soil erosion studies because it has gained widespread popularity due to their global coverage and free availability to researchers. LS factor values range between 0 and 4,228 (Figure 7(c)). The lowest LS factor values are found in the southern part of the study area along the River Swat because the slope is very low. The highest LS values are seen in the northern part of the district which is why severe soil erosion areas are mostly found in the northern part, which indicates that erosion is highly influenced by the LS factor in the study area.

Knowledge of LULC has enormous potential for providing valuable insights into the distribution of various land use classes such as agricultural land, forests, barren land, built-up land, and water bodies. This understanding is crucial for developmental planning and erosion studies. In this study, the LULC map shows that the study area consists of different land cover classes such as forest (27.98%), barren land (23.92%), agriculture (21.05%), glacier and snow (16.92%), rangeland (6.31%), and water bodies (1.11%). The accuracy assessment shows that the accuracy of the classified image ranged from 77.9 to 84.8% while the user's accuracy ranges from 82.4 to 86.3% and the producer's accuracy ranged from 73.4 to 81.3%. The C factor values ranged from 0 to 1. The higher C values are found in the northern part of the district due to the presence of bare rocks, while the middle part of the study area has low C values as mostly it is covered with forest. The lower portion of the study has slightly high C values due to the plantation activities, as this area is relatively plain and has fertile soil. The zonal statistic tool of ArcGIS tool result shows a mean and standard deviation of 0.28 and 0.31, respectively. Nearly 42% of the study area has a C value of more than 0.5 (Figure 7(d)). This shows that these areas are more vulnerable to erosion because of the direct impacts of rainfall.

The support practice factor and the cover management factor are related in that they both aim to minimize soil erosion caused by human activities. However, the P factor focuses specifically on the management of runoff by altering its direction, pattern, and speed through management practices like strip cropping, contour tillage, and terraces, while the C factor considers the impact of land use and management practices on soil cover. Like the C factor, the P factor map derived from classified LULC and slope map for the study area has values ranging between 0 and 1 (Figure 7(e)). In the study, the higher values (P = 1) are assigned to the non-agricultural areas, while the agricultural areas, which occupy 21.05% of the study region, are assigned different P values from standard Table 4 based on the slope gradient. As most of the agriculture is done in the southern plain areas of the district that is why this area has lower P values than the other parts of the district. The higher values are found in the northern part of the district because no agriculture is practiced here.

Soil erosion poses a significant threat to the sustainable utilization of land resources in District Swat. In this particular research, a combination of rainfall data, soil data, DEMs, and satellite imagery was employed to create a soil erosion hazard map that highlights areas of potential risk. The findings indicate that the annual soil loss in the study area amounts to 598,504 tons/ha/year, with an average value of 878.94 tons/ha/year. Approximately 53.7% of the region experiences low erosion, while high and very high erosion affects 39.3% of the study area. A comparison between the potential and actual soil erosion maps shows that potential soil erosion exceeds the actual soil erosion due to the absence of cover management and support practices. These findings enable decision-makers to understand the maximum erosion that could occur in a given area and develop scenarios involving land use changes to enhance conservation management.

The study area is exposed and vulnerable to soil erosion hazards because of frequent and high intensity rainfall events, cutting of trees and removal of green cover, texture of soil, steep slopes, land cover conditions, and soil conservation practices. The economic status of the farmers also plays a role in cultivation on slopes in the terraced fields. Sometimes, they are not able to maintain the retaining walls of the fields, which exposes the soil to the action of rainfall and surface runoff. It is the first time effort to undertake spatial estimation of soil loss and delineation of the study area into the different risk zone, no previous study has been conducted to estimate soil loss in the study area. The results of the study need to be validated and as a result the soil erosion map was correlated with the landslide location map established through frequency ratio-based statistical analysis. The relationship between soil erosion and landslides occurrence shows that very high erosion zone has higher probability of landslides, correspondingly, for very low erosion zone, has very low probability of landslide occurrences. Islam et al. (2022) found that steep slopes have higher landslide occurrence zones in the District Swat, similar to our study where high erosion zones are found in these regions.

Similarly, Markose & Jayappa (2016) made projections indicating a higher erosion rate in the Kali River basin located in Karnataka, India. The identification of high-risk soil erosion areas in mountain slopes and steep terrains emphasizes the significant role of topographic factors in influencing soil water erosion. In another study, Thomas et al. (2018) documented severe soil erosion in the steep slopes of the mountainous region within the Muthirapuzha River basin in India. Zerihun et al. (2018) estimated substantial soil loss in the steep slope areas of Northwestern Ethiopia. Kayet et al. (2018) discovered that the steep slope regions of West Singbhum district in Jharkhand, India, experienced high rates of soil erosion. Consequently, the findings of this research align with various studies conducted in different geographical regions with similar conditions.

The findings of this study hold significant importance and value for decision-makers, as they provide valuable insights for identifying the regions most susceptible to soil erosion hazards. Consequently, this research can assist in the development of regional conservation plans, while also offering the potential for expansion to a national level for the formulation of comprehensive soil conservation plans. Soil erosion assessment requires careful consideration of the factors that contribute to soil erosion and the methods used to assess them. It also involves complex models that require significant computational resources and expertise to develop and run. The accuracy of these models can be affected by uncertainties in input data and parameterization. The availability and quality of data on soil erosion processes and factors affecting them can be a significant challenge in soil erosion assessment. Therefore, up-to-date, sufficient data and careful consideration of the factors for the estimation of the RUSLE parameters are essential for the reliable assessment of soil loss. In addition, the effectiveness of different erosion prediction models in accurately predicting soil erosion rates is also required in future studies. This could involve comparing the predictions of different models with actual field measurements of soil erosion rates.

Assessing erosion risk and its spatial distribution in a watershed remains a significant challenge for the scientific community. This study sought to address this issue and investigated the quantitative assessment of soil erosion by water in District Swat Eastern Hindu Kush, Pakistan. The results indicate that the area is significantly vulnerable to soil erosion and if left unchecked, the current rate of soil erosion will likely lead to severe land degradation. Steep slopes, erodible soils, and intense rainfall aggravate the severity of this situation. The analysis concludes that around 34% of the study area is at high risk of soil erosion. The high-risk areas have steep slopes and experience high surface runoff with high stream density. The comparison of the actual and potential soil loss assists in evaluating the impacts of vegetation and support measures in reducing soil erosion, as a substantial increase in vegetation can greatly reduce the amount of soil erosion. The validity of the model was assessed by establishing the relationship between the landslide susceptibility map and soil erosion in the Swat River Basin. As a result, a reasonable correlation between the soil erosion intensity map and data on landslide events has been found. The study concludes that incorporating the RUSLE model into a geographic information system offers numerous benefits. It enables the rational management of qualitative and quantitative data related to various soil degradation factors and facilitates the creation of a comprehensive map displaying the distribution of sensitivity levels to erosion in different areas of the watershed.

This study is part of my PhD research work. The authors would like to thank anonymous reviewers for their useful comments and suggestions.

The paper is not submitted to any other journal.

All authors contributed equally to the research work, analysis, and development of the manuscript.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.