The increasing population, deforestation and conversion of agricultural land to the built-up areas are putting pressure on land resources. Moreover, among land degradation, soil loss is one of the common issues that has had adverse consequences for natural ecosystems, thus affecting livelihoods. The Panjkora River Basin is selected as the study area due to its very fragile soil and having shown regular soil loss activity. In the study area, the scientific communities are consistently insisting on monitoring the LULC changes and exploring the extent of soil loss. To achieve the stated objectives, the RUSLE approach was applied to generate maps of soil loss for the years 1990, 2005 and 2020. The analysis revealed that during the past three decades (1990–2020), the built-up areas have been increased by 20%. Contrary to this, a decrease of 3% in barren land, 2% in area under water, 3% in snow cover and 13% in area under vegetation have been recorded. The analysis further revealed that the maximum actual annual soil loss consistently increased from 5,195 tons/ha/year in 1990 to 6,247 tons/ha/year in 2005 and 8,297 tons/ha/year in 2020. This research implies that geospatial technologies are effective tools for modeling the erosion of soil.

  • The soil provides a diverse array of essential ecosystem services.

  • However, the degradation of land and soils has escalated significantly.

  • In the study area, erosion occurs at a relatively high intensity.

  • This study is unique in its assessment of soil erosion in this area, as no previous assessment of this kind has been conducted.

  • This paper will serve as a foundation for future research in the entire basin.

Since the beginning of human history, numerous civilizations have used the surface of the Earth extensively, resulting in changes. The rate of erosion is predicted to grow in the 21st century due to intensified anthropogenic activities like changes in LULC, which is linked to a higher potent hydrological process characterized by extreme rainfall amplitude (Panagos et al. 2012, 2017). Water-induced erosion of soil is influenced by a variety of natural and human processes, and as a result, it has an indirect or direct effect on natural ecosystems (Leh et al. 2013; Erol et al. 2015). In addition to the natural processes, human alterations have increased erosion rates above normal levels. Boardman et al. (2003) and Tağıl (2007) emphasized that the disturbance of human history (such as deforestation, industrialization, urbanization, and agriculture) has been crucial because of the effects of the surface conditions of the soil. Human activities are therefore the main cause of the degradation of the land (Williams 1991; Morris 2020).

After the population explosion, the second-biggest issue for the world's environment is considered to be erosion-induced soil degradation (Kiassari et al. 2012; Nikkami 2012). The issue has become rather serious as a result of improper and excessive agricultural area usage resulting in soil erosion (Pradhan et al. 2012). This soil erosion results in loss in agriculture productivity and ultimately land destruction (Pimentel 2006; Kidane & Alemu 2015). The ability of soil to retain water is changed by soil erosion. Eroded soil is transported by water flow and dumped in the river system where the flow of water is interrupted (Ganasri & Ramesh 2016; Rawat et al. 2016). The carrying capacity of the river is reduced by runoff from the erosion of top land that is dumped there, which causes siltation in irrigation canals, and damage to turbines and water management systems. Changes in land use and land cover (LULC) and ineffective changes in LULC management strategies have all contributed to an increase in erosion (Sharma et al. 2011).

Numerous research studies have shown that the erosion of soil affects the changes in LULC at various spatial and temporal balances (Sharma et al. 2011; Tadesse et al. 2017). A higher risk of landslides exists in the region with highly erodible soil (Abidin & Abu Hassan 2005; Khosrokhani & Pradhan 2014). Therefore, it is vital to understand LULC and its impact on the dynamics of erosion of soil. To manage LULC more effectively and lower the danger of a major landslide, it would be helpful to understand the pattern and dispersion of change in LULC and the risk of soil erosion.

The soil loss models can be divided into empirical, regression and conceptual models. While models vary their sole purpose is to assess water-induced erosion of soil. The USLE and RUSLE are examples of empirical or regression models (Alkharabsheh et al. 2013). Among the many models available, the RUSLE model is frequently used to access and recognize soil loss due to its low requirement of data easy understanding (Perović et al. 2013). The model is helpful in determining the soil loss distributed spatially over a broad area (Ganasri & Ramesh 2016). It can determine the risk of erosion on a regular grid (Shinde et al. 2010). Integration of geospatial techniques and these models have produced good results when used to assess soil loss (Renard et al. 1991; Uddin et al. 2016). To address the problem of soil erosion, this integration has made evaluation and forecasting simpler and more effective (Tadesse et al. 2017).

This study holds importance as it offers insights into a region characterized by elevated erosion rates. These insights could prove helpful to urban planners and conservationists in their efforts to mitigate and analyze soil erosion. The mitigation and regulation of soil erosion have the potential to enhance agricultural yield and prolong the functional longevity of dams and reservoirs by preventing sediment accumulation. The discoveries of this research will serve as a valuable resource for soil experts, policymakers, irrigation divisions and the soil survey department within the Panjkora Basin, aiding them in implementing effective watershed and sediment control measures.

Environmental setting of the study area

The Panjkora River Basin (PRB) is located in the Hindu Kush Mountain range in the Khyber Pakhtunkhwa province of Pakistan. The Kohistan, Barawal, Jandool, Dir, Maidan and Usherai rivers are its major tributaries. It flows from north to south and joins the Swat River near the village of Qalangi. It is approximately 113 km long and has a total catchment area of 5,905 km2. The PRB lies from 34°38′59″ to 35°48′00″ North latitude and 71°12′58″ to 72°21′54″ East longitude (Figure 1). The range of elevation is from 568 m in the south to about 5,773 m in the north. Southward, the elevation declines to below 1,950 m above mean sea level. The maximum point height at the eastern watershed boundary is 4,173 m above mean sea level. The area experiences hot summers, with temperatures ranging from 15.9 to 33.1 °C and winter in this region is cold, with temperatures dropping below freezing at night. Snowfall is common in the higher elevations of the region, and the mountain passes may be closed due to snow and ice (GOP 2000). In addition, the River Panjkora receives glacier milk (meltwater) from 82 small and large glaciers (Ives et al. 2010). Rainfall occurs throughout the year in the study area. Winter rainfall is caused by the Western Depression, whereas summer rainfall is caused by the Monsoon. Rainfall varies annually between 823 and 2,149 mm (Mahmood et al. 2016; Rahman & Khan 2011).
Figure 1

