Soil erosion and sediment yield estimation in a tropical monsoon dominated river basin using GIS-based models

Abstract

The increasing soil erosion (SE) and the associated problems for society, economy, and environment sparked a lot of interest in estimating and mapping SE at different basin scales. The estimation of SE exhibits that SE ranges from 10 to 50 t ha⁻¹yr⁻¹, with a mean SE of 20 t ha⁻¹ yr⁻¹. The very steep slopes account for 54.21% of total soil loss. The SRB areas where soil loss rates are >10 t ha⁻¹yr⁻¹ are considered the target areas which account for 27% of the study area and 96% of the soil loss). The high SY is concentrated only in the first-order basins located in a higher slope zone in the northern part of the river. Besides, basin morphometry (basin shape, relative relief) and anthropogenic activities (agricultural land) are retained in the PSLR model as significant factors contributing to SY in the entire river basin.

Keywords:

• RUSLE model
• sediment delivery ratio
• sediment yield
• sediment production rate
• partial least square regression

1. Introduction

In the twenty first century, the exposure of soils to water erosion is increasing, directly leading to soil erosion (SE), sediment yield (SY) and degradation of land (Renschler and Harbor 2002; Napoli et al. 2016; Singh and Panda 2017; Zhang et al. 2017; García-Soriano et al. 2020; Quesada-Román et al. 2021). Climate change, coupled with direct and indirect human innervations, has made the SE phenomenon a serious environmental challenge across the world. Rapid land use and land cover (LULC) changes especially in the forms of deforestation and expansion of agricultural and pasture lands are responsible for acute SE (García-Ruiz et al. 2015; Gessesse et al. 2015; Quesada-Román and Mata-Cambronero 2021; Quesada-Román and Campos-Durán 2022). Bhattacharyya et al. (2015) and Pal (2016) mentioned that ∼1964.4 Mha of global soils have faced human-induced erosion, while 1903 Mha have been exacerbated by the work of water. Consequently, the Anthropocene era has observed negative effects of accelerated SE on a variety of topics, including crop production (Balthazar et al., 2013), hydro-morphology (Thomas et al. 2018), reservoir capacity (John et al. 2022), landslides (Kirkby 1987), economic loss (Van Pelt et al. 2017; Zerihun et al. 2018), ecosystem health (Fujaco et al. 2016), and many more.

A traditional approach involves estimating the rate of soil erosion (SE) and comprehending their nature, processes, and consequences at a river basin scale. The basic basin morphometric characteristics like shape, size, and relief are explored to understand their direct and indirect influence on the nature of erosion (Das 2014; Kadam et al. 2019). Previous studies (Jose and Das 1982) developed a number of empirical models based on the geomorphological parameter for measuring soil erosion. Using morphometric criteria, Chaudhary and Sharma (1998) assessed and prioritized erosion hazards at sub-watershed level. Numerous empirical techniques, such as the Revised Universal Soil Loss Equation or RUSLE (Wischmeier and Smith 1978), and the sediment yield index (Lu et al. 2003), are widely utilized in watershed prioritization at the basin level. Garde and Kothari (1987) created an empirical equation for determining average annual sediment output considering catchment area, catchment slope, drainage density, plant cover factor, and yearly precipitation. The equation was based on data from 50 catchments in the plain region of India.

There are many numerical models for estimating soil loss. Though European SE Model (EuroSEM) (Morgan et al. 1990), Erosion Productivity Impact Calculator (EPIC), Limburg SE Model (LSEM) (Hudson 1981), Areal Nonpoint Source Watershed Environment Response Simulation (ANSWERS) (Beasley et al. 1980), Chemicals, Runoff and Erosion from Agricultural Management System (CREAMS) (Knisel 1980) are popular and highly capable models for estimating the soil loss, they are rarely found to fit them with the small-scale river basin. Contrarily, the empirical models, for example, the universal soil loss equation (USLE) and its revised version (RUSLE) are well-suited for the estimation of soil loss of small-scale river basins and therefore these models are popularly used across the world (Shamshad et al. 2008; Brooks et al. 2014; Shivhare et al. 2018; Lin et al. 2020). Since the USLE and RUSLE techniques do not account for the transportation and deposition of the eroded sediments, they only reflect the overall rate of soil erosion. Therefore, both the sediment delivery ratio (SDR) and the USLE/RUSLE concepts can be coupled to assess soil erosion/deposition and sediment yield (Jain and Das 2010; Thomas et al. 2018). This coupling approach has been well-acknowledged in recent studies for an in-depth investigation of soil erosion at different spatiotemporal scales. For example, Demir et al. (2022) estimated the soil loss of the Capakcur watershed (Turkey) by applying the RUSLE method and compared the results with the predicted SY and SDR. The watershed is estimated to have lost 96916.20 tons of soil annually, whereas the Capakcur stream transports 68656.09 tons of soil annually. A value of 0.78 was determined for the SDR.The morphometric, hydrologic, soil and LULC factors of a drainage basin have a major role in controlling the processes of soil erosion and sediment delivery. To comprehend the relative influences of these factors, several multivariate regression approaches are being applied. Though Principal Component Analysis (PCA) and Multiple Linear Regressions (MLR) are sophisticated techniques in this arena, PLSR has been most popularly used because it can handle the high inter-predictor multicollinearity (Beusen et al. 2005).

Moreover, geospatial techniques are also widely used to address the issue of basin-scale SE (Rahman et al. 2009; Aher et al. 2014; Chatterjee et al. 2014; Gayen et al. 2019). Research on multi-criteria decision-making (MCDM) methods for handling complicated decision-making has received a lot of attention in recent years. Instead of utilizing only one criterion, multiple relevant criteria are utilized in this decision-making process to optimize the model, which improves the accuracy of the decision-making process (Chitsaz and Banihabib 2015). Using several morphometric characteristics, the application of the MCDM techniques of TOPSIS, SAW, and VIKOR, aids in determining the sub-basin prioritization in a drainage basin based on the degree of susceptibility to SE (Ameri et al. 2018; Bhattacharya et al. 2020).

The humid subtropical region has been experiencing the acute problem of SE and researchers from diversified disciplines are working relentlessly to understand the dynamics of SE and develop management options. For example, Senanayake et al. (2022) estimated an increase in SE from 9.08 Mg ha⁻¹ yr⁻¹ to 11.08 Mg ha⁻¹ yr^-1 during 2000–2019 in Sri Lank and linked it to land degradation and sustainable food production. The National Bureau of Soil Survey (NBSS) estimated that ∼45% (146.8 Mha) area of India exhibits higher susceptibility to SE causing the loss of nutrients of ∼74 million tons every year correspondent to the annual loss of ∼68 billion INR (India national rupees) due to lower crop productivity (Bhattacharyya et al. 2015; Lal 2015). Eastern India is also experiencing high SE in the wake of the rapid expansion of agricultural land and deforestation. For example, Rajbanshi and Bhattacharya (2020) studied the Konar River basin located in eastern India and observed mean soil loss as 4.42 t ha⁻¹yr⁻¹. Das et al. (2022) studied the representative catchment of the Chota Nagpur Plateau located in eastern India and estimated and predicted the range of mean erosion rates for cultivated lands (26–35 t ha⁻¹yr⁻¹), built-up areas (43–48 t ha⁻¹yr⁻¹), and bare lands (30–36 t ha⁻¹yr⁻¹) during 2001–2040. Bag et al. (2022) detected 14.48% of the Sobha watershed areas undergoing severe erosion. The study area, located in eastern India, exhibits the problems of SE-induced land degradation. The SRB is principally an agriculture-dominated basin with more than 70% of its total population engaged in agrarian pursuits. Moreover, the intensive agricultural practices without proper management coupled with the accelerated rate of human activities, i.e. deforestation, settlement expansion, and other construction works have exacerbated the recent SE of the SRB. However, the river basin is completely unexplored from the standpoint of SE, land degradation and management. So far, no earlier studies have been focused on the coupling relationship between SE rate and SY in a tropical plateau fringe river basin. Thus, it would be a novel attempt to identify the morphological characteristics of each sub-basin order as well as their unique responses to SE in this small river basin. Analysing soil erosion in small basins offers crucial insights into the environmental, agricultural, hydrological, and economical effects of erosion. The management of water resources, strategies for preventing flooding, land-use planning, and actions to adapt to climate change are all influenced by this comprehension. The study would help the regional planners and other stakeholders to better manage the soil erosion of the agro-based region based on the findings of the severely susceptible SE-prone basins. Therefore, the objectives of the current investigation are (1) to estimate the sub-basin scale SE and SY using the RUSLE-SDR, and sediment production rate (SPR) models, and (2) to find out the underlying processes controlling the SE and sediment production using the PLSR model. The present study would reinstate the central argument that SE and SY decrease with increasing basin order. The novelty is that the current investigation represents a ground-breaking effort to characterize basin order-wise SE potential, and SY using an integrated novel approach.

