This study was conducted in the mid- and high-mountain regions of the Tamor River basin (TRB), one of the major sub-basins of the Koshi River (Figure 1). TRB originates from the eastern Himalayan regions of the Kanchanjunga range and joins the main Koshi River at Tribeni, Nepal. The study area is situated between 87.08°-88.06°E and 26.51°-27.42°N, which drains an area of nearly 4189 km². Its elevation ranges from 93 to 4647 m with a mean elevation of 1814 m. The mean annual rainfall in the study area is approximately 1587 mm yr^-1. Northeasterly and southwesterly wind systems that occur in the winter and summer, respectively, are the major sources of rainfall (Dahal et al., 2020). Summer (June-September) is the major rainy season, providing nearly 80% of the annual rainfall (Dahal et al., 2020). According to the FAO/UNESCO Digital Soil Map of the World, the dominant soil types in the study area are dystric cambisols (66.8%), lithosols (23.6%), and humic acrisols (9.7%) (FAO/UNESCO, 2007). Forest land, agricultural land, grassland, shrubland, barren land, water, and built-up areas are the major LULC in the study area (Uddin et al., 2016). Forest land covers the largest area, followed by agricultural land, while built-up areas constitute the smallest share of the basin (Sinha et al., 2019). Forest, shrub, and grasslands predominantly dominate the steep and high elevation areas of the basin, while agricultural lands and built-up areas are concentrated in relatively less steep slopes and moderate elevation areas along the side/banks of rivers. As part of the Himalayan region, the TRB is one of the key SL and SE hotspot sites in the world because of its inherent unconsolidated and fragile geology, steep slopes, torrential rainfall, and intensive anthropogenic activities (Jain et al., 2001; Golosov and Walling, 2019). This basin covers approximately 11% of the upper Koshi River basin and contributes around 16% (16 Mt yr^-1) of the total sediment load to the basin (Sinha et al., 2019). The basin is generally the most sensitive and fragile environment. A small disturbance to its natural components significantly affect its overall sediment dynamics, making this basin the preferred choice for the current study.

For this study, rainfall data for a period of 30 years (1987-2017) at 16 gauge stations were obtained from the Department of Hydrology and Meteorology (DHM), Nepal. Similarly, monthly scale discharge and sediment concentration data from 2006 to 2013 for the validation of our SE simulation results were obtained from DHM, Nepal. In addition, a 12.5 m spatial resolution digital elevation model (DEM) and a digital soil map (scale 1:5,000,000) were used to generate the topographic and soil erodibility (K) factor maps of the study area, respectively. The LULC maps of the study area for the years 1990, 2002, 2010, and 2021 were generated using Landsat 5/7/8 images from the United States Geological Survey (USGS) archive available in the GEE cloud platform. Detailed data types and the sources used in this study are listed in Table S1. The overall study process from data acquisition to the final SL and SE analysis is shown in Figure 2.

In this study, the InVEST model was chosen for SL and SE estimation because it is more user-friendly and viable for data-scarce basins than other models. It also enables the evaluation of sediment export, which is often considered one of the major limitations of the Revised Universal Soil Loss Equation (RUSLE) model. The InVEST model calculates the SL and SE for a specific pixel based on the concepts of the RUSLE (Renard et al., 1997) and hydrological connectivity index (Borselli et al., 2008) models. The model requires input datasets of geospatial layers consisting of the DEM, LULC maps, K factor, and rainfall erosivity (R) factor maps in a raster file format, a study area polygon, and a look-up table containing cover management (C) and conservation practice (P) factors in a CSV file format (Figure S1). The parameters and equations used to generate input datasets for this model are discussed below.

A long-term average of the measured sediment data should be imperative to calibrate and validate the InVEST model (Sharp et al., 2018). The model results of this study were calibrated using average observed sediment load data measured at the final outlet of the study area (at Mulghat station) from 2006 to 2013. As there are no plot-based experimental studies and long-term observed sediment data in the study area, the simulation results were validated by comparison to the results of other studies conducted in KRB using the radionuclide techniques and RUSLE model. In this study, we adopted the approach proposed by Quilbe et al. (2006) to manage the missing sediment load data. Based on the observed sediment load data together with the results of previously published studies in KRB (Sinha et al., 2019; Yigez et al., 2021; Yuan et al., 2021), it was possible to confirm that the model and parameters are suitable for simulating the spatiotemporal variation of SL and SE in this case study area. The deviation between the simulated data in this study and the reported SL and SE in the previous studies varies from 1.67 to 46.84 t ha^-1 yr^-1 and 0.06 to 14.32 t ha^-1 yr^-1, respectively (Table 1).