Location of the study area of the Panjkora River Basin, Eastern Hindu Kush.

Figure 1

Location of the study area of the Panjkora River Basin, Eastern Hindu Kush.

Close modal

Acquisition of data

The district boundary was acquired from a survey of Pakistan and the same administrative boundary was used for the extraction of the PRB from the satellite images. Three temporal Landsat images were acquired from the United States Geological Survey (USGS) Earth Explorer website. Images of the years 1990, 2005, and 2020 were chosen. The rainfall data was obtained from the Pakistan Metrological Department (PMD). Soil data was collected from the soil survey and agriculture department of Khyber Pakhtunkhwa. For slope map and topographical analysis, SRTM data of 30 m were acquired from USGS (Figure 2; Table 1).
Table 1

Datasets and its sources

YearMonthDatasetsSpatial resolutionSources
1990 12 June LANDSAT 5 TM 30*30 USGS 
2005 21 June LANDSAT 7 ETM 30*30 USGS 
2020 29 June LANDSAT 8 OLI/TIRS 30*30 USGS 
– – SRTM DEM 30*30 USGS 
– – Rainfall – Pakistan Metrological Department 
– – Soil Type – Soil Survey Department 
YearMonthDatasetsSpatial resolutionSources
1990 12 June LANDSAT 5 TM 30*30 USGS 
2005 21 June LANDSAT 7 ETM 30*30 USGS 
2020 29 June LANDSAT 8 OLI/TIRS 30*30 USGS 
– – SRTM DEM 30*30 USGS 
– – Rainfall – Pakistan Metrological Department 
– – Soil Type – Soil Survey Department 
Figure 2

Flowchart of research methodology.

Figure 2

Flowchart of research methodology.

Close modal

Image analysis

Image classification is a vital tool for processing and analysis of satellite imagery. Acquired Landsat images were classified.

Five classes, including the vegetation area, barren land, built-up area, snow cover and water body training samples, were taken. More than 200 training samples of each class were taken. The classified images were assessed for reliability by performing an accuracy assessment in which the classified images were compared with the other source of data that was recognized as reliable, validation data and data retrieved by interpreting satellite imagery.

The estimated outcome offers overall map accuracy as well as the map's accuracy for each class. It can be calculated using the formula:
formula

The assessment included calculating:

  • (i)

    Commission error: The number of incorrectly classified pixels in a row.

  • (ii)

    Omission error: The incorrectly classified pixels in a column.

  • (iii)
    Producer's accuracy: The number of pixels correctly classified in a particular class as a percentage of the total number of pixels actually belonging to that class in the image. Producer's accuracy measures how well a certain area has been classified. It was calculated using the following formula:
    formula
  • (iv)
    User's accuracy: The sum of correctly classified pixels is divided by the total number of pixels. It is calculated by using the following formula:
    formula

Data acquisition and analysis of soil erosion

Application of the RUSLE model

The RUSLE approach provides a good methodology for calculating the erosion of soil, and its contributing causes, which maintains its factors and formulae (Renard et al. 1994). These factors are rainfall erosivity, erodibility of soil, length of slope and steepness, management cover and practice support. It has been computerized to help with the process of soil erosion (Renard & Ferreira 1993). It is mathematically represented as
formula
(1)
where A is the rate of average annual soil erosion (tons/ha/year), R is rainfall erosivity (MJ mm/ha/h/year), K is soil erosivity (tons ha/J/mm), LS is the sloping length and steepness, C is the cover management and P is the support practice factor. In this equation, dimensionless factors are P, C and LS.

Rainfall erosivity (R)

The important initiating water erosion factor is the intensity and quantity of rainfall (Foster et al. 2002). According to EI30 (the sum of energy of kinetic and the rainfall of 30 min intensity), erosion would increase as rainfall intensity and amount increase (Renard 1997). For the R factor calculation, continuous data of long-term rainfall are required; however, in most countries, these EI30 data are not accessible. It is a challenging and time-consuming task even if there are still enough data available. To calculate the factor of R from the data of the annual average rainfall, certain straightforward approaches have been applied in different regions of the world. The availability and dependability of monthly rainfall data is the key benefit of this streamlined method (Lee & Heo 2011). Numerous experts have recognized a strong association between monthly data from various regions of the world and rainfall erosivity (Teh 2011).

For this study, the erosivity factor of R was computed from the annual average rainfall data. Data on rainfall were obtained from the Pakistan Meteorological Department, Peshawar for the two met stations close to Panjkora Basin (Timergara and Dir metrological stations) to estimate the erosivity factor R. This data was used to create a point map, which was then interpolated in ArcMap.

The metrological station data was computed for rainfall of annual average and erosivity of rainfall data. Equation (2) was used which was developed by Singh et al. (1981):
formula
(2)
where R is the erosivity of rainfall and P is the annual average rainfall. Singh et al. (1981) were used because no particular equation exists in Pakistan for estimating the R factor. Due to Pakistan and India's rainfall patterns being comparable, the equation was applied to the research area.

Erodibility of soil (K)