2. Materials and methods

2.1. Study area

The Sankha River basin (SRB) is situated in eastern India, especially in the Purulia district of West Bengal (Figure 1a). The basin is also a part of the Chota Nagpur Plateau region, which is known for its undulating topography and numerous plateaus. The river originates from Ajodhya Hill (669 m above mean sea level) and joins with the Subarnarekha River at the Chandil Reservoir (23°01′58ʺ N and 86°01′21ʺ N). The river is 35 km long covering basin area 194 km² (Figure 1b). The tropical humid climate prevails in SRB where the maximum rainfall is received in the monsoon months (June to September) and the average rainfall is observed as 240 mm per year. The geo-physical setting of the basin, i.e. the steep and unstable slope, and less vegetation cover collectively lead to the development of a poor soil profile that has low resistance to external forces like the action of running water. As the basin in its catchment area also maintains high drainage density with the first order streams, therefore, the geo-physical settings provide the ideal environment for immature soil to be eroded (Mahala 2018). The LULC dynamics of the basin portray a declining trend of vegetation (55.62 km² or 28.6% in 2021 from 68.13 km² or 35.11% in 2011); cultivation land 26.79% in 2021 from 38.30% in 201; and waterbodies (0.2%) (Figure 1c and d). In contrast, the fallow land has increased by 13%, from 2011 (12.15%) to 26.14% in 2021. Overall, the basin has experienced rapid deforestation, development of settlements, and transportation networks during the last decade. Therefore, all the anthropogenic activities operating in this basin are escalating the SE process.

Figure 1. Location of the SRB and dynamics of land use and land cover. A. Location of the SRB within India, b. Drainage system of the SRB on DEM, c. LULC of the SRB for 2011, d. LULC of the SRB for 2021.

2.2. Database

Morphometric analysis of the SRB was carried out for 236 first-order, 65 s-order, 17 third-order and 4 fourth-order sub-basins (Figure 5a–d). The morphometric parameters considered for the present investigation and computed using ArcGIS software (v.10.6) are presented in Table 1. First, data from the Shuttle Radar Topography Mission (SRTM) with a resolution of 30 m were obtained from the USGS Earth Explorer (USGS 2022). Subsequently, the basin area was demarcated from the DEM and the drainage networks were extracted after the Strahler method. The results have been validated with the Survey of India (SOI) toposheet (73I/4) which has a 1:50,000 scale. In addition, rainfall grid data encompassing 0.5° of geographical resolution from the Climate Research Unit (CRU) were extracted and resampled on 30 m for preparing the SE map based on the RUSLE model (CRU 2022). The Food and Agriculture Organisation of the United Nations (FAO) produced a digital soil map of the world at a scale of 1:5.000.000, which was used to collect soil-related data, such as the amount of sand, silt, clay, and organic matter (FAO 2022). The LANDSAT TM (November 2011) and OLI (November 2021) satellite images (30 m spatial resolution), acquired from http://earthexplorer.usgs.gov, were utilized to create the LULC map of the basin. Using maximum likelihood algorithms and supervised image classification, the basin’s LULC map was created.

Table 1. Morphometric measurements of the river basin.