To evaluate the spatiotemporal association between the change in SL and SE with LULCC, we calculated the mean SL and SE change rates in each LULC type and compared them with the corresponding LULC dynamics from 1990 to 2021 in KRB. An average of the annual rainfall data from the year 1987 to 2017 was used as model input to isolate the impact of the temporal variation in rainfall on SL and SE rates and independently examine the impacts of LULCC on SL and SE dynamics (Gashaw et al., 2018). In this study, we used a correlation analysis to present the association between simulated SL, SE, and the proportion of different LULC classes in the study area. The zonal statistics and fishnet analysis tools in the ArcGIS 10.2 environment were used to generate the statistical samples. Finally, IBM SPSS Statistics 20 (IBM Corp., Armonk, NY, USA) was employed for the correlation analysis in this study. In this study, the normality of the data was initially tested, and the variables that failed to meet the normality criteria were logarithmically transformed.

To examine the severity of SL in the study area, the SL rate was further categorized into six classes, as proposed by Singh et al. (1992) (Table 2). To evaluate the variability of SL along the slope gradients in this study, we categorized the slope gradient into six steepness classes (i.e., 0°-5°, 5°-8°, 8°-15°, 15°-25°, 25°-35°, >35°) (Wang et al., 2016; Teng et al., 2019). The variation in SL rate analysis among these categories was computed using the zonal statistics tools in the ArcGIS environment.

The LULC maps of the study area in the years 1990, 2002, 2010, and 2021 are depicted in Figure 3. The overall accuracy and Kappa coefficient values of each of these maps range from 90 to 93% and from 0.86 and 0.91, respectively (Tables 3 and S3). This indicated that the classified maps had a very good agreement with the reference data (Monserud and Leemans, 1992).

The proportion and spatiotemporal dynamics of each LULC type at the four time frames are shown in Table 4. The change detection result depicted that the LULC in the study area has undergone a considerable change during the last three decades (1990-2021), since nearly 21% of the study area has experienced the LULC transition (Figure 3; Table S4). The largest area loss was detected for agricultural land. The area of forest land, shrubland, grassland, and built-up areas showed a continuous increase during the study period at the expense of agricultural and barren lands. From 1990 to 2021, the amount of agricultural land decreased from 28% to 18% (-432 km²), which is an average loss of approximately 14 km² yr^-1. Conversely, the area of forest land, grassland, and shrubland increased by ~4% (173 km²), 2.4% (154 km²), and 3.7% (99 km²), respectively (Table 4). Out of the total agricultural land loss, 69%, 19%, and 10% were converted to forest land, grassland, and shrubland, respectively. The largest proportion of change was observed for built-up areas (> six-fold from the initial state) and the smallest change was identified for forest lands (Table 4).

Compared to the three study periods (time intervals) (1990-2002, 2002-2010, and 2010-2021), the smallest and largest changes in area proportion were detected in the second and third study periods, respectively (Table 4). During the third time interval, the extent of grassland (+1.87%), shrubland (+1.5%), and built-up (+0.12%) areas substantially expanded over the other two periods (Table 4). Meanwhile, agricultural land experienced a drastic loss (-5.46%) during this time frame. The major hotspot areas of LULCC were mostly observed in sloping areas and areas adjacent to natural forests and far from settlement centers and river networks (Figures 3a-3d).

The rate of SL and SE from 1990 to 2021 at a basin and sub-watershed (SW) scale are shown in Figures 4 and 5 and Tables 3 and S5. Overall, the mean rates of SL and SE in the study area steadily decreased over the study period. It can be seen that the mean SL and SE rates in the study basin dropped by approximately 48% and 60%, respectively, between 1990 and 2021. The mean SL decreased from 45.21 t ha^-1 yr^-1 in 1990 to 32.18, 30.31, and 23.33 t ha^-1 yr^-1 in 2002, 2010, and 2021, respectively. Similarly, the corresponding mean rate of SE gradually decreased from 7.49 t ha^-1 yr^-1 in 1990 to 4.78, 4.20, and 2.98 t ha^-1 yr^-1 in 2002, 2010, and 2021, respectively. As a result, the total SL in the study area dropped from 18.94×10⁶ t in the base year of 1990 to 13.48×10⁶, 12.69×10⁶, and 97.73×10⁶ t in 2002, 2010, and 2021, respectively, which is equivalent to an average decline of 2.9×10⁵ t yr^-1. Likewise, the corresponding total SE decreased from 3.14× 10⁶ in 1990 to 2.0× 10⁶, 1.76 × 10⁶, and 1.25 ×10⁶ t yr^-1 in 2002, 2010, and 2021, respectively.