Erodibility is the erosion of soil resistance caused by runoff and rainfall. The vulnerability of different types of soil to erosion of water exhibits varying degrees (Thomas et al. 2018). Numerous soil physical and chemical characteristics have an impact on it. However, the RUSLE model only takes into account the physical characteristics of the soil, such as its structure, particle size, organic content, and permeability, which are the key factors affecting soil erodibility. The soil type map of the PRB was obtained from the soil survey in Pakistan, and soil types were categorized into four classes, namely, Glaciers, Haplic Xerosols, Lithosols and Eutric Cambisols. This was then used to calculate the factor of erodibility by giving the values of K obtained from various sources (Table 2).

Table 2

Panjkora River Basin, soil type and K factor values of the study area

Soil typesK values
Glaciers 
Haplic Xerosols 0.19 
Lithosols 0.25 
Eutric Cambisols 0.34 
Soil typesK values
Glaciers 
Haplic Xerosols 0.19 
Lithosols 0.25 
Eutric Cambisols 0.34 

Length of slope and steepness (LS)

The LS is a mixture of two topographical variables: slope length (L) and slope steepness (S), which greatly influences the erosion of the soil process (Gashaw et al. 2018). On a steep slope, the water rushes faster, causing increased pressure on the surface and, as a result, increasing the transport of many sediments. Slope length, which refers to the distance from the point where streamflow originates to either the point where the slope drops and sedimentation takes place or the point where the water can flow into discrete channels, also contributes to erosion (Haile & Fetene 2012).

Nowadays, DEM with 30 m resolution was used to compute the LS factor in this study, which was obtained from the Earth Explorer website, using Bizuwerk et al. (2003) established equation in the GIS environment for the LS factor calculation (Equation (3)):
formula
(3)

The slope percentage suggests that the value of exponent m ranges from 0.2 to 0.5 (Table 3). Because the majority of the area in the study region has a steeper slope of 5%, 0.5 was taken as the m value from Table 3 for Equation (3).

Table 3

Panjkora River Basin, values of m for different slopes

m valueSlope (%)
0.2 <1 
0.3 1–3 
0.4 3–5 
0.5 >5 
m valueSlope (%)
0.2 <1 
0.3 1–3 
0.4 3–5 
0.5 >5 

Management cover (C)

Factor C takes into account the effects of the combination of cover and activity of management on the loss of soil. Vegetation can significantly slow down the runoff thus protecting soil from erosion. Anthropogenic activities are the major factors of change (Karaburun et al. 2010). Plant cover significantly decreases soil erosion because it intercepts rainwater, slows down rainfall, runoff and speeds up infiltration (Nedd et al. 2021). Remote sensing provides up to date and accurate data of earth surface on a large scale, which is extremely useful for understanding the earth surface dynamics (Karaburun et al. 2010).

Normalized Difference Vegetation Index (NDVI) is found useful in providing knowledge about the vegetation cover (Fujaco et al. 2016). Different algorithms have been used by a large number of researchers using NDVI to determine the C factor (Karaburun et al. 2010). In this study, Durigon et al. (2014) suggested that Equation (4) was applied to calculate the factor of C:
formula
(4)
whereas
formula
(5)

Practice support (P)

The factor of P differs from the factor of C in that it describes the impact of strategic planning on the runoff by changing its direction, pattern and speed (Panagos et al. 2015). Some scientists suggested that the factor of P values is largely influenced by the gradient's slope (Nagendra et al. 2004), while others suggested calculating the value of factor P using agricultural production. The factor P can be calculated in various ways, such as by directly examining the type of land usage in the field and by identifying particular farming practices that are especially time- and money-consuming. Landsat-classified images were used in this study to develop land cover classes and compute the practice support factor (Table 4).

Table 4

Panjkora River Basin, factor P values

Land use land coverSlope (%)P values
Agriculture 0–5 0.10 
5–10 0.12 
10–20 0.14 
20–30 0.19 
30–50 0.25 
50–100 0.33 
Other lands All 
Land use land coverSlope (%)P values
Agriculture 0–5 0.10 
5–10 0.12 
10–20 0.14 
20–30 0.19 
30–50 0.25 
50–100 0.33 
Other lands All 

Spatio-temporal analysis of LULC

The maximum likelihood of the supervised classifier was used to generate LULC maps. The Basin of the Panjkora River area is approximately 5,905 km2. About 600 training samples were taken from the image and provided to the classifier. The imagery of the study area was divided into five different LULC classes: water bodies, vegetation, built-up area, snow cover and barren area. Thus, the images were classified into five classes as shown in Figure 3.
Figure 3

Panjkora River Basin classified images for the years 1990, 2005 and 2020.

Figure 3

Panjkora River Basin classified images for the years 1990, 2005 and 2020.

Close modal

The results show that, in the last three decades, LULC in the study area changed significantly from 1990 to 2020 (Table 5). The built-up areas increased up to 19% while a decrease of 3% in barren land and 13% in vegetation was recorded in the study area.

Table 5

LULC changes detection (area)

Land coverArea in sq. km
LULC 1990LULC 2005LULC 2020
Barren land 1,012.85 910.36 854.58 
Built-up area 374.45 980.44 1,528.6 
Snow cover 574.93 473.69 365.62 
Vegetation 3,760.74 3,404.01 2,949.26 
Water 182.21 136.7 107.03 
Total 5,905 
Land coverArea in sq. km
LULC 1990LULC 2005LULC 2020
Barren land 1,012.85 910.36 854.58 
Built-up area 374.45 980.44 1,528.6 
Snow cover 574.93 473.69 365.62 
Vegetation 3,760.74 3,404.01 2,949.26 
Water 182.21 136.7 107.03 
Total 5,905 

From the analysis of Landsat images, sharp changes occurred in the LULC of the PRB from 1990 to 2020. For the purpose of analysis, the whole PRB was classified into five major LULC classes, namely, barren land, built-up area, snow cover, water bodies and vegetation. The analysis reveals that snow and water in the study area have somewhat varied over the study period. In the same way, vegetation and barren areas have also changed over the study period. It is significantly reduced. The built-up area has increased, which means that rapid urban expansion is taking place in the study area as shown in Figure 4.
Figure 4