Parameters	Description	Reference(s)
Stream order (u)	Geometric ordering of streams according to the scale magnitude of a stream in networks	Strahler (1952)
Circularity ratio (Cr)	Rc = 4π Ab/Pb²; where, π =3.14, Ab = basin area, Pb = basin perimeter	(Miller, 1953)
Elongation ratio (Er)	Er = 2 * [(Ab)0.5/ (Dl*(π)0.5)]; where, Ab = area of basin, Dl = maximum basin length	(Schumm, 1956)
Form factor (Ff)	Ff = Ab / Dl2; where, Ab = Area of basin, Dl = maximum basin length	Horton (1932)
Compactness Coefficient (Cc)	Cc = Pb / (2 (π Ab)0.5); where, Pb = basin perimeter, Ab = basin area	(Gravelius, 1914)
Length of overland flow (LOF)	Basin CCM/2; where, CCM = constant of channel maintenance	Chorley (1969)
Relative relief (RR)	Highest relief–lowest relief	(Smith 1950)
Average slope (As)	As= (tan^−1 ((total length of contours x contour interval) / area)	Horton (1945)
Drainage density (DD)	DD = Total length of streams within a basin / area	Horton (1945)
Stream frequency (SF)	SF = Total number of streams within a basin / area	Horton (1945)
Drainage texture (DT)	DT = Basin SF × basin DD	Smith (1950)
Ruggedness index (RI)	RI = (relative relief of basin × drainage density of basin) / 1000	Horton (1945)
Dissection index (DI)	DI = (Relative relief / Maximum elevation of basin)
Hypsometric integral (HI)	HI = $\left\{\sum \text{Xi}\mathrm{(}\text{yi + 1}\mathrm{)}\sim \text{yi}\mathrm{(}\text{xi + 1}\mathrm{)}\right\}$/2 Where xi is the proportion of the basin’s area (Xi = a/A) and yi is the proportion of the basin’s height (Yi = h/H). Once more, h = cumulative relative height and H = basin relative height; a = cumulative area; A = total basin area; and h = basin relative height	Strahler (1952)

2.3. Methods of data processing

Figure 2 depicts the current study, which was carried out using a sound methodology. The first section pertains to the delineation of the watershed area along with numbers, orders, and lengths of streams. The second section covers the morphometric measurements of the river basin (Table 1). These measurements lead to an understanding of the hydro-morphological pattern and behaviour of the sub-watersheds and, of course, the nature, mechanism and magnitude of SE in the respective sub-watersheds. Moreover, in the third section, an attempt has also been made to apply the RUSLE model for SE potential followed by SDR and SY models for comparison with the SPR. The next section is concerned with the application of Jose and Das’s (1982) morphometry-based experimental model for quantifying the SPR of the respective sub-watersheds. Finally, the nature of SE of the entire basin has been studied dealing with the underlying processes triggering the SE in varying magnitude.

Figure 2. Methodological flow chart representing the processes followed for the present investigation. The process starts with building the RUSLE model database. Following that, the model’s outcome was validated. The PLSR method has been used to discover probable causes of soil erosion if the results accurately reflect reality. In the event that the validation fails to accurately reflect reality, the RUSLE database has been examined and improved.

2.4. SE estimation by RUSLE model

The model is efficiently used for estimating soil loss of a specific region based on five basic model parameters as presented in Equation (1)(1) $$\mathit{\text{SE}}=R\times K\times L\times S\times C\times P$$(1) . (1) $$\mathit{\text{SE}}=R\times K\times L\times S\times C\times P$$(1)

Where SE is denoted as estimated average SE potential (t ha⁻¹yr⁻¹), R as a factor of rainfall erosivity measured in MJ mm ha⁻¹ h⁻¹ yr⁻¹, K as a factor of soil erodibility measured in t h M J⁻¹ mm⁻¹, LS as a factor of dimensionless slope length and steepness, C as a factor of land use specific cover management, and P as conservation support practice factor.

The degree of uncertainty in the model depends on the selection of the parameters as model inputs. DEM is one of the important factors that, together with rainfall and soil erodibility parameters, is likely to induce uncertainty in predictions of soil erosion. When estimating soil loss using the RUSLE model, the SRTM-DEM-derived parameters performed better than the Cartosat-DEM and ASTER-DEM generated values and comparable results have been seen in many river basins of India (Mondal et al. 2017; Pandey et al. 2021). Soil erosion estimates based on SRTM-DEM are thus expected to have less uncertainty. Particularly, the total energy of the rainfall and the maximum 30-minute rainfall intensity (I30) of storm events are rarely available in developing countries, which raises the uncertainty of the conclusions of the model. However, for India’s climatic zones, soil erosion estimates utilising the correlation between the annual rainfall and EI30 values (a rainfall erosivity index based on 30-minutes rainfall intensity) proposed by Singh et al. (1981) produced better findings. Therefore, the same has been followed for deriving the R-factor of SRB. Thus, uncertainties in the present computation tend to be minimal.

2.4.1. Derivation of RUSLE parameters

The derivations of the model parameters followed for the present study have been illustrated in the following sections.

2.4.1.1. Rainfall erosivity factor (R)

R-factor indicates the likelihood of erosion capability of rainfall which is entirely dependent upon the physical characteristics of rainfall e.g. quantity and intensity of rainfall, size and terminal velocity of the raindrop, and kinetic energy (Behera et al. 2020). However, the lack of the above mention parameters in the rainfall dataset leads to the adoption of an indirect method for calculating the R-factor (Diodato and Bellocchi 2007). Equation (2)(2) $$R=\sum _{i=1}^{12}1.735\times {10}^{(1.5\times \mathrm{\text{log}}\hspace{0.17em}\left(\frac{P_i^2}{P}\right)-0.08188)}$$(2) suggested by Wischmeier and Smith (1978) and later on revised by Arnoldus (1980) has been considered for the R factor. The method is widely applied for estimating soil erosion of the small (Prasannakumar et al. 2012) and large river basins (Ganasri and Ramesh 2016) located particularly in the tropical region. Therefore, SRB located in the tropical region is also suited to this method. (2) $$R=\sum _{i=1}^{12}1.735\times {10}^{(1.5\times \mathrm{\text{log}}\hspace{0.17em}\left(\frac{P_i^2}{P}\right)-0.08188)}$$(2) where Pi and P, respectively, stand for the monthly and yearly rainfall (mm). For the present study, monthly time series of rainfall data has been used from 1990 to 2021.

2.4.1.2. Soil erodibility factor (K)

It represents the potentiality of the intrinsic erodibility of soil particles by the means of separation and transit of soil particles. Therefore, K-factor is derived considering the physicochemical properties of the soil (top layer) of the studied region. In the RUSLE model the value of K-factor is derived followed by the equation of Wischmeier and Smith (1978). In this method the value of K is derived using the percentage of sand, silt, clay and organic matter present in every soil type (Table 2). This method is widely used in the tropical region, predominately in India (Chatterjee et al. 2014; Mahala 2018). Therefore, we have considered the same methods for estimating this factor (Equation 3–4). (3) $$K=0.1317\times \frac{0.00021\times M^{1.14}\times \left(12‐\mathrm{\text{OM}}\right)+3.25\left(c_{\mathrm{\text{soilstr}}}-2\right)+2.5(c_{\mathrm{\text{perm}}}-3)}{100}$$(3) where M is the textural factor computed using Equation (4)(4) $$M=(\mathrm{\%}{m}_{\mathit{\text{silt}}}+{\mathrm{\%} m}_{\mathit{\text{vfs}}})\times (100-{\mathrm{\%} m}_c)$$(4) . (4) $$M=(\mathrm{\%}{m}_{\mathit{\text{silt}}}+{\mathrm{\%} m}_{\mathit{\text{vfs}}})\times (100-{\mathrm{\%} m}_c)$$(4) where OM represents the soil organic matter (%); m_silt is the percentage fraction content of silts (0.002–0.5 mm), m_vfs is the percentage fraction content of fine sand (0.05–0.1 mm), m_c is the percentage fraction content of clay, cs_oilstr is soil structure code given by USDA, c_perm is profile permeability class of USDA.

Table 2. K factor values of soil of SRB.

Soil unit symbol	Sand (%)	Silt (%)	Clay (%)	OM (%)	K-factor
LF	74.6	9.6	15.9	0.39	0.151112
ND	38.9	17.6	43.6	1.57	0.11027
I	58.9	16.2	24.9	0.97	0.140413

2.4.1.3. Slope length factor (L)

The L factor denotes the linear distance measured between the initiation of overland flow and the point where the slope is imperceptible leading to the depositional process. Measurement of the L factor in the entire region is both costly and time-consuming. Therefore, remote sensing as an indirect method is widely used for this measurement. Desmet and Govers (1996) first suggested DEM-based measurement for the USA condition. The present study also followed the remote sensing approach and applied Equation (5)(5) $$L_{\mathit{\text{ij}}}=\frac{({A_{j,i}+D^2)}^{m+1}-A_{i,j-\mathit{\text{in}}}^{m+1}}{D^{m+2}\times x_{i,j}^m\times {22.13}^m}$$(5) for the derivation of the L factor. (5) $$L_{\mathit{\text{ij}}}=\frac{({A_{j,i}+D^2)}^{m+1}-A_{i,j-\mathit{\text{in}}}^{m+1}}{D^{m+2}\times x_{i,j}^m\times {22.13}^m}$$(5) where L_ij is the L factor having the coordinates of i and j. A_ij is the area contributed to the inlet of the respective grid cell with their coordinates (i, j) measured as square metre, and D is the cell size of the respective grid. X_i,j = sin a_i,𝑗 + cos a_i,𝑗, where the a_i,𝑗 indicates the aspect direction of the respective grid cell (i,j), and m is the ratio of β representing the rill to inter-rill erosion and calculated using Equations (6)(6) $$m=\frac{\beta }{\beta +1}$$(6) and (7)(7) $$\beta =\frac{\frac{\mathrm{\text{sin}}\theta }{0.0896}}{[0.56+3\times {\left(\mathrm{\text{sin}}\theta \right)}^{0.8}]}$$(7) . (6) $$m=\frac{\beta }{\beta +1}$$(6) (7) $$\beta =\frac{\frac{\mathrm{\text{sin}}\theta }{0.0896}}{[0.56+3\times {\left(\mathrm{\text{sin}}\theta \right)}^{0.8}]}$$(7) where θ represents the slope angle in degrees. The m has a range of 0 to 1, and it gets close to 0 when the ratio of rill to inter-rill erosion is almost zero.

2.4.1.4. Slope steepness factor (S)

The empirical findings of McCool et al. (1987) work confirmed faster soil loss with a slope steeper than 9%. Following this observation, Renard (1997) adopted an algorithm for computing the S factor for RUSLE (Equations 8(8) $$S=10.8\times \mathrm{\text{sin}}\theta +0.03$$(8) and 9(9) $$S=16.8\times \mathrm{\text{sin}}\theta -0.5$$(9) ). (8) $$S=10.8\times \mathrm{\text{sin}}\theta +0.03$$(8) (9) $$S=16.8\times \mathrm{\text{sin}}\theta -0.5$$(9) where θ represents the slope gradient measured in degrees. Moreover, Equation (7)(7) $$\beta =\frac{\frac{\mathrm{\text{sin}}\theta }{0.0896}}{[0.56+3\times {\left(\mathrm{\text{sin}}\theta \right)}^{0.8}]}$$(7) is applied when the slope is <0.09 whereas Eq. has been fitted for the slope ≥ 0.09. The ‘L’ and ‘S’ factors jointly represent the topographic factor that has a direct influence on the process of SE.

2.4.1.5. Crop cover management factor (C)

The ratio of soil loss under a particular crop to that of the base soil is represented by the cover and management factor (Morgan et al. 1990). The C-factor indicates how different factors, such as crops, different soil covers, availability of subsurface bio-mass, etc., affect the rate of SE Kouli et al. (2009) and Zhou et al. (2008). Since plot-level identification of crop types with course resolution satellite images is becoming less useful, NDVI has been employed in many research as an alternative for determining C-factor (De Jong et al. 1999; Zhou et al. 2008; Kouli et al. 2009). By lowering the energy of the raindrops during rainfall, vegetation cover stops soil erosion. Nonetheless, the types of vegetation and their density determine how much energy is absorbed by the vegetation. Therefore, NDVI-based C-factor has been determined using (Equation 10(10) $$\mathrm{C}=\mathrm{\text{exp}}[-\alpha \times \frac{\mathit{\text{NDVI}}}{\left(\beta -\mathit{\text{NDVI}}\right)}]$$(10) ). (10) $$\mathrm{C}=\mathrm{\text{exp}}[-\alpha \times \frac{\mathit{\text{NDVI}}}{\left(\beta -\mathit{\text{NDVI}}\right)}]$$(10) where α and β denote the shape parameters used to determine the shape of the NDVI-C curve where α and β are considered as 2 and 1, respectively, according to Van der Knijff et al. (2000).

2.4.1.6. Conservation practices factor (P)

The ratio of soil loss from a field with a certain conservation technique to one where none is used is known as the support practice or P-factor. According to Morgan et al. (1990), the value of P is 1.0 in the absence of conservation measures. By altering the concentration and velocity of runoff, conservation strategies aim to reduce the likelihood of erosion (Renard 1997). Therefore, the P-factor values are derived using information about the slope steepness and soil management practice of SRB. It is pertinent to mention here that the land use map has been prepared using the Landsat OLI images (2021) processing while the slope classes have been prepared by ASTER DEM processing. Then, different weights have been assigned to different slope classes and management strategies. The management strategies that significantly lessen SE have been given higher weights, and vice versa as proposed by USDA-SCS (1972). Similarly, different weights are assigned to different slope classes, with higher weights assigned to those reducing soil erosion (Behera et al. 2020) (Table 3). Finally, in ArcGIS 10.6, a raster overlay process has been used to build the P-factor map integrating the land use and slope maps. The P value ranges from 0 to 1, with values close to 1 indicating areas more vulnerable to erosion (Chatterjee et al. 2014).

Table 3. P factor values of different land-use categories of SRB.

	Slope (%)	P-factor
Cultivated land	0-5	0.1
	5–10	0.12
	10-20	0.14
	20-30	0.19
	30-50	0.25
	>50	0.7
Waterbody	0–100	0
Settlement	0–100	0.05
Vegetation cover	0-5	0.03
	5–10	0.05
	10-20	0.1
	31-50	0.2
	>50	0.5
Fallow land	0-5	0.1
	5–10	0.13
	11-20	0.15
	21-30	0.2
	30-50	0.4
	>50	0.75

2.5. Estimation of SDR and SY

All the individual RUSLE parameters are computed and developed as thematic layers with the help of the ArcGIS 10.6 software. Afterward, all the layers are multiplied using the ‘raster calculator’ to obtain the output the model. The model portrays the spatial variation of the SE map of the SRB. However, it is more important to compute the actual amount of soil leaving the catchment (SY). Therefore, many researchers have already addressed this point and attempted to quantify SDR-based SY in the Indian context following the recommendation of USDA SCS (Vigiak et al. 2012; Behera et al. 2020; Rajbanshi and Bhattacharya 2020). USDA SCS recommendation for SY estimation is very popular, reliable, and simple to address (Behera et al. 2020; Rajbanshi and Bhattacharya 2020). Therefore the present study has followed the same for estimating the SDR and SY (Equations 11(11) $$\mathit{\text{SDR}}=0.51\times A^{-0.11}$$(11) And 12(12) $$\mathit{\text{SY}}=\mathit{\text{SDR}}\times \mathit{\text{SE}}$$(12) ). (11) $$\mathit{\text{SDR}}=0.51\times A^{-0.11}$$(11) (12) $$\mathit{\text{SY}}=\mathit{\text{SDR}}\times \mathit{\text{SE}}$$(12) where A denotes the order-wise area of the SRB (km²).

2.6. Estimation of SPR

Many researchers have studied drainage basin morphometric analysis to understand the complexity of the erosion process mechanisms (Kadam et al. 2019; Kulimushi et al. 2021). For the present study, Jose and Das (1982) morphometry-based experimental model has been followed to estimate the capacity and actual sediment production of the SRB. The model was developed based on the empirical relationship involving specific morphometric parameters such as form factor (Ff), Circularity ratio (Rc), and Compactness coefficient (Cc) (Kulimushi et al. 2021). The reasons behind adopting this method are (i) it is suitably applied across the world to a small-scale (Kadam et al. 2019; Sofi et al. 2021) and large-scale (Kulimushi et al. 2021) river basin for estimating SPR, and (ii) the parameters involved in this method can be measured easily. Moreover, as Jose and Das’s (1982) morphometry-based experimental model is very popular in India, the present study also applied it to every stream order (first to fourth) of the SRB (Equation 13(13) $$\mathrm{\text{Log}}\left(\mathrm{\text{SPR}}\right)=4919.80+48.64\{\mathrm{\text{Log}} (100+\mathrm{\text{Ff}}\left)\right\}-1337.77\{\mathrm{\text{Log}} (100+\mathrm{\text{Rc}}\left)\right\}-1165.65\{\mathrm{\text{Log}} (100+\mathrm{\text{Cc}}\left)\right\}$$(13) ). (13) $$\mathrm{\text{Log}}\left(\mathrm{\text{SPR}}\right)=4919.80+48.64\{\mathrm{\text{Log}} (100+\mathrm{\text{Ff}}\left)\right\}-1337.77\{\mathrm{\text{Log}} (100+\mathrm{\text{Rc}}\left)\right\}-1165.65\{\mathrm{\text{Log}} (100+\mathrm{\text{Cc}}\left)\right\}$$(13) where Ff indicates the form factor, Rc represents the circularity ratio, and Cc represents the compactness coefficient (Table 1).

2.7. Estimation of runoff rate

The amount of water that is transported via rivers from their sources to other regions is measured as run-off. It was reliant on the rates of precipitation, evapotranspiration, and percolation, each of which was regulated by lithology, relief, climate, and other factors. Because of their low relief, gentle slopes, and consequently good terrain permeability, plain areas typically have low run-off rates, whereas hilly areas typically have high run-off rates because of their high relief as well as steep slopes. The run-off rate of the SRB has been using the Jose and Das (1982) model, which is represented by Eq. 14(14) $$\mathrm{\text{Log}}\left(\mathrm{Q}\right)=2238.43+22.12\{\mathrm{\text{Log}} (100+\mathrm{\text{Ff}}\left)\right\}-608.28\{\mathrm{\text{Log}} (100+\mathrm{\text{Rc}}\left)\right\}-530.02\{\mathrm{\text{Log}} (100+\mathrm{\text{Cc}}\left)\right\}$$(14) . The method followed uses logarithmic functions and variables like Ff, Rc, and Cc, and the availability and accuracy of data for these variables influence the choice of the method. (14) $$\mathrm{\text{Log}}\left(\mathrm{Q}\right)=2238.43+22.12\{\mathrm{\text{Log}} (100+\mathrm{\text{Ff}}\left)\right\}-608.28\{\mathrm{\text{Log}} (100+\mathrm{\text{Rc}}\left)\right\}-530.02\{\mathrm{\text{Log}} (100+\mathrm{\text{Cc}}\left)\right\}$$(14)

2.8. Investigating the influencing factors of erosion using partial least square regression (PLSR)

Partial Least Squares Regression (PLSR) is a statistical technique commonly used for analyzing the relationships between a set of independent variables and a dependent variable, especially when there are high collinearity or multicollinearity issues (Beusen et al. 2005). The advantages of principal component analysis (PCA) and multiple linear regression (MLR) are combined in PLSR. (Zhang et al. 2011). Therefore, PLSR has been applied to identify the influencing factors of soil erosion. The PLS tool from OriginPro software (v. 2016) was used for the present investigation. Therefore, a total of 26 variables under morphometric, LULC, climate and soil characteristics are considered as the predictors of sediment yield at every sub-basin order of SRB. Finally, PLSR determines which independent variables (predictors) have a significant influence on the dependent variable (erosion) by considering their contributions to both the dependent and independent variance.

3. Results and discussion

3.1. Estimation of SE potential using RUSLE

3.1.1. Spatiality of RUSLE parameters

Figure 3 illustrates the regional variation in the total amount of rainfall (rainfall depth) and its erosivity for the subbasins of the SRB. The rainfall erosivity value correlates directly and linearly with the magnitude of rainfall depth for the majority of the places showing the greatest value. In Figure 3a, an inverse distance weighting (IDW) approach to interpolation has been used to create a map showing the spatial distribution of the R-factor using the point rainfall erosivity data. The mean reading of the rainfall erosivity (measured in MJ mm ha⁻¹ h⁻¹yr⁻¹) is found to be 738.76 with a standard deviation (SD) of 2.58 and a range of 733.39 to 746.21. The highest value is restricted to the basin’s lower part and while the minimum value is recorded in the upper part of the basin. The spatial variation of the K-factor (measured in t ha MJ⁻¹ ha⁻¹ mm⁻¹) is portrayed in Figure 3b. The mean concentration of the K-factor is recorded as 0.13. The upper part of the SRB exhibits an elevated concentration of the K-factor (0.15) and while it is low in the lower portions of the basin (0.11). AS the SRB is situated on the fringe of the Chotonagpur Plateau and has a tropical humid environment, it is the ideal place for the formation of lithosol soil, which makes up 85% of the study area. The little variation in K-factor values is caused by the uniform geo-climatic conditions of the small river basin. The L and S factors were originally calculated separately, and the topographic factor (LS-factor) was generated by multiplying the LS-factor. The spatial pattern of the LS-factor shown in Figure 3c varied from 0.00 for level terrain to 167.13 for steep slopes. Furthermore, the LS-factor, which has values between 1 and 5, covers more than 80% of the catchment area, suggesting that the topography of the basin is not as variable. The spatial pattern of the C-factor exhibited in Figure 3d ranged from 0.27 recorded in the lower part of the basin to 0.58 located in the upper part of the basin. The P-factor relates to the decreased SE risk brought on by the protective management practices or structures implemented in a particular area. Therefore, the spatial variation of the P-factor was created by incorporating these data into the raster calculator (Figure 3e). The mean P-factor value was revealed to be 0.12 with a range between 0 and 0.75, with 0 representing locations with the lower SE potential as a result of adequate adaptation measures and 0.75 denoting extremely erosion-prone regions as a consequence of a lack of appropriate management approaches.

Figure 3. Spatial pattern of RUSLE parameters: (a) R factor, (b) K factor, (c) LS factor, (d) C factor, and (e) P factor.

3.1.2. RUSLE based-SE analysis

3.1.2.1. SE analysis of different erosion severity classes

We employed the raster calculator of the ArcGIS software (v. 10.1) to overlay the six components of RUSLE to calculate the average SE. The annual SE for the entire basin is shown in Figure 4 as a spatial map that shows the variation in SE at the pixel level for the entire river basin. The mean SE is 20 t ha⁻¹yr⁻¹ ranging from 10 t ha⁻¹yr⁻¹ to 50 t ha⁻¹yr⁻¹. Additionally, SE values have been classified into 5 categories – very low erosion area (<10 t ha⁻¹yr⁻¹), low (10–20 t ha⁻¹yr⁻¹), moderate (20–30 t ha⁻¹yr⁻¹), high (30–40 t ha⁻¹yr⁻¹) and very high (>40 t ha⁻¹yr⁻¹). In detail, SE < 10 covers 89% of the SRB while SE > 40 covers only 3%. Moreover, 1.1% of the SRB falls into the category of moderate SE risk, while 2.3% falls into the category of low SE risk.

Figure 4. SE potential severity of the Sankha river basin.

3.1.2.2. Order-wise SE analysis

The SE has also been computed order-wise in t ha⁻¹yr⁻¹. Regarding the first-order basin, the mean SE has been observed as 173 t ha⁻¹yr⁻¹ with the range: 0–1792; standard deviation (SD): 280. Out of 236 first order river basins, 31,44,30, 13, and 118 basin record SE in the categories of very low with SE <10 t ha⁻¹yr⁻¹ (mean: 6.41 t ha⁻¹yr⁻¹); low SE: 10-20 (mean: 14.35 t ha⁻¹yr⁻¹); moderate SE: 20-30 t ha⁻¹yr⁻¹ (mean: 24.81 t ha⁻¹yr⁻¹); high SE: 30-40 t ha⁻¹yr⁻¹ (mean: 33.32 t ha⁻¹yr⁻¹), and very high SE> 40 t ha⁻¹yr⁻¹ (mean: 330.31 t ha⁻¹yr⁻¹), respectively. This classification is made based on the rage of SE values and following the classification of SE done by others researches in the Indian context and adjacent to our study area (Dutta et al. 2015; Samanta et al. 2016; Behera et al. 2020) For the classification of SE, A total of 3, 6, 8, and 48 second order subbasins, respectively, exhibit low (mean: 15.27), moderate (mean: 25.91), high (mean: 35.32), and very high (mean: 440.32) SE. Furthermore, only one basin falls into the high (30-40 t ha⁻¹yr⁻¹) category, while the rest fall into the very high (SE > 40 t ha⁻¹yr⁻¹, mean: 816.94) category, while one basin falls into the very low (SE <10 t ha⁻¹yr⁻¹) and moderate (SE: 20-30 t ha⁻¹yr⁻¹) categories, while the other two fall into the very high (SE> 40 t ha⁻¹yr⁻¹) category.

3.1.2.3. SE analysis of different slope classes and LULC types

The classes of very gentle and gentle slopes make up about 84% of the basin’s total surface. The average SE from very gentle slopes is 2.22 t ha⁻¹yr⁻¹ with an SD of 7.08 while the average soil loss rate from gentle slopes is 5.38 t ha⁻¹yr⁻¹ with an SD of 12.79. Moreover, both slope classes account for 1.63% of the total soil loss. Furthermore, very steep slopes account for about 2.5% of the total basin area (Table 4) and mainly concentrate in the north portion of the river basin. The very steep slopes account for 54.21% of total soil loss (mean: 267.43 t ha⁻¹yr⁻¹, SD: 249.21), whereas steep slopes (5.1% of the total SRB) account for 35.40% of total soil loss (mean: 174.63 t ha⁻¹yr⁻¹, SD: 185.25). Therefore, very steep and steep slopes jointly contribute the lion’s share (71%) of erosion of the basin. Besides, the moderate slopes mainly concentrate in the southern portion of the river basin (7.5% of the total SRB) accounting for 8.75% of the total soil loss (mean: 43.16 t ha⁻¹yr⁻¹, SD: 72.42) (Table 4). The forest land accounts for around 28% of the total basin area followed by agricultural land (27%), fallow land (26%), settlement (13%), and waterbody (4%) (Table 4). In the forest land area of the SRB higher soil loss is recorded as high. The increased slope steepness may be the key cause of the higher soil loss rate in the forest. Furthermore, mean soil loss follows the same pattern, with forest land accounting for 58.39%, fallow land accounting for 11.97%, settlement accounting for 8.55%, and agricultural land accounting for 4.06%. The contrasting observations are commonly recorded that human interventions like agricultural land intensify high rates of soil erosion (Rahman et al. 2009; Farhan and Nawaiseh 2015; Mukanov et al. 2019). For example, Ganasri and Ramesh (2016), Nethravathi Basin, India, observed that due to an increase in agricultural area, the gross quantity of erosion had increased to 488,012.6 t yr⁻¹ from 473,339 t yr⁻¹. However, the observation similar to our study has also been recorded by others. For examples, Nyesheja et al. (2019) studied Congo Nile Ridge region of Rwanda and observed that forestland had the highest soil loss rate of 104.93 t ha⁻¹yr⁻¹, compared to the agricultural land 19.10 t ha⁻¹yr⁻¹. They also argued that elevated terrain with steep slopes causes the high soil erosion in the forest area.

Table 4. Erosion potential classes and their distribution based on the LULC and slope categories.

Classes		Erosion potential
		Mean (t ha^−1 yr^−1)	Area (%)
Basin specific	Very low	6.30	89.28
	Low	14.41	4.28
	Moderate	25.59	2.34
	High	34.94	1.13
	Very high	398.28	2.95
Land use specific	Forest	11.97	28.6301
	Agriculture	4.06	26.93
	Fallow lands	8.55	26.00
	Settlement	58.39	4.96
Slope specific	Very gentle	2.22	13.46
	Gentle	5.83	52.50
	Moderate	43.16	32.15
	Steep	174.63	5.16
	Very steep	267.43	2.67

3.1.3. RUSLE model findings and other similar studies

Soil erosion models try to approximate the reality, however, there may be some degree of uncertainty in these models. Even though uncertainties related to the model are not shown to be larger than some of the complicated process-oriented models, soil erosion predictions from SRB need to be validated. The outcomes of research related to the study region conducted by others considering similar geo-climatic variables have been used to validate our findings. The estimated rate of soil erosion in the SRB in comparison with other river basins in India with similar geo-climatic conditions indicates a good fit (70–85%). For example, Rajbansi and Bhattacharya assessed SE in the Konar River basin located in a similar geo-climatic setting (physiography, temperature, soil, and vegetation) with human interventions (land use patterns) that closely correspond to the SRB. Their estimate of 4.42 t ha⁻¹ yr⁻¹ for the mean SE rate represents parallel nature of soil erosion. Samanta et al. (2016) used the RUSLE model to study SE along the Subarnarekha River, a distributary of the Sankha River. They found that 26% of the basin area belonged to the very low SE class (less than 5 t ha⁻¹ yr⁻¹), followed by 20% of moderate erosion (between 10 and 20 t ha⁻¹ yr⁻¹). Their findings concur with those of the current study with a similarity of ∼80%. A similar model was used by Mahala (2018) on the Kangsabati River basin (adjacent to SRB) and found that 82% of the basin area experienced low SE with SE 2.5 t ha⁻¹ yr⁻¹. Besides, a large portion of the Kangsabati River basin falls into the low (0.57-0.83 mt ha⁻¹) and medium (0.84–1.08 mt ha⁻¹) SY classes, according to Bhattacharya et al. (2020). Another comparable study conducted by Dutta et al. (2015) on the Sanjal watershed of the Kharkai River sub-basin (in the West Singbhum district of Jharkhand and near the Sankha River) explained that 49 of the basin area fell into the 0–5 t ha⁻¹ yr⁻¹ class. Therefore, their estimation supports the validity of our findings.

3.2. SDR and SY analysis

The SDR values varied from 0.44 to 0.74 (mean: 0.60, SD: 0.05) for first order, from 0.41 to 0.63 (mean: 0.51, SD: 0.04) for second order, from 0.38 to 0.51 (mean: 0.42, SD: 0.04) for third order, and from 0.36 to 0.39 (mean: 0.37, SD: 0.01) for fourth order. Therefore, value indicates SDR in SRB in general for all sub-basin order. However, it is classified into high, moderate and low based on the range of SRD. The high SDR values of first-order basins have been found in hillslope zones, the northern part of SRB (Figure 5a). In 100 basins of the 236 first-order river basins, the SDR was observed high (SDR > 0.62 and above) with a mean 0.65. besides moderate SDR (0.54-0.62) was observed for 111 first-order subbasins scatted in the western part of the basin. Low SDR (<0.54) was recorded for 25 first-order subbasins (mean= 0.51) scattered in the central part of the basin. Regarding second order basins, 18 basins exhibit high SDR (mean: 0.63) located in the northern part of SRB. In contrast, moderate SDR (mean: 0.55) has been recorded for 16 basins located in the north-western portion of SRB. Additionally, 31 basins show the lowest SDR (mean: 0.47) located in the eastern and central parts of SRB (Figure 5b). Additionally, the entire third order (mean: 0.42) and fourth order (mean: 0.37) sub-basins exhibit very low to low SDR (Figure 5c–d). The SDR values in the Brahmani River basin, which is adjacent to SRB, were discovered to be unevenly distributed and to range from 0.22 to 0.54. A SDR value ranging from 0.22 to 0.32 included roughly 95% of the basin. It was discovered that the minimum basin area was below the 0.32–0.42 SDR range.

Figure 5. Sub-basin order-wise estimation of sediment delivery ratio of the Sankha river. (a) SDR of first order basin, (b) SDR of second order basin, (c) SDR of third order basin, and (d) SDR of fourth order basin.

The values of SY in SRB are classified into six classes: very low (<10 t ha⁻¹yr⁻¹), low (10-20 t ha⁻¹yr⁻¹), moderate (20-40 t ha⁻¹yr⁻¹), moderately high (40-80 t ha⁻¹yr⁻¹), high (80–160 t ha⁻¹yr⁻¹), and very high classes (>160 t ha⁻¹yr⁻¹). Figure 6a depicts the spatial variation of SY of 236 first-order sub-basins. The first-order subbasins have an average SY of 102.37, with 64, 45, 24, 60, 32, and 55 basins falling into the very low (mean: 5.96 t ha⁻¹yr⁻¹), low (mean: 14.14 t ha⁻¹yr⁻¹), moderately low (mean: 27.6 t ha⁻¹yr⁻¹), moderately high (mean: 60.06 t ha⁻¹yr⁻¹), high (mean: 118.3 t ha⁻¹yr⁻¹ 1), and very high (mean: 322.64 t ha⁻¹yr⁻¹) SY classes. The second-order subbasins have an average SY of 169.79, with 3, 15, 13, 15, 7, and 23 basins falling very low (mean: 8.51 t ha⁻¹yr⁻¹), low (mean: 15.37 t ha⁻¹yr⁻¹), moderately low (mean: 24.47 t ha⁻¹yr⁻¹), moderately high (mean: 57.53 t ha⁻¹yr⁻¹), high (mean: 113.44 t ha⁻¹yr⁻¹), and very high (mean: 409.99 t ha⁻¹yr⁻¹) SY classes (Figure 6b). The third-order subbasins have an average SY of 87.89 t ha⁻¹yr⁻¹. The SY variation indicates that 2, 1, 1, 2, and 11 basins fall into low, moderate, moderately high, high, and very high classes (Figure 6c). The fourth-order sub-basins exhibit relatively low SY recording the average 12.08 t ha⁻¹yr⁻¹ where 2 located in the northern part of the SRB are exhibiting relatively hight SY (Figure 6d). As water and sediment move downstream, they often encounter larger and deeper channels that allow for sediment deposition. Higher order sub-basins may have more significant sediment storage areas, which can reduce the overall sediment yield.

Figure 6. Sub-basin order-wise estimation of SY (t ha⁻¹year⁻¹) of the Sankha river. (a) SY of first order basin, (b) SY of second order basin, (c) SY of third order basin, and (d) SY of fourth order basin.

3.3. SPR analysis

The value of SPR of the study basin has been computed by applying the Jose and Das (1982) method which combined the effects of form factor, compactness coefficient, and circulatory ratio in sediment production. As a result, before performing SPR analysis, the variation of the factors was investigated across the entire basin. The elongated basin generally exhibits lower Ff with low peak flow for a longer duration. For example, the Shivganga watershed (part of the Western Ghat region of Maharashtra state, India) has Ff value of 0.42, which indicates that it is a sub-circular to the somewhat elongated basin with low peak flows for a long duration (Magesh and Chandrasekar, 2016). For the present study, the mean values of Ff at first, second, third and fourth order basins are observed as 0.46 (range: 0.048-4.50), 0.40 (range: 0.18-0.78), 0.45 (range: 0.21-0.91) and 0.39 (range:0.32-0.48) respectively. Watersheds with low form factors, or 0.14, exhibit their elongated shape with streams running for longer durations and are thus least prone to SE. Besides, the high form factors indicate high peak flows of shorter duration and consequently high erosion. Generally, a basin looks circular in shape when Cc = 1. The basin deviates from its circular shape when Cc > 1. (Altaf et al. 2014). For the present study, the mean values of Cc at first, second, third and fourth order basins are observed to be 1.21 (range:1.02–1.71), 1.24 (range: 1.04–1.82), 1.33 (ranger:1.07–1.61) and 1.45 (range:1.26–1.65) respectively. Moreover, the mean values of Cr at first, second, third and fourth order basins are observed as the 0.71 (range:0.34-0.96), 0.66 (range: 0.30-0.92), 0.58 (ranger:0.38-0.86) and 0.48 (range:0.36-0.62) respectively. The inverse relationship between the compactness coefficient and circularity ratio arises because the metrics are essentially measuring opposite aspects of shape. As a shape becomes more compact (less elongated) and approaches a circular form, the compactness coefficient decreases, but the circularity ratio increases (Gravelius, 1914; Miller 1953). Conversely, as a shape becomes more elongated and deviates from a circular form, the compactness coefficient increases, but the circularity ratio decreases.

The erosion rate for the first-order basin ranges from 0.48 to 15.75 hm/100 km²/yr with an average value of 1.66 hm/100 km²/yr and an SD of 1.22. The first-order streams distributed in the upper portions of the basin portray a high erosion (Figure 7a). Besides, the moderate erosion zone was found in the middle sections of the basin. The SPR values are recorded as high in 4 subbasins (mean: 8.6 hm/100 km²/year), and moderate in 79 subbasins (mean: 2.18). and low in 153 subbasins (mean: 1.21 hm/100 km²/yr). The spatial variation of the erosion rate at the second-order sub-basin is shown in Figure 7b. With an average value of 1.68 hm/100 km²/yr and an SD of 0.54, the second-order basin’s erosion rate ranges from 0.26 to 2.44 hm/100 km²/yr. The upper part of the basin’s second-order streams shows a high potential for erosion. Additionally, the centre portions of the basin had a moderate erosion potential zone. In 33 subbasins, the likelihood of SE is high (mean: 2.12 hm/100 km²/yr), in 14 subbasins it is moderate (mean: 1.59 hm/100 km²/yr), and in 18 subbasins it is low (mean: 0.94 hm/100 km²/yr). The spatiality of the SE at the third-order sub-basin is shown in Figure 7c. Among the 17 third-order sub-basins, the chance of soil erosion is high in 8 subbasins (mean: 2.38 hm/100 km²/yr), moderate in 2 (mean: 2.04 hm/100 km²/yr), and low in 7 subbasins (mean: 1.32 hm/100 km²/yr). High soil erosion at third-order subbasins are mostly located in the north-eastern portion whereas moderate soil erosion are found in the western portion and low erosion are found in the southern portion. The four fourth-order subbasins exhibit an average erosion rate of 1.70 ranging from 0.76 to 2.77 (Figure 7d). In the SRB, the subbasin with the highest rate of SE is found in the western section, while the subbasin with the lowest rate is found in the eastern part.

Figure 7. Sub-basin order-wise estimation of the SPR of the Sankha river. (a) SPR of first order basin, (b) SPR of second order basin, (c) SPR of third order basin, and (d) SPR of fourth order basin.

3.4. Factors controlling soil erosion

The SE potential of the SRB was calculated using the RUSLE model with six commonly used parameters. The catchment area of the basin exhibits high SE with a low R factor, and vice-versa. This signifies that the other factors are strongly influencing the SE potential. Regarding the soil types of SEB the entire basin is dominated by two soil types: lithosol (I) and dystic nitosol (ND). The lithosol (sand: 58.9%; silt:16.2%; clay:24.9%; OC: 0.97% at topsoil sand: 65%; silt: 17%; clay:27%; OC: 0.4% at sub-soil) covers 85% area of the SRB. The soil type ND (sand: 38.9%; silt:17.6%; clay:43.6%; OC: 1.57% at top-soil and sand: 31.9%; silt:13.8%; clay:54.4%; OC: 0.44% at sub-soil) covers 14% of the SRB and highly concentrated in the southern part. Besides, LF is observed in a very negligible portion of the basin. The presence of FL on hills and hill slopes is a common occurrence in SRB. The soil type LF (sand:74.6%; silt:9.6%; clay:15.6%; OC:0.39% at top-soil and sand:67.7%; silt:8.9%; clay:23.4%; OC:0.25% at sub-soil) covers 1% of the SRB.

Lithosol having a shallow nature and limited soil development, often has poor structural stability and is prone to erosion, especially if the rock fragments provide pathways for water to move quickly over the surface. On moderately sloping stable land formations, ND develops. This soil containing high clay content indicates high soil erodibility and resists erosion. While it is inherently highly erodible, the erodibility can still be influenced by factors such as slope, vegetation cover, and land management practices. The FL is characterized by clay illuviation, indicating the movement of clay from upper horizons to lower horizons. The presence of clay can enhance the soil’s ability to retain water and nutrients. However, the erodibility of can still be influenced by factors such as slope steepness and land management practices.

Additionally, the run-off characteristics also support the erosion variability in the magnitude of SE among the sub-basins. For example, among the subbasin, the first-order basin has recorded a low average runoff rate (1.49 cm/km²) compared to the other sub-basins (1.50 cm/km² for second-order, 1.53 cm/km² for third order and 3.91 cm/km² for fourth order basin). These opposite findings indicate that dynamic land use and land cover (especially the expansion of settlement and fallow land and reduction of waterbody area) increase the runoff rate causing significant SE. The majority of the basin’s agricultural land is situated on gentle to moderate slopes. Poor farming practices have an impact on the soil’s surface and structure, especially in agrarian areas with mild to moderate slopes. The influence of the C-factor on SE seems to be of significant relevance due to the greater concentration of agriculture in the basin.

The present investigation has also focused on the relationship between the determinants (independent variables) and SY (dependent variable) and the identification of leading factors at every basin order (first to fourth) using the PLSR model (Figures 8a–d and 9a–d). Soil erosion potential was first thought to be predicted by 26 independent variables under 5 major categories (8 under basin geometry, 5 under drainage, 4 under topography and climate, and 5 under LULC and 4 under soil). Finally, influencing variables have been retrieved for each dependent variables depending on the proportion of variance explained. Cross-validation analysis, loading PLS factors, XY scores, X weights, and variable importance for the projection (VIP) are among the PLS outcomes. Therefore, this section explains the PLS model’s empirical findings in each basin order. Four elements in the PLSR model of the first order basin have an explained variance (jointly by F1 and F2) greater than 70% (A: 82%, P: 94%, RR: 89%, DI: 89%). Furthermore, we discovered that six elements: P (1.02), CR (1.38), CR (1.39), RR (2.54), DI (2.46), and vegetation (1.047) have VIP greater than 1, implying an influential impact on SY (Table 5). Besides, only the variables with weights greater than 0.5 are taken into account when interpreting X weights.

Figure 8. PLSR analysis for identifying the leading factors of SY. (a) Variable importance for the projection (VIP) of SY for the first-order basin, (b) PLSR loading plots for dimension reductions of the factors responsible for SE yield for the first-order basin, (c) Variable importance for the projection (VIP) of SY for second-order basin and (d) PLSR loading plots for dimension reductions of the factors responsible for SE yield for second-order basin.

Figure 9. PLSR analysis for identifying the leading factors of SY. (a) Variable importance for the projection (VIP) of SY for third-order basin, (b) PLSR loading plots for dimension reductions of the factors responsible for SE yield for third-order basin, (c) Variable importance for the projection (VIP) of SY for fourth-order basin and (b) PLSR loading plots for dimension reductions of the factors responsible for SE yield for fourth-order basin.

Table 5. Order-wise factor loading score of different parameters in the PSLR model.

Broad categories	Basin orders/ independent variables	First order		Second order		Third order		Fourth order
		Factor (1)	Factor (2)	Factor (1)	Factor (2)	Factor 1	Factor 2	Factor 1	Factor 2
Basin geometry	A	0.096	0.301	0.094	0.224	0.146	0.066	0.176	0.503
	L	0.148	0.261	0.032	0.122	0.067	0.110	0.525	0.229
	P	0.175	0.331	0.034	0.112	0.083	0.054	0.159	0.506
	Cr	−0.522	−0.504	0.135	0.212	0.215	0.034	−0.110	−0.351
	EL	−0.009	0.066	0.113	0.128	0.228	−0.194	−0.196	0.247
	Ff	0.017	0.083	0.112	0.117	0.505	−0.176	−0.193	0.259
	Cc	0.540	0.507	−0.097	−0.254	−0.539	0.320	0.100	0.298
	ELL	0.043	−0.007	−0.110	−0.162	−0.598	0.271	0.205	−0.208
Drainage	DF	−0.108	−0.266	0.052	−0.210	−0.026	−0.327	−0.205	0.103
	DD	−0.062	−0.177	−0.099	−0.152	−0.112	−0.255	−0.206	0.122
	CCM	0.034	0.225	0.142	0.202	0.198	0.177	0.516	−0.041
	LOF	0.034	0.225	0.142	0.202	0.198	0.177	0.516	−0.041
	DT	−0.094	−0.220	0.040	−0.175	0.000	−0.350	−0.197	0.166
Topography and climate	RR	0.585	−0.273	−0.554	0.333	−0.574	0.327	−0.191	0.020
	DI	0.570	−0.216	−0.531	0.331	−0.557	0.293	−0.196	0.009
	HI	0.010	−0.041	0.056	−0.209	0.050	−0.254	−0.148	0.214
	Rain	−0.015	−0.016	0.046	−0.030	0.185	0.230	0.222	0.048
LULC	Cul	−0.072	0.278	0.177	0.187	0.256	0.060	0.227	−0.022
	Fallow	−0.020	0.296	0.157	0.219	0.205	0.097	0.237	0.047
	Sett	−0.019	0.319	0.209	0.102	0.168	0.070	0.228	0.102
	Veg	0.343	0.081	−0.275	0.103	−0.125	−0.006	−0.157	0.336
	Wat	−0.069	0.035	0.093	−0.059	0.172	−0.012	0.214	0.060
Soil	Sand	0.029	0.038	−0.151	−0.226	−0.161	−0.227	−0.200	−0.024
	Silt	−0.029	−0.038	0.151	0.226	0.161	0.227	0.200	0.024
	Clay	−0.029	−0.038	0.151	0.226	0.161	0.227	0.200	0.024
	OM	−0.029	−0.038	0.151	0.226	0.161	0.227	0.200	0.024

A = basin area, L = basin length, P = basin perimeter, EL = elongation ratio, Cr = circularity ratio, Ff = form factor, Cc = compactness coefficient, ELL = ellipticity index, DF = Drainage frequency, DD = drainage density, CCM = constant of channel maintenance, LOF = length of overland flow, DT = drainage texture, RR = relative relief, DI = dissection index, HI = hypsometric integral, Rain = Rainfall, Cul = cultivation land, Fallow = fallow land, Sett = settlement, Veg = vegetation, Wat = waterbodies, OM = soil organic matter.

The following are the PLS factor 1 weights for the dependent variable: RR (0.585), DI (0.570), CC, and (0.54) (Table 5). This implies that topography plays a significant role in SE. Additionally, the preservation of CR (0.504) and CC (05.07) in PLS factor 2 explains the significance of the basin geometry factor. Likewise, for second order basin, explained variances of more than 70% are retained by RR: 92% and, DI: 94%; VIP of more than 1 is recorded for RR (2.46), DI (2.35) and vegetation (1.27); and X significant weights are RR: −0.55, DI: −0.53 at F1. Therefore, it can be argued that relative relief and dissection index jointly have a significant influence on erosion at the second-order basin scale the correlation coefficient values also confirm the importance of the RR and DI and CC in SY of SRB (Figure 10). Moreover, for third order basin, explained variances of more than 70% are retained by DF: 86%, DD: 79%, and CCM: 73% LOF: 73%, DT: 87% DI: 71% Cul: 74% and, Fallow: 76%; VIP greater than 1 are recorded for CR (1.02), EL (1.09), FF (1.12), ELL (1.41), DI (1.47) and Cul (1.06); and X significant weights are FF: 0.50, CC: −0.53, ELL: −0.59, RR: −0.57, and DI: −0.55 at F1.

Figure 10. Heatmap explaining the correlations among the dependent and independent variables of SE (L – basin length; A – basin area; P – basin perimeter; Cr – circularity ratio; RR – relative relief; DD – drainage density; DF – drainage frequency; Er – elongation ratio; Ff – form factor; CC – compact coefficient; Di – dissection index; HI – Hypsometric integral; SE – soil erosion; SY – sediment yield; SDR – sediment delivery ratio).

Therefore, it can be argued that basin shape, drainage characteristics, and land use (cultivation pattern) are the key factors contributing to the SE at the third order basin scale. Finally, for fourth order basin, explained variances of more than 70% are retained by L: 80%, FF: 76%, ELL: 78% ELL: 83%, DF: 79% DD: 81% CCM: 84% LOF: 83%, DT: 76, Cul: 93% and, Fallow: 98%, Sett: 88, sand: 70, silt: 70, clay: 70; VIP more than 1 are recorded for L (1.14), EL (1), ELL (1), DF (1), DD (1), EL (1), CCM (1), LOF (1), DT (1), and Cul (1.06); and X significant weights are FF: 0.50, CC: −0.53, ELL: −0.59, RR: −0.57, and DI: −0.55 at F1. and RR: −0.55%, DI: −0.53 at F2. The PLSR approach investigated how the number of components increases as the basin order increases. More importantly, manmade variables such as various land uses have evolved into influencing factors for SE in the higher basin order. In particular, the land use-specific analysis shows that compared to other LULCs, agricultural and fallow land contribute more to SE.

3.5. Implications of the study for land management

The findings of the current study will help the stakeholders for erosion control in the SRB. Selecting severe erosion-prone locations is a pre-requisite strategy for the efficient allocation of scarce resources because conserving all the basin areas suffering from soil loss of different magnitudes is neither economically feasible nor technically possible (Zhou et al. 2019). To enhance the health of the landscape ecology by boosting ecosystem diversity and productivity, many studies across the world have estimated the rate of SE at different spatio-temporal scales and focused on the regions that are most vulnerable to SE. For example, Elnashar et al. (2021) proposed an urgent policy impact for the areas of the Blue Nile Basin where SE is occurring at a rate of more than 10 t ha⁻¹yr⁻¹. Furthermore, it has been demonstrated that countries in Asia, Europe, Africa, and America focus their policy recommendations on areas with erosion rates higher than 20 t ha⁻¹yr⁻¹ (Aneseyee et al. 2020; Majhi et al. 2021). On the other side, it has been observed that in tropical India, where SE rates exceed 20 t ha⁻¹yr⁻¹ in large river basins, attention is being paid to SE mitigation through civil and bioengineering approaches. Small river basins with SE rates of greater than 20 t ha⁻¹yr⁻¹, however, are planned as priority sites for SE control. SE is posing challenges to the small SRB in the tropical area. Therefore, in the areas of SRB where soil loss rates are more than 10 t ha⁻¹yr⁻¹, we recommend urgent, focused implementation. These target regions make up 27% of the research area and represent 96% of the soil loss rate in the SRB.

To reduce soil loss in such target areas, appropriate land management techniques should be put into place. As the northern portion of the basin experiences high soil erosion due to steep slopes physical measures such as soil bunds, stone bunds and percolation ditches may be implemented. Besides, agronomic practices (changes in cropping and land management practices such as contour ploughing, and terracing can adopted to minimize soil erosion. Moreover, improved agronomic practices are recommended for the areas located in the middle and western parts of the SRB that have low to moderate SE so that future erosion may be checked. Therefore, to successfully address erosion issues, conservation planning might be improved in the SRB using the results of the current study. This could create chances for collaboration between the local government and the local population.

4. Conclusions

The basin-scale SE study of the Sankah River reveals that most of the river basin has low SE, and the upper portion of the SRB is more susceptible to erosion. With a range of 10 t ha⁻¹yr⁻¹ to 50 t ha⁻¹yr⁻¹, the mean SE is 20 t ha⁻¹yr⁻¹. In particulars, low erosion potential (SE < 10 t ha⁻¹yr⁻¹) covers 89% of the SRB, whereas ow erosion potential (SE > 40 t ha⁻¹yr⁻¹) only accounts for 3% of the SRB. The forest land has experienced high soil loss due to steep slopes, accounting for 58.39% of the total. Besides, the results of sediment delivery ratio and sediment yield have also been integrated with SE. The average SDR for first, second, third and fourth order basins are 0.60, 0.63, 0.42, and 0.37 respectively while the average SY of the basin varies from 102.37 to 332.4 t ha⁻¹yr⁻¹. The PSLR models reveal that circularity ratio, compactness coefficient, relative relief, and dissection index are the main factors controlling soil erosion for the first-order river basin. Besides, relative relief, dissection index, and vegetation cover for the second-order River basin. circularity ratio, ellipticity index, form factor, compactness coefficient, basin elongation, relative relief, dissection index, cultivation for third order and basin length, basin elongation, drainage frequency, drainage density, constant of channel maintenance fallow land, cultivation land, settlement for fourth order are the significant factors for promoting the soil erosion.

Therefore, the explored method shows that increasing the basin order increases the number of factors. More importantly, anthropogenic factors such as different land uses have become influential factors for SE in the higher basin order. Precisely, the land use-specific study indicates that fallow and agricultural land contribute more to SE than other LULCs. Besides, the slope-specific study indicates that the steeper slope, often found in the upper (north) part of the basin, is more prone to SE. Therefore, the present study clearly demonstrates that different basin orders require different management techniques to prevent SE. Especially, the lower-order basin demands slope stability through aforestation, systematic agricultural practice, controlling the headward erosion of gullies, etc. Overall, the higher-order basins demand controlled anthropogenic activities and sustainable restoration of degraded lands.

The present effort involves research to characterize soil erosion in light of the present-day situation. This approach aims to depict spatiotemporal variations more accurately, fostering the integration of basin-scale erosion dynamics with land management strategies. Furthermore, future land management would better be framed if social perspectives especially social needs and demand could be integrated. Another dimension emanating from the present investigation is the application of sophisticated geospatial and geostatistical models to simulate the future trend of erosion dynamics in different geo-climatic, physiographic environments.

Disclosure statement

There is no conflict of interest.

Data availability statement

Data available on request from the authors.