At the SW scale, the severe SL and SE hotspots were concentrated in the middle and mid-western parts of the river basin. The maximum mean SL rates were all observed in SW5 in the years 1990, 2002, 2010, and 2021 and their SL rates were 88.66, 62.13, 55.03, and 52.72 (t ha^-1 yr^-1), respectively. The corresponding maximum SE rates in 1990, 2002, 2010, and 2021 were observed in SW5, SW5, SW14, and SW5, respectively, while the minimum rate was observed in SW7, SW2, SW7, and SW18, respectively (Table S5). These maximum SL and SE hotspot SWs were all located in the middle part of the study area, as they were largely dominated by sloping farmlands (~40%, 58%, 30%, and 28%, respectively) and less vegetation coverage (Figure 2).

The study also presented a temporal variation in SL and SE rates among the severity classes (Table 1 and Figure 6). During the study period, the proportion of slight, moderate, and high SL severity classes steadily increased (by 8.22%, 1.28%, and 0.56%, respectively) (Figures 6a and 6c, Table 6) The remaining three classes (very high, severe, and very severe) gradually decreased by 0.43%, 1.96%, and 7.67% for SL and 2.98%, 1.98%, and 1.09% for SE from 1990 to 2021 (Figure 6, Table 4 and S6). More specifically, the two dominant classes, slight and very severe, experienced the largest areal variation (an increase of 8.22% and a decrease of 7.67%, respectively) (Figure 6a and Table 6). It can be seen from Figure 6 that approximately 31% of the study area in the base year of 1990 was found under the category of high to very severe SL rate. However, the areal proportion under this category dropped to 21% in 2021. These intense SL severity classes were predominantly concentrated in the middle and mid-western parts of the study basin comprising SW5, SW8, SW12, SW13, SW14, SW16, and SW17 (Table S5 and Figures 4a-4d). Despite these classes covering a small percentage of the basin area, their contribution to the total basin SL was very large (> 88%). Conversely, the north and northeastern parts of the study basin (SW-13, SW6-7, and SW9-10) had relatively low rates of severity (Table S5). Since SE is strongly determined by the SL rate, the spatial dynamics of SE in this basin generally showed a similar spatial distribution pattern to the SL rate (Figures 4 and 5a-5d).

During the last three decades (1990-2021), remarkable changes in LULC were experienced in the study area (Figures 3a-3d). Along with these changes, considerable dynamics in the rates of SL and SE were also exhibited in this basin (Figures 4 and 5a-5d). Table 5 summarizes the rate of SL and SE from each LULC type as a function of time. We found that the mean SL and SE rates varied with LULC types in decreasing order from agriculture > barren land > built-up > grassland > shrubland > forest land over the study periods. The change in SL and SE rates in the study area was principally linked with the dynamics in the proportion of the two dominant LULC (agricultural and forest lands) (Figure 3 and Table 5). The SL and SE rates detected in agricultural land were much higher than that of all other LULC, ranging from 95 to 130.77 and 12.65 to 22.25 (t ha^-1 yr^-1), respectively (Table 5). This rate is more than six times larger than the maximum SL tolerance rate of the Himalayan region (15 t ha^-1 yr^-1) and four times larger than that of the mean SL rate estimated in the current study (22.22 t ha^-1 yr^-1) for the study area. On the contrary, the minimum SL and SE rates were observed in forest land, varying from 6.36 to 10.24 t ha^-1 yr^-1 and 0.57 to 0.75 t ha^-1 yr^-1, respectively.

The correlation analysis also showed a clear linkage between LULCC with SL and SE dynamics in the study area (Table 7). It can be seen that SL and SE showed a significant (p<0.05) negative correlation with the proportion of forest land with R² values ranging from -0.45 to -0.55, and -0.47 to -0.57, respectively, from 1990 to 2010. On the other hand, the SL and SE rates had a significant (p<0.01) positive correlation with the proportion of agricultural land change with R² values ranging from 0.64 to 0.71 and 0.63 to 0.77, respectively (Table 7). Moreover, the relationship between grassland with SL and SE was significant (p<0.05) and moderately correlated in the year 2010. Generally, this analysis suggested that SWs dominated with a high proportion of vegetated LU types were less subjected to SL and SE, and vice versa.