LULC for the years 1990, 2005 and 2020.

Figure 4

LULC for the years 1990, 2005 and 2020.

Close modal

Accuracy assessment

The accuracy assessment and Kappa coefficient calculation for the years 1990, 2005 and 2020. The 1990 image showed producer's and user's accuracy scores above 80% in all categories. Overall accuracy and Kappa coefficient index were calculated at 81.16% and 0.81, respectively. In the 2005 image, producer's and user's accuracy scores exceeded 85% in all categories. Overall accuracy and Kappa coefficient index were calculated at 85.91% and 0.85, respectively. The 2020 classified image showed producer's and user's accuracy scores above 90% in all categories with a high accuracy index, with overall accuracy recorded at 90.12% and a Kappa coefficient index of 0.91, respectively, as shown in Table 6.

Extent and evaluation of soil erosion

In the PRB, all five factors of the RUSLE model were used to evaluate the potential annual soil loss and average annual soil loss (Table 6). These factors and model results were in the raster format with 30 m. The model is run using the GIS capacity to analyze spatial data and perform overlay analysis. In this study, two approaches have been applied in estimating soil erosion: potential annual erosion of soil and calculation of actual annual erosion of soil. It can be derived from the following expressions:
formula
(6)
formula
(7)
Table 6

Accuracy assessment

Landsat images (Year)Overall accuracy (%)Kappa coefficient
1990 81.16 0.81 
2005 85.91 0.85 
2020 90.12 0.91 
Landsat images (Year)Overall accuracy (%)Kappa coefficient
1990 81.16 0.81 
2005 85.91 0.85 
2020 90.12 0.91 

Calculation of the rainfall erosivity (R) factor

In the PRB, there are two meteorological stations where rainfall data is regularly acquired. In this regard, for the factor of R computing, the past 30 years' rainfall data of Dir and Timergara were collected from PMD. The factor of R is extremely sensitive in calculating the erosion of soil risk. Using the inverse distance weighted (IDW) method of interpolation, the rainfall map was created, which was then used to calculate the rainfall and erosivity factor (R). As a result, an R factor raster map for the PRB was created using the annual rainfall data of 1990, 2005 and 2020. It can be analyzed from the R factor of 1990 to 2020 that the southeastern part of the Panjkora receives maximum rainfall while the northeastern part of the study area receives less rainfall. The rainfall erosivity values ranged between 289 and 369 MJ mm/ha/h/year in 1990, which increased to 292 and 371 MJ mm/ha/h/year in 2005 and further increased to 294 and 384 MJ mm/ha/h/year in 2020 (Figure 5).
Figure 5

Panjkora River Basin, R factors for the years 1990, 2005 and 2020.

Figure 5

Panjkora River Basin, R factors for the years 1990, 2005 and 2020.

Close modal

Calculation of the soil erodibility (K) factor

The soil type map of the PRB was obtained from the soil survey in Pakistan, and soil types were categorized into four classes, namely, Glaciers, Haplic Xerosols, Lithosols and Eutric Cambisols. This was then used to calculate the factor of erodibility by giving the values of K obtained from various sources of literature. Resistance of soil to erosion is called erodibility, which offers resistance. The soil erodibility K factor values range from 0 to 0.25 tons ha h/ha/MJ/mm. In soil type, the erodibility of Glacier values are 0, Haplic Xerosols are 0.19, Lithosols are 0.25 and Eutric Cambisols are 0.34 and which is a high erodibility value. Lithosols type of soil is found in the Northern, while the Eutric Cambisols soil type is found in the Southern area of the Panjkora (Figure 6).
Figure 6

Panjkora River Basin, soil type and K factor.

Figure 6

Panjkora River Basin, soil type and K factor.

Close modal

Calculation of the slope length and steepness (LS) factor

The LS is a mixture of two topographical variables: slope length (L) and slope steepness (S), which greatly influences the erosion of the soil process. When the length of slope increases gradient, it will accelerate flow velocity and hence maximum the erosion of soil increases the factor of LS. The length of slope and steepness are proportional to the rate of erosion of soil. The length of slope rises as a result of water being deposited into the slope; it accelerates. It was computed from the Earth Explorer website 30 m DEM. This map depicts that elevation ranges between 568 and 5,773 m above mean sea level. The northern part in the mountainous area has high altitude while the southern part has low altitude in the study area, while the slope ranges between 0 and 429%. The erosion risk increases as the slope becomes longer and steeper. Because of the greater accumulation of flow, water erosion of the soil rises with a steep slope. Water rushes faster down a steep slope, causing increased pressure on the surface and, as a result, increasing the transport of many sediments. The LS factor is higher in the northern area. I think it is the effect of a steep slope in that area (Figure 7).
Figure 7

Panjkora River Basin, DEM, slope and LS factor for the years 1990, 2005 and 2020.

Figure 7

Panjkora River Basin, DEM, slope and LS factor for the years 1990, 2005 and 2020.

Close modal

Calculation of the cover management (C) factor

The relative accuracy and low cost of vegetation indicators like NDVI and LULC categorized maps are now preferred over traditional methods. The C factor was calculated in this study using Landsat 5, 7, and 8 satellite imagery downloaded for the years 1990, 2005 and 2020. Higher NDVI values indicate greater vegetation cover, while lower values indicate sparse or no vegetation. The C factor values are opposite to the NDVI. With low C factor values and high NDVI values, the lands have more resistance to soil erosion. The C factor value in the study area ranges from 0.12 to 0.76 in 1990, 0.14 to 0.68 in 2005 and 0.18 to 0.62 in 2020 Due to the sparse vegetation in the northern part of the study area, the C values are high, while they start to decline as you move southward (Figure 8).
Figure 8