The spatial distribution of SL and SE for the years 1990, 2002, 2010, and 2021 across topographic factors was evaluated by overlaying the spatial distribution map of SL and SE with the slope aspect and slope gradient maps of the study area (Table 8; Figure7a and 7b). The mean SL and SE in different slope gradients generally showed an increasing trend with increasing slope gradients up to 25°, and then gradually decreased with increasing slope gradients beyond this gradient in the study area (Table 8). The mean SE rate showed a decreasing trend from 1990 to 2021 across all slope classes. Concurrently, the mean SL rate also showed a decreasing trend from 1990 to 2021 under slope gradient classes beyond 25°. In the current study, SL and SE also varied with slope aspects (Figures 7a and b). Sun-facing aspects, which receive better solar radiation (SE, S, E, and SW), showed higher SL and SE rates than shady slopes in the study basin. On the contrary, shady slopes comprising flat, NW, and N aspect slopes exhibited the lowest rates of SL and SE in the study area. We also found that the SL and SE rates under sun-facing and semi-sunny slopes experienced a decreasing trend from the base year of 1990 to the final year of 2021. A slope gradient from 8°-35° and sunny and semi-sunny slope aspects including S, SE, E, and SW accounted for nearly 91% and 57%, respectively, of the total SL and SE in the final year of this study.

In environmentally fragile regions, LULCC has a significant impact on sediment flux (Bakker et al., 2008; Xu et al., 2011). This study investigated the responses of SL and SE to LULCC in the study area of Nepal for a time window of 32 years (1990-2021). The results of the current study showed a substantial decrease in the average rates of SL (48%) and SE (60%) from 1990 to 2021. Concomitantly, the corresponding total SL and SE decreased by 91.66×10⁵ and 18.89×10⁵ t, respectively. The study suggested that there has been a considerable shift in SL severity classes (about 10% of the high to very severe classes transforming to moderate to slight classes) over the study period in the study area. The significant drop in SL and SE rates was attributed to the continuous conversion of intensified (agricultural land) to de-intensified (forest, grass, and shrublands) LULC types associated with the implemented government and community-based forest development, protection, and other soil and water conservation interventions in the late 1980s and early 1990s (such as the Master Plan for the Forestry Sector in 1988 and the Forest Protection Law in 1993) in the study area (Gautam et al., 2004; Li et al., 2017; Chalise et al., 2019). According to Gautam et al. (2004), the formalization of local forestry institutions by forest user groups after the implementation of the community forestry program also contributed positively to the increase in vegetation cover. In addition to the above factors, migration-triggered farmland abandonment and associated stimulation in LULCC also had positive implications on the decline in SL and SE in the study area. During the study period, other LULC classes also showed both gradual declining (in barren land) and increasing (in built-up and water area) trends, but they had no significant effect on the SL and SE rates because of their small area.

Although the mean rate of SL in this study showed a continuous decreasing trend from 1990 to 2021, the mean SL rates of all four study years (about 21% of the study area in 2021) are already beyond the maximum SL tolerance of 15 t ha^-1 yr^-1 in the Himalayan region (Jasrotia and Singh, 2006) and almost double the global tolerance rate of 10 Mg ha^-1 yr^-1 (Borrelli et al., 2017). The estimated rate of SL in the study area is in line with other sub-basin and watershed-scale study results; however, it is higher (nearly 89%) than previous results reported by Uddin et al. (2016) in the entire KRB, as this study utilized coarse resolution DEM (90 m) and incorporated the Terai region of KRB (an area of large flat topography), which led to the low rates of SL and SE. Moreover, the different models used for estimation might be the other possible reason for the variation.

The type of land use and management practice applied to a specific watershed strongly determines the variability in vegetation cover and, thereby, the process of SL and SE (Zhou et al., 2019). LULC with higher vegetation cover has higher soil conservation services (Bakker et al., 2008). For example, a study by Zuo et al. (2016) found that reforestation decreased the rate of SE by more than 40% in the Huangfuchuan watershed, China. Similarly, vegetation coverage greater than 60% has the potential to reduce annual SE by up to 60% (Chen et al., 2021). In the current study, the strongest soil conservation function was observed in forest land use followed by shrubland and grassland, since these LULC categories can provide better ground cover against erosive agents than agricultural and barren lands.