Panjkora River Basin, C factors for the years 1990, 2005 and 2020.

Figure 8

Panjkora River Basin, C factors for the years 1990, 2005 and 2020.

Close modal

Calculation of the support practice (P) factor

The P factor is used for water and erosion speed reduction. It is calculated from the map of land use of the research area obtained from Landsat image classification from 1990, 2005 and 2020 based on agricultural systems on various slopes. The Landsat satellite imagery was classified into five classes of different LULC, e.g., built-up area, water bodies, vegetation, barren land and snow cover using supervised classification techniques. The slope layer was overlaid on the land cover layer and allocated various slope values to the vegetation based on the proposed values from the slope map (Table 4), while all non-vegetation land use classes were given a P-value of 1. Where the value of 1 denotes a non-anthropogenic resistance erosion facility and the value 0 indicates that anthropic erosion resistance is very good (Figure 9). Due to the steep slopes in the study area, farming practices have evolved. Generally, both the P and C factors are related because they are both used to decrease the erosion of soil. However, the factor of P differs from the factor of C in that it describes the impact of strategic planning on the runoff by changing its direction, pattern and speed.
Figure 9

Panjkora River Basin, P factors for the years 1990, 2005 and 2020.

Figure 9

Panjkora River Basin, P factors for the years 1990, 2005 and 2020.

Close modal

Potential annual soil loss

This is the ability to erode an area per year. Erosion that might occur without any type of protection is referred to as potential erosion. It provides details about the worst situation that could happen whereas actual erosion of soil represents the current condition, taking into account support practices and management that reduce erosion risk significantly. In this study, the rate of maximum potential loss of soil was 5,853 tons/ha in 1990, which consistently deteriorated in the year 2005 and reached 7,048 tons/ha/year. After analysis of 2020, it was found that the rate of soil loss further multiplied and reached an alarming figure of 9,162 tons/ha/year. The amount of potential erosion of soil is high in the South Eastern area because of the effect of more rainfall, soil type and less vegetation (Figure 10).
Figure 10

Panjkora River Basin potential soil loss from 1990, 2005 and 2020.

Figure 10

Panjkora River Basin potential soil loss from 1990, 2005 and 2020.

Close modal
The calculated erosion of soil rate was divided into five severity categories, namely, free from soil loss (<1 tons/ha/year), low (1–25 tons/ha/year), moderate (25.1–50 tons/ha/year), high (50.1–75 tons/ha/year) and very high (>75 tons/ha/year), as shown in Table 7. The result shows that, in the last three decades of potential soil loss, the area under the free from soil loss category is 50.65%, low at 30.46%, moderate at 9.76%, high at 7.57% and very high at 1.56% in 1990, which increased to 55.16, 30.02, 10.06, 3.12 and 1.65% in 2005 and further increased to 45.69, 30.44, 11.93, 9.36 and 2.60% in 2020 (Figure 11).
Table 7

Panjkora River Basin potential annual soil loss (area)

Erosion categoriesPotential loss (tons/ha/year)Area in sq. km
Area %age
Potential soil loss 1990Potential soil loss 2005Potential soil loss 2020Potential soil loss 1990Potential soil loss 2005Potential soil loss 2020
Free from soil loss <1 2,990.87 2,967.42 2,697.75 50.65 50.25 45.69 
Low 1–25 1,798.77 1,772.49 1,797.32 30.46 30.02 30.44 
Moderate 25.1–50 576.21 593.99 704.27 9.76 10.06 11.93 
High 51–75 447.09 473.97 552.70 7.57 8.03 9.36 
Very high >75 92.26 97.44 153.38 1.56 1.65 2.60 
Total 5,905 100 
Erosion categoriesPotential loss (tons/ha/year)Area in sq. km
Area %age
Potential soil loss 1990Potential soil loss 2005Potential soil loss 2020Potential soil loss 1990Potential soil loss 2005Potential soil loss 2020
Free from soil loss <1 2,990.87 2,967.42 2,697.75 50.65 50.25 45.69 
Low 1–25 1,798.77 1,772.49 1,797.32 30.46 30.02 30.44 
Moderate 25.1–50 576.21 593.99 704.27 9.76 10.06 11.93 
High 51–75 447.09 473.97 552.70 7.57 8.03 9.36 
Very high >75 92.26 97.44 153.38 1.56 1.65 2.60 
Total 5,905 100 
Figure 11

Panjkora River Basin potential annual soil loss for the years 1990, 2005 and 2020.

Figure 11

Panjkora River Basin potential annual soil loss for the years 1990, 2005 and 2020.

Close modal

Actual annual soil loss

There are several hotspots of erosion in the Panjkora Basin. The hotspot areas are maximum in the headwater regions, whereas they are comparatively low in downstream areas. These sub-basins are further divided into a very high-risk zone, high-risk zone, moderate-risk zone, low-risk zone and free from soil loss. The maximum annual soil loss was 5,195 tons/ha/year in 1990 which increased to 6,247 tons/ha/year in 2005 and further increased to 8,297 tons/ha/year in 2020. In the South East area, the amount of actual annual erosion of soil is high because of the effect of more rainfall, soil type and less vegetation (Figure 12).
Figure 12

Panjkora River Basin, annual soil loss from 1990, 2005 and 2020.

Figure 12

Panjkora River Basin, annual soil loss from 1990, 2005 and 2020.

Close modal
The calculated erosion of soil rate was divided into five severity categories, namely, free from soil loss (<1 tons/ha/year), low (1–25 tons/ha/year), moderate (25.1–50 tons/ha/year), high (50.1–75 tons/ha/year) and very high (>75 tons/ha/year), as shown in Table 8. The result shows that, in the last three decades of annual soil loss, the area under the free from soil loss category is 58.89%, low at 25.33%, moderate at 10.97%, high at 3.96% and very high at 0.82 in 1990, which increased to 53.01, 30.70, 12.46, 2.68 and 1.15% in 2005 and further increased to 47.46, 32.10, 11.81, 6.53 and 2.09% in 2020 (Figure 13).
Table 8