The intensity of SL and SE was considerably higher in the middle and mid-western parts of the study area in all four study years, since these areas are dominated by a high proportion of sloping farmlands and other LULC with low vegetation coverage. Conversely, a lower intensity of SL and SE was observed in the north and northeastern parts of the study area. In addition to agricultural land proportion, the inherent lower inter-annual ground/vegetation cover variability of forest, grass, and shrubland might also be another possible reason for the relatively lower rates of SL and SE in the north and northeastern part of the study area, despite the topographic and K factor potential (steep slope and high elevation) of these regions being very prone to soil erosion.

In this study, the variability of SL, SE, and LULCC among the three study periods was also evaluated. Accordingly, the period from 1990 to 2002 was highly characterized by the largest drop in SL and SE over the other two periods. During this period, the extent of agricultural land drastically decreased (a drop of nearly 33.32%) at the expense of forest, grassland, and shrubland expansion. Concurrently, the total amount of SL and SE substantially decreased by approximately 32% and 39%, respectively (Table 5). Comparable results have also been reported in other studies in different regions (Bakker et al., 2008; Sun et al., 2013; Yan et al., 2018; Teng et al., 2019; Zhou et al., 2019). The LULCC and corresponding decrease in SL and SE during this period were mostly attributed to the implementation of different forest, soil, and water conservation interventions such as the community, leasehold, and private forestry programs introduced in the late 1980s and early 1990s. Reduced pressure on the natural landscape encourages natural regeneration, which, in turn, enhances the conversion of some shrublands into mature forest trees (Bakker et al., 2008; Sun et al., 2013; Li et al., 2017). In this regard, the application of improved agricultural technologies and practices that enable farmers to produce high yields (such as agroforestry practices, the promotion of improved seed/vegetative varieties, the expansion of irrigation systems, and pesticide and fertilizer application) in recent decades could also be another reason for the increase in vegetation cover, resulting in a decrease in SL and SE in the study basin (Gilani et al., 2017; Li et al., 2017). Comparatively, the magnitude of all three variables (SL, SE, and LULCC) in the period from 2002 to 2010 was lower than both the first and third study periods. The decrease in mean SL and SE was only 1.87% and 0.58%, respectively (Table 4). Differing from the other periods, the extent of shrubland in this period decreased by 0.28% compared with that of 2002. This could be due to the shift in human pressure from forest land to shrubland. The major sources of construction and fuelwood for rural communities in developing countries including Nepal are mostly forests. However, forest resources are now under better management conditions since the launch of the community forestry program (in the early 1990s) in Nepal, and thus, the pressure on forests for the abovementioned consumptions have drastically decreased. This is likely to have resulted in a shift in pressure from forest land to shrubland to compensate for the demand for wood obtained from community forest resources. This might lead to the selective clearing of shrubs and scant tree cover in this land use type, leading to areas with no shrubs that are dominated by grass cover. This could be one reason for the change in shrubland to grassland and the decrease in its coverage in the second period of this study. However, the intensity of LULCC from 2010 to 2021 was more substantial than both of the previous two periods. This might be due to not only the implementation of the national forest conservation measures (community forestry, afforestation, and reforestation), but also the other socioeconomic drivers such as the increasing rate of internal migration between 2010 and 2021 (Bhawana et al., 2017; Chidi, 2017; Paudel et al., 2020). In recent decades, there has been a high outflow of the active labor force from the fragile hill and mountain regions to plains and urban areas of the country in search of better infrastructure, fertile land, and access to off-farm income activities, etc. (Maharjan et al., 2020; Paudel et al., 2020). This high rate of outmigration creates a considerable shortage in terms of an active labor force for on-farm activities and, in turn, leads to an increased de-intensification and abandonment of agricultural land. The subsequent abandonment of farmland in marginal areas of the study basin has enhanced vegetation succession, replacing historically highly erosive cereal fields with grasses, dense shrubs, and forest communities, thereby reducing SL and SE. In addition to natural vegetation colonization, farmland abandonment in the study area also encourages farmers to grow trees such as uttis (Alnus nepalensis), chiraito (Swertia chirayita), loth salla (Taxus wallichiana), and cardamom (Elettaria cardamomumon) on these abandoned lands, thereby increasing vegetation cover and decreasing the rate of SL and SE in recent years (Pandey et al., 2016). Similar experiences of landscape transitions following depopulation/migration and a reduction in the rate of SL and SE have been reported in other parts of the world (Li et al., 2016; Bhawana et al., 2017; Kolecka et al., 2017; Xu et al., 2019).