Actual annual soil loss (area)

Erosion categoriesAnnual loss (tons/ha/year)Area in sq. km
Area %age
Annual soil loss 1990Annual soil loss 2005Annual soil loss 2020Annual soil loss 1990Annual soil loss 2005Annual soil loss 2020
Free from soil loss <1 3,477.86 3,080.77 2,802.95 58.89 52.16 47.46 
Low 1–25 1,495.99 1,813.13 1,895.55 25.33 30.70 32.10 
Moderate 25.1–50 648.10 665.54 697.38 10.97 11.27 11.81 
High 50.1–75 234.07 278.03 385.65 3.96 4.71 6.53 
Very high >75 48.52 67.83 123.16 0.82 1.15 2.09 
Total 5,905 100 
Erosion categoriesAnnual loss (tons/ha/year)Area in sq. km
Area %age
Annual soil loss 1990Annual soil loss 2005Annual soil loss 2020Annual soil loss 1990Annual soil loss 2005Annual soil loss 2020
Free from soil loss <1 3,477.86 3,080.77 2,802.95 58.89 52.16 47.46 
Low 1–25 1,495.99 1,813.13 1,895.55 25.33 30.70 32.10 
Moderate 25.1–50 648.10 665.54 697.38 10.97 11.27 11.81 
High 50.1–75 234.07 278.03 385.65 3.96 4.71 6.53 
Very high >75 48.52 67.83 123.16 0.82 1.15 2.09 
Total 5,905 100 
Figure 13

Panjkora River Basin actual annual soil loss for the years 1990, 2005 and 2020.

Figure 13

Panjkora River Basin actual annual soil loss for the years 1990, 2005 and 2020.

Close modal

Comparison of soil loss and LULC

The hotspot areas were identified and presented in the form of maps from 1990 to 2020, using the soil erosion inventory. Year-wise hotspots were identified and it was found that, with the LULC change, the rate and erosion of soil area have also increased. The maximum LULC change class occurred from 1990 to 2020 in the snow cover, built-up area, vegetation and barren land. The area of snow, vegetation and barren land has been drastically reduced very rapidly in the study area while the built-up area is increased. In the Panjkora Basin, the hotspot areas of very high erosion of soil increased in 2020 as a comparison to 1990 and 2005 (Figure 14).
Figure 14

Panjkora River Basin, LULC and high-risk zone of soil erosion of 1990, 2005 and 2020.

Figure 14

Panjkora River Basin, LULC and high-risk zone of soil erosion of 1990, 2005 and 2020.

Close modal

Discussion

Extreme erosion of soil decreases not only the fertility and land productivity but also the capacity and effectiveness of dams and reservoirs by adding a significant amount of sediment to them. It takes a lot of time and effort to calculate and map the spatial extent of erosion of soil hazards, but the approach of RUSLE and GIS combination is a highly beneficial tool for measuring and mapping the soil erosion of a region. This research is significant because it provides personal knowledge of an area with high rates of erosion, which could assist planning experts and conservationists in reducing and examining the erosion of soil. The prevention and control of the erosion of soil will contribute to a rise in cultivated land output, and the removal of sediments will lengthen the lifespan of dams and reservoirs.

A study conducted by Maqsoom et al. (2020) employed the RUSLE model to assess the annual erosion of soil within the Chitral River Basin. The findings indicated a computed loss of soil of 58 tons/ha/year from the Chitral River Basin. Notably, the study site corresponds to the northern border of the PRB. Another study by Batool et al. (2021) in Pakistan's Potohar region suggested an average sediment yield of 4.3 tons/ha/year due to soil erosion. Moreover, a separate study undertaken by Ullah et al. (2018) in Punjab province's Chakwal district estimated an annual average soil removal of 18 tons/ha/year. This estimation pertained specifically to areas adjacent to rivers and hilly terrain.

Gilani et al. (2022) conducted an assessment of the yearly soil erosion in Pakistan using the RUSLE model. The study provided estimations for annual soil loss across various administrative divisions of the country. According to their findings, the annual soil loss in Khyber Pakhtunkhwa is approximately 11.78 ± 21.97 tons/ha/year. In a separate study, Nasir et al. (2023) determined that the average annual soil loss from the PRB is 10.25 tons/ha/year. It's noteworthy that the Panjkora Basin, the study area, falls within the geographical scope of Khyber Pakhtunkhwa. The current study's evaluation of annual soil loss stands at 16.1 tons/ha/year, which closely corresponds to the figure reported by Gilani et al. (2022) for Khyber Pakhtunkhwa and the 10 tons/ha/year estimation for the Panjkora Basin as suggested by Nasir et al. (2023).

The PRB is selected as the study area due to its fragile soil and having shown regular soil erosion. The increasing population, deforestation and conversion of agricultural land to the built-up areas are causing continuous pressure on land resources. Soil loss has posed adverse consequences to natural ecosystems and sources of livelihood. This study holds importance as it offers insights into a region characterized by high erosion rates. This study will prove helpful to urban planners, agriculturists and conservationists to mitigate soil erosion, aiding them in implementing effective watershed and sediment control. These measures will enhance agricultural yield and prolong the functional longevity of dams and reservoirs by preventing sediment accumulation. The finding reveals that in the study area water occupied 3.09% in 1990 while in 2020 it has reduced to 1.85%. In the same way, vegetation and barren areas have significantly reduced; in 1990, the vegetation covered 63.68% of the study area while, in 2020, it has reduced to 50.89%. Like the vegetation cover, the barren land in 1990 was 17.14%, and it decreased to 14.75% in 2020. The area under snow cover was 9.74% in 1990, which reduced to 6.31% in 2020. The built-up area has increased from 1990 to 2020 from 6.34 to 26.38%. The maximum actual annual soil loss per hectare was 5,195 tons in 1990, which gradually increased to 6,247 tons in 2005. Similarly, in the study region, during 2020, the actual annual soil loss was further accelerated and marked the figure of 8,297 tons/ha/year, whereas the average annual loss of soil per hectare is also increasing. It was found from the analysis that, in 1990, the average soil loss per hectare was 9.2 tons, which gradually increased to 12.4 tons in 2005 and in 2020 it was 16.1 tons/ha/year. This demands that government line agencies and policymakers start watershed management and conservation in the basin. Many of the limitations and uncertainties in our research could be attributed to existing RUSLE formulations, including uncertainty associated with the model's simple empirical nature and numerous subcomponents; data availability issues; and the model's inability to account for soil loss due to gully erosion, mass sliding events, or the prediction of soil erosion. For large-size basins, these results are significant for baseline data preparation but for accurate and comprehensive analysis, it is recommended to use high-resolution satellite images and DEM. The sub-basin study will also lead toward more appropriate results.

The paper is not submitted to any other journal.

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

The authors declare there is no conflict.