Although the LULCC from 2010 to 2021 exhibited a relatively higher magnitude than the previous two periods, the changes in SL and SE rates from 2010 to 2021 were not as strong as in the first interval. This could be due to the concentration of the changing area in less SL- and SE-sensitive areas (steep slope, high rainfall, or areas with high erodible soils) compared to the previous two periods, since the severely eroded areas have already been changed to less intensified land use types (better vegetation coverage) during the first and second study periods. Moreover, the lower SL and SE rates exhibited during the third interval were attributed to the high change area proportion of non-vegetated LULC types (agricultural land) to less erosion-sensitive vegetated classes compared to the other two periods. Comparable results have also been reported in previous studies (Bakker et al., 2008; Zhu et al., 2018; Chen et al., 2021) stating that de-intensification (such as afforestation and abandonment) leads to a more favorable landscape pattern with respect to the reduction of SL and SE. On top of the above rationale, the plot size of the changed area also has a significant impact on the rate of SL and SE by influencing runoff and sediment connectivity (Wang et al., 2018; Fang, 2020). Due to the rugged nature of the study area, most of the agricultural plots that were abandoned during the third interval were mostly concentrated in high elevation and steep areas and they were small in size and fragmented. Thus, the response in SL and SE might not be as significant as the change in large-sized agricultural land plots.

In the study area, SL and SE varied with different slope aspects and gradients. Generally, the mean rates of SL and SE gradually increased with increasing slope gradients up to 25° and then decreased beyond this gradient. As can be seen from the spatial distribution LULC maps (Figure 3) in this study, the decreasing trend in the rate of SL and SE under slope gradients >25° might be due to the presence of high vegetation coverage in steep slope areas, since steep slope areas are mostly recommended for the use of forest, shrubland, or grassland in the study area. This result suggested that the increase in vegetation recovery/coverage can effectively weaken the impact of slope gradient on SL and SE rates. In addition to the observed variation in SL and SE across slope gradients, they also varied within the same slope gradient classes during the study periods, with a general decreasing trend from 1990 to 2021. This might be attributed to the improvement and expansion of vegetation coverage in steep and high-elevation areas associated with the recently implemented conservation measures. The slope aspect also showed a clear variation in mean SL and SE rates in this basin. High SL and SE were observed in S-, SE-, E-, and SW-facing slopes than in the north-facing counterparts. The slope aspect controls the solar angle and determines the rate of SL and SE through changing the local climate, and consequently impacts vegetation coverage and human activities (Chidi, 2017). Most crops need sufficient sunlight for proper growth, flowering, and fruiting, especially in areas where there is high rainfall; thus, sun-facing slopes, mostly south-facing in Nepal, are more suitable for most crop types than north-facing slopes (Chidi, 2017). This indicates that south-facing slopes are highly subjected to intensive human activities such as agricultural practices. This could be the reason for the high rate of SL and SE in sunny slopes in the study area. On the contrary, shady slopes are more suitable for non-farming practices. Chidi (2017) and Yan et al. (2018) found that the rate of farmland abandonment and vegetation recovery is higher on shady slopes than on sunny slopes. The result thus suggests that the ecological restoration and conservation measures should be continued and should focus more on sunny slopes, since this slope aspect is one of the hotspot areas of SL and SE in this basin.

Sediment sources in most mountainous basins such as the study area are very complex. Unmanaged road construction, gullies, landslides, and debris flow in mountainous basins contribute a significant amount of sediment load to the river channel. Thus, the sediment contribution from these sources should be accurately quantified using other modeling approaches in the future for better planning and management interventions. The temporal dynamics of structural soil and water conservation measures (P factor) play a major role in reducing SL and SE (Chen et al., 2021); however, they were not taken to account in this study due to reliable data limitations. Remote sensing data may not effectively detect existing tree crop areas, which have widely expanded in recent decades in the study site, such as tea plantations and extensive cardamom (Elettaria cardamomumon) and chiraito (Swertia chirayita) farming. Therefore, rigorous field verification and secondary data at the local level will be very helpful in addressing this problem. Although this study has the abovementioned limitations, they do not influence the general trend of the results and the conclusions drawn.

Text 1 Equations (EPEC) used to compute the soil erodibility factor (Sharpley and Williams, 1990) (Equations S1-S5):

where K is the soil erodibility factor; SAN, SIL, CLA, and C denote the percentage of sand, silt, clay, and organic carbon contents in the soil, respectively; SN1 = 1 - (SAN/100). The value of 0.1317 is the unit conversion factor; f_csand, f_cl-si, f_orgc, and fhisand denote the functions of highly coarse sand, clay, silt, organic carbon, and high sand contents in the soil, respectively.