Abidin
R. Z.
&
Abu Hassan
Z.
2005
‘ROM’ scale for forecasting erosion induced landslide risk on hilly terrain
. In:
Landslides: Risk Analysis and Sustainable Disaster Management
(Sassa, K., Fukuoka, H., Wang, F. & Wang, G., eds.)
.
Springer
,
Berlin, Heidelberg
, pp.
197
202
.
Alkharabsheh
M. M.
,
Alexandridis
T. K.
,
Bilas
G.
,
Misopolinos
N.
&
Silleos
N.
2013
Impact of land cover change on soil erosion hazard in northern Jordan using remote sensing and GIS
.
Procedia Environmental Sciences
19
,
912
921
.
Batool
S.
,
Shirazi
S. A.
&
Mahmood
S. A.
2021
Appraisal of soil erosion through RUSLE model and hypsometry in Chakwal Watershed (Potwar-Pakistan)
.
Sarhad Journal of Agriculture
37
(
2
),
594
606
.
Bizuwerk
A.
,
Taddese
G.
&
Getahun
Y.
2003
Application of GIS for Modeling Soil Loss Rate in Awash River Basin, Ethiopia
.
International Livestock Research Institute (ILRI)
,
Addis Ababa
,
Ethiopia
, pp.
1
11
.
Boardman
J.
,
Parsons
A. J.
,
Holland
R.
,
Holmes
P. J.
&
Washington
R.
2003
Development of badlands and gullies in the Sneeuberg, Great Karoo, South Africa
.
Catena
50
(
2–4
),
165
184
.
Durigon, V. L., Carvalho, D. F., Antunes, M. A. H., Oliveira, P. T. S. & Fernandes, M. M. 2014 NDVI time series for monitoring RUSLE cover management factor in a tropical watershed. International Journal of Remote Sensing 35 (2), 441–453.
Erol
A.
,
Koşkan
Ö.
&
Başaran
M. A.
2015
Socioeconomic modifications of the universal soil loss equation
.
Solid Earth
6
(
3
),
1025
1035
.
Foster
G. R.
,
Yoder
D. C.
,
Weesies
G. A.
,
McCool
D. K.
,
McGregor
K. C.
&
Bingner
R. L.
2002
User's Guide – Revised Universal Soil Loss Equation Version 2 (RUSLE 2)
.
USDA–Agricultural Research Service
,
Washington, DC
.
Fujaco
M. A. G.
,
Leite
M. G. P.
&
Neves
A. H. C. J.
2016
A GIS-based tool for estimating soil loss in agricultural river basins
.
REM-International Engineering Journal
69
,
417
424
.
Gilani
H.
,
Ahmad
A.
,
Younes
I.
&
Abbas
S.
2022
Impact assessment of land cover and land use changes on soil erosion changes (2005–2015) in Pakistan
.
Land Degradation & Development
33
(
1
),
204
217
.
GOP (Government of Pakistan)
2000
District Census Report
.
Pakistan Bureau of Statistics Government of Pakistan
,
Islamabad
.
Haile
G. W.
&
Fetene
M.
2012
Assessment of soil erosion hazard in Kilie catchment, East Shoa, Ethiopia
.
Land Degradation & Development
23
(
3
),
293
306
.
Ives
J. D.
,
Shrestha
R. B.
&
Mool
P. K.
2010
Formation of Glacial Lakes in the Hindu Kush-Himalayas and GLOF Risk Assessment
.
ICIMOD
,
Kathmandu
, pp.
10
11
Karaburun
A.
,
Demirci
A.
&
Suen
I. S.
2010
Impacts of urban growth on forest cover in Istanbul (1987–2007)
.
Environmental Monitoring and Assessment
166
,
267
277
.
Kiassari
E. M.
,
Nikkami
D.
,
Mahdian
M. H.
&
Pazira
E.
2012
Investigating rainfall erosivity indices in arid and semiarid climates of Iran
.
Turkish Journal of Agriculture and Forestry
36
(
3
),
365
378
.
Kidane
D.
&
Alemu
B.
2015
The effect of upstream land use practices on soil erosion and sedimentation in the Upper Blue Nile Basin, Ethiopia
.
Research Journal of Agriculture and Environmental Management
4
(
2
),
55
68
.
Lee
J. H.
&
Heo
J. H.
2011
Evaluation of estimation methods for rainfall erosivity based on annual precipitation in Korea
.
Journal of Hydrology
409
(
1–2
),
30
48
.
Leh
M.
,
Bajwa
S.
&
Chaubey
I.
2013
Impact of land use change on erosion risk: An integrated remote sensing, geographic information system and modelling methodology
.
Land Degradation & Development
24
(
5
),
409
421
.
Maqsoom
A.
,
Aslam
B.
,
Hassan
U.
,
Kazmi
Z. A.
,
Sodangi
M.
,
Tufail
R. F.
&
Farooq
D.
2020
Geospatial assessment of soil erosion intensity and sediment yield using the revised universal soil loss equation (RUSLE) model
.
ISPRS International Journal of Geo-Information
9
(
6
),
356
.
Nagendra
H.
,
Munroe
D. K.
&
Southworth
J.
2004
From pattern to process: Landscape fragmentation and the analysis of land use/land cover change
.
Agriculture, Ecosystems & Environment
101
(
2–3
),
111
115
.
Nasir
M. J.
,
Alam
S.
,
Ahmad
W.
,
Bateni
S. M.
,
Iqbal
J.
,
Almazroui
M.
&
Ahmad
B.
2023
Geospatial soil loss risk assessment using RUSLE model: A study of Panjkora River Basin, Khyber Pakhtunkhwa, Pakistan
.
Arabian Journal of Geosciences
16
(
7
),
440
.
Nikkami
D.
2012
Investigating sampling accuracy to estimate sediment concentrations in erosion plot tanks
.
Turkish Journal of Agriculture and Forestry
36
(
5
),
583
590
.
Panagos
P.
,
Borrelli
P.
,
Meusburger
K.
,
Yu
B.
,
Klik
A.
,
Jae Lim
K.
&
Ballabio
C.
2017
Global rainfall erosivity assessment based on high-temporal resolution rainfall records
.
Scientific Reports
7
(
1
),
1
12
.
Perović
V.
,
Životić
L.
,
Kadović
R.
,
Đorđević
A.
,
Jaramaz
D.
,
Mrvić
V.
&
Todorović
M.
2013
Spatial modelling of soil erosion potential in a mountainous watershed of South-eastern Serbia
.
Environmental Earth Sciences
68
(
1
),
115
128
.
Pimentel
D.
2006
Food and environmental threat of soil erosion
.
Journal of the Environment, Development and Sustainability
8
,
119
137
.
Pradhan
B.
,
Chaudhari
A.
,
Adinarayana
J.
&
Buchroithner
M. F.
2012
Soil erosion assessment and its correlation with landslide events using remote sensing data and GIS: A case study at Penang Island, Malaysia
.
Environmental Monitoring and Assessment
184
(
2
),
715
727
.
Rawat
K. S.
,
Mishra
A. K.
&
Bhattacharyya
R.
2016
Soil erosion risk assessment and spatial mapping using LANDSAT-7 ETM + , RUSLE, and GIS – A case study
.
Arabian Journal of Geosciences
9
(
4
),
1
22
.
Renard
K. G.
,
Foster
G. R.
,
Weesies
G. A.
&
Porter
J. P.
1991
RUSLE: Revised universal soil loss equation
.
Journal of Soil and Water Conservation
46
(
1
),
30
33
.
Renard
K. G.
&
Ferreira
V. A.
1993
RUSLE model description and database sensitivity
.
Journal of Environmental Quality
22
(
3
),
458
466
.
Renard
K. G.
,
Foster
G. R.
,
Yoder
D. C.
&
McCool
D. K.
1994
RUSLE revisited: Status, questions, answers, and the future
.
Journal of Soil and Water Conservation
49
(
3
),
213
220
.
Renard
K. G.
1997
Predicting Soil Erosion by Water: A Guide to Conservation Planning with the Revised Universal Soil Loss Equation (RUSLE)
.
United States Government Printing
.
Sharma
A.
,
Tiwari
K. N.
&
Bhadoria
P. B. S.
2011
Effect of land use land cover change on soil erosion potential in an agricultural watershed
.
Environmental Monitoring and Assessment
173
(
1
),
789
801
.
Shinde
V.
,
Tiwari
K. N.
&
Singh
M.
2010
Prioritization of micro watersheds on the basis of soil erosion hazard using remote sensing and geographic information system
.
International Journal of Water Resources and Environmental Engineering
2
(
3
),
130
136
.
Singh
G.
,
Chandra
S.
&
Babu
R.
1981
Soil loss and prediction research in India. Central Soil and Water Conservation Research Training Institute. Bulletin No T-12 D, 9, 1981
.
Tadesse
L.
,
Suryabhagavan
K. V.
,
Sridhar
G.
&
Legesse
G.
2017
Land use and land cover changes and soil erosion in Yezat Watershed, North Western Ethiopia
.
International Soil and Water Conservation Research
5
(
2
),
85
94
.
Tağıl
Ş.
2007
Land degradation risk assessment for Tuzla Creek basin (Biga Peninsula) using a GIS-based RUSLE model
.
Ecology
17
,
11
20
.
Teh
S. H.
2011
Soil Erosion Modelling Using RUSLE and GIS on Cameron Highlands, Malaysia for Hydropower Development
.
Doctoral dissertation, Master's Thesis
,
The School for Renewable Energy Science, University of Iceland
.
Ullah
S.
,
Ali
A.
,
Iqbal
M.
,
Javid
M.
&
Imran
M.
2018
Geospatial assessment of soil erosion intensity and sediment yield: A case study of Potohar Region, Pakistan
.
Environmental Earth Sciences
77
,
1
13
.
Williams
J.
1991
Search for sustainability: Agriculture and its place in the natural ecosystem. Paper presented to the Australian Institute of Agricultural Science Symposium ‘Agriculture and the Ecosystem’, Townsville, August 1990
.
Agricultural Science (Melbourne)
4
(
2
),
32
39
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).