Himalayan watersheds in Nepal record high soil erosion rates estimated using the RUSLE model and experimental erosion plots

The rising unpredictability in the food supply chain in many parts of the world is related to soil loss and poor agricultural output. The Revised Universal Soil Loss Equation (RUSLE), widely used for estimating soil loss, was applied in the western mid-hills in Nepal, with steep slopes and fragile geology. This region is at high risk for rapid soil erosion and mass wasting. To estimate soil loss, this study utilized the RUSLE model with experimental erosion plots in the Aadhikhola and Tinahukhola watersheds, capturing real-time erosion in the field. The annual soil loss for the Aadhikhola watershed is estimated at ∼41.4 tons ha−1 yr−1. In contrast, in the Tinahukhola watershed, soil loss is low (∼24.1 tons ha−1 yr−1). Although annual rainfall showed an increasing trend in both watersheds, the change in soil loss was statistically insignificant. The high erosion rates from the experimental plots in both watersheds support the model outputs. Results from the experimental plots recorded the rate of soil erosion for different land use as: irrigated agricultural land > rainfed agricultural land > forests. The trends highlight the role of human activities in enhancing soil erosion in these mountainous terrains in terms of medium to long-term perspectives. Therefore, sustainable agriculture practices in these terrains must investigate alternate ways to decrease soil erosion to support people's livelihoods.


Introduction
Precipitation is the most prominent agent fostering soil erosion along young and fragile mountain landscapes dominated by monsoons [1]. Such episodic monsoon showers affect soil productivity and agricultural output. For example, 9% of crop productivity has declined in Africa [2]. Similarly, in Europe, the erosion rate was ~2.46 tons ha − 1 , 1.6 times more than the rate of soil formation [3]. In Asia, ca. 663 million hectares of land are affected by soil erosion, which is ca. 20 tons ha − 1 yr − 1 . The main hotspots are China, India, and Indonesia. In the last 20-100 years, the rapid increase in human population and agriculture has worsened the situation in countries like China, Nepal, Lebanon, etc. [4].
Previous studies indicate that factors like rainfall and surface runoff, agricultural activities, vegetation cover, and slope result in differential soil loss patterns [5]. In addition, rapid growth in the human population and non-sustainable land-use practices enhance soil loss in the watersheds. Reports suggest that one-third of the world's total arable land has been lost due to soil erosion in the last 40 years, and it continues at a rate of 10 million hectares annually [6,7]. In particular, converting natural forests to agricultural croplands cause severe negative impacts leading to high soil erosion due to reduced vegetation cover [8,9]. These assessments predict that the world's topsoil could become unproductive within 60 years [9], which has grave consequences.
Nepal is facing unregulated and haphazard construction and land-use changes in the hills, making them more vulnerable to erosion and landslides; it destabilizes the fragile geology and alterations in local hydrological conditions [10]. Moreover, intense rainfall, mainly from June-September, has made the steep slopes vulnerable to high soil erosion and collapse [11,12]. According to the Ministry of Forests and Environment report [12], the average rainfall is projected to increase from 2.1 to 4.5% by 2045 and increase to 10.7-23% by 2100. The increasing precipitation trend, coupled with conventional tilling practice in mountainous watersheds, can quickly sweep away topsoil from the steep slopes [13]. This makes the fragile slopes of the Nepal mid-hills highly vulnerable to soil erosion and substantial loss.
Agricultural land in Nepal constitutes 28% of land cover of which 21% is cultivated, and 7% remains uncultivated [14]. Land degradation, mainly soil erosion affects agriculture and is a major challenge in Nepal [13]. Conventional tillage practice, observed in arable lands in Nepal, leads to significant loss of topsoil during the rainy season due to steep and fragile landscapes situated on sloping hills [15]. Consequently, the amount of soil loss from conventional tillage is 5.5 Mg ha − 1 (1 Mg = 1000 kg) more than the reduced tillage system, and it is accompanied by a decline in nutrients due to enhanced soil loss [15]. Previous studies indicate that the degradation of arable land in Nepal includes the structural deterioration in soil caused by the increase in bulk density and runoff, and a decrease in porosity [13]. The use of tractors and harvesters for tillage results in soil compaction. Similarly, pastures are used for free grazing of animals, resulting in enhanced soil compaction and removal of vegetation and leaf litter. As a result, this soil is more susceptible to mass wasting and soil loss. In 2006, it was reported that 0.2 million ha of land was affected by soil degradation, accounting for 1.4% of arable land in Nepal [13]. Poor land management led to a decline in soil quality and crop production, and realistic plans are urgently needed to enhance local participation to decrease land degradation [14].
The Universal Soil Loss Equation (USLE), developed by Wischemeir and Smith in 1965 and further improved in 1978, is the most widely used empirical model [16]. RUSLE accounts for erosion due to rainfall and comprises five main factors, which are rainfall erosivity (R), soil erodibility (K), slope length (L), steepness (S), cover factor (C), and support practice factor (P) [16]. This model has been frequently used to estimate soil loss in small watersheds, including those in Nepal [17]. Although empirical models like RUSLE are simple and helpful in assessing soil erosion, these models require additional practical/on-site verification processes for broader applicability and better predictability [18]. One of the significant limitations of this model is the precise estimation of deep-seated gully erosion. As a result, erosion rates predicted by RUSLE have a low threshold in watersheds exposed to gully and large-scale mass movement, e.g., the active occurrence of landslides. Hence, runoff plots for erosion have been used widely for a better assessment of topsoil loss and management practices and to verify model outputs [19,20]. Several parameters are weighed carefully while constructing these plots, e.g., precipitation, amount of soil transported, surface evaporation, runoff, infiltration, moisture, fertility, and agricultural output, providing a robust assessment [19]. Most existing studies on soil erosion consider the average precipitation values derived from sensors for modeling. The main drawback is that sensors lack the flexibility to mimic the dynamism of natural rainfall [21]. Similarly, laboratory studies cannot fully simulate the complicated interaction of precipitation with surface features in any watershed.
The mid-hills in Nepal are prone to extensive soil erosion [22,23]. However, regional studies on soil erosion have used models without any field plot methods or corroboration of the proposed inferences [22] based on models. Utilizing both schemes involving soil erosion models and erosion plots for estimating soil loss at the watershed level is a more robust and accurate practice for long-term predictions and implementing soil conservation measures at the watershed level [15]. Therefore, in this study, we aim to estimate the soil erosion rate using empirical models in conjunction with field experiments in watersheds from the mid-hills of Nepal. We observed moisture, temperature, and electrical conductivity with the help of soil sensors installed 15 cm below the ground level and taking continuous measurements of various physical parameters in both watersheds. In addition, we examined the effect of precipitation on erosion under different land use practices in these mountainous terrains. A better understanding of erosion in these watersheds has broad interests due to its influence on diverse geological landscapes, soil types, vegetation, and climate change effects within a narrow geographic terrain in the Himalayan mid-hills sensitive to natural and anthropogenic changes [20]. In particular, the modeling output combined with the field experimental data from these mountainous terrains could help better assessment of sudden or incremental precipitation and land-use practices in these catchments.

Study area
We selected two watersheds, Tinahukhola and Aadhikhola (Fig. 1), representing the typical mid-hill terrains in Nepal. Tinahukhola watershed stretches along the Rupandehi and Palpa districts in Lumbini province, with an area of 1081 km 2 [24]. The altitudinal range varies from 85 to 1940 m above sea level. The watershed is drained by the Tinahukhola River, which is a monsoon-dominated perennial river originating from the southern slopes of the Mahabharata Range in the Palpa district [25]. The significant elevation difference has resulted in diverse climatic conditions ranging from sub-tropical warm to humid conditions. The average annual rainfall ranges from 1500 mm on the northern side and 2400 mm on the southern side of the watershed [26]. The major soil types within this watershed are fertile eutrophic cambisols and fluvisols [27]. Geologically, this watershed is mainly dominated by sedimentary and metasedimentary rocks upstream, while thick sandstone beds (1-45 m) are observed in the mid-section. The downstream Terai plains have rich alluvial sediments containing boulders, gravel, silt, and clay [28]. The hilly region of the Tinahukhola watershed is dominated by forests (62%), followed by agricultural land (33%); the floodplains in the Terai region have urban settlements with patches of croplands [20].
On the northern side of Tinahukhola lies the Aadhikhola watershed, with an area of 195 km 2 [24], located in the Syangja district of the Gandaki province. The rainfall in this region ranges from 2800 to 3200 mm; monsoons contribute around 80% of the precipitation [29]. Aadhikhola watershed includes soil types dominated by fluvisols, cambisols, and minor gleyic luvisols [27]. The temperature ranges from 5 to 32 • C [29]. This watershed lies in the Lesser Himalayan Zone in Nepal. The geology is dominated by non-fossiliferous, sedimentary, and meta-sedimentary rocks such as slate, phyllite, schist, quartzite, limestone, and dolomite, ranging from Precambrian to Eocene ages [28]. Like the Tinahukhola watershed, the Aadhikhola watershed has a dominant forest cover. This is followed by shrubland and cropland, covering nearly 27% of the area, whereas the built-up area, grassland, and water bodies cover only 7% of the total area [30].

Estimation of soil loss from RUSLE
The RUSLE model (Equation I) was used to determine the annual soil loss induced by slope runoff, specific cropping, and management practices [31]. The schematic framework of this assessment is shown in Fig. 2.
where, A denotes the annual soil loss, and provides value for the spatial and temporal average soil loss (tons ha − 1 yr − 1 ), R (MJ mm ha − 1 h − 1 yr − 1 ) is the rainfall-runoff erosivity factor, K (tons h MJ − 1 mm − 1 ) is the soil erodibility factor, LS (unitless) indicates slope length, gradient, C (unitless) represents land cover, and P (unitless) is conservation practices.

Rainfall erosivity factor (R-factor)
Rainfall erosivity represents the effect of rainfall volume and intensity as a main driving factor for water-induced soil loss. The rain erosivity factor was obtained based on the precipitation data from the Department of Hydrology and Metrology (DHM) using equation II [23]: where P = mean annual precipitation. This data was extrapolated to the whole watershed using the inverse distance weightage method [32].

Soil erodibility factor (K-factor)
The soil erodibility factor measures the susceptibility of the soil particles to detach and transport due to rainfall and runoff [22]. The erosion rate in a standard plot with a slope length of 22.13 m and a 9% slope gradient is the slope erodibility factor reflecting the soil loss based on the rainfall erosivity index [33]. In our study, we used the empirical relationship with soil properties from different soil samples collected from watersheds to calculate the soil erodibility factor [33]. The textural composition was obtained using the Bouyoucos hydrometer method [34]; soil organic matter was obtained using the dry combustion method [35]. Then, following the methods described earlier [36], the soil erodibility factor was obtained using Equation (3).
which can further be divided into the following equations: where, SAN= Sand percentage, SIL= Silt percentage, CLA = percentage of Clay, C = Soil Organic Carbon, SN1=Sand content subtracted from 1 divided by 100. The data obtained were extrapolated to the whole watershed using the Inverse Distance Weighing technique [37,38,39].

Slope length and steepness factor (LSfactor)
The digital elevation map was obtained from the ASTER platform. The spatial resolution of this image was 15 m, which was converted to 30 m resolution to match the cell size with the LANDSAT image (used for C-factor). The slope length denoted by L (equation IV) and slope gradient represented by S (equation V) were obtained previously for the mid-hills of Nepal [40]. The slope length and steepness map were obtained based on previous literature [41].

Land management factor (C-factor)
LANDSAT images (<10% cloud cover) from November 2018 were processed. While using the satellite images, we opted for less cloud cover (<10%) and synchronized other relevant data (e.g., rainfall, soil properties), which are used for calculating the RUSLE factors. Various corrections were included to account for atmospheric conditions, radiometric calibration, and geometric correction. The images were classified with the maximum likelihood algorithm. Various land-use classes were assigned to different C-factor values provided earlier and shown in Table 1 [3].

Conservation practice factor (P-factor)
The support practice factor provides the rate of soil erosion based on various cultivated practices in the study region. The value of the P-factor depends on the type of cropping practices done in the watershed, whether it is contouring, terracing, strip cropping, etc. Based on the land use practices, the slope values and the P-factor is assigned, ranging from 0 to 1, 0 being the least soil erosion and 1 being the highest rate of soil erosion [42]. In our study area, we mainly observed the practice of contouring. Thus, respective values were assigned as given in Table 1.

Observation of response of soil erosion to the change in rainfall
Based on the annual precipitation obtained from various sites within the study area, the average rainfall for five years was calculated from 1978 to 2017. The annual precipitation was used to calculate erosivity due to rainfall (R-factor) for the specific time interval. The R-factor was used in the RUSLE equation keeping the other factors (K-factor, LS-factor, C-factor, and P-factor) constant (same as in 2018), which yielded a change in soil loss due to variation in rainfall.

Experimental plots
We established experimental plots in three land-use settings (forests, rain-fed agricultural land, and irrigated agricultural land) to assess the impacts of rainfall on topsoil erosion (Table 2). Irrigated agricultural land is intensively utilized (3-4 crops a year) compared to rainfed agricultural land (1-2 crops a year). Therefore, data obtained from these plots are beneficial for comparing the erosion rates [19,43] calculated from the RUSLE model. Three replicate plots of dimension 6 × 3 m 2 were constructed in each land-use setting in both watersheds [44]. Manual funnel-type rain gauges were placed on each plot to measure rainfall. In addition, a gutter was placed, so the overland flow from the plots was collected into metallic drums to track the net discharge.
The data collection was done for the 2019 monsoon season. The measurements included precipitation amount and the overflow into drums after each rainfall event. The sample taken from the drum was transported to the laboratory to determine suspended solids. The relation between rainfall and the amount of soil transported during each rainfall event was selected from the specific land-use class for both watersheds. Further soil samples from each replicate and land-use plots from both watersheds were collected for bulk density estimation. The samples were collected at two depths (0-15 cm and 15-30 cm) to include top and subsoil profiles. Studies have been published to document the variability in topsoil erosion due to the differences in soil bulk densities [45]. In this study, we focused on rainfall and soil loss during the peak monsoons (~2 months), which coincides with the major duration for topsoil erosion in the mid-hills of Nepal.

Precipitation, overland flow, and soil loss
The funnel rain gauge measured the precipitation in the experimental plots. The manual readings were further verified with data obtained from DHM, Nepal. Precipitation contributes to subsurface flow (percolation and infiltration) and overland flow. The overland flow is responsible for soil transport; hence tanks were used to collect the overland flow. The water height in the tanks was used to estimate the overland flow volume.
Additionally, samples from the tanks were taken to estimate the total suspended load. Finally, the amount of sediment transported (lost) during each rainfall event was calculated from the suspended load and the overland flow volume.

Total soil loss = suspended load * volume of overland flow
After each sampling event, the drums were emptied and cleaned. All the drums, after cleaning, were ready to collect fresh runoff water from the plots during rainfall.

Soil water content, conductivity, and temperature
We installed EM-50 sensors in locations representing major land-use classes in the Tinahukhola watershed. These sensors were installed following guidelines provided by the manufacturer for recording soil moisture (cm 3 cm − 3 ), electrical conductivity (mS cm − 1 ), and temperature ( ⁰ C). These sensors use the dielectric method for measuring moisture, a thermistor for temperature, and an electromagnetic field for electrical conductivity. The sensors were programmed to take measurements every 6 h during the monsoons (see Supplementary data; Table 4). In addition, three sensors were placed in the selected land-use settings for automated measurements.

Revised universal soil loss equation (RUSLE) at plot level
To validate the RUSLE model output in these watersheds, we compared this with data generated from the experimental plots. The precipitation value obtained from the rain gauges installed within plots was used for the R-factor. Similarly, for the inputs of the Kfactor, soil samples were collected from the plots and analyzed in the lab; the values were used in equation II. The slope length factor (L) for these plots was 6, whereas the slope steepness was measured with a clinometer, and the value was used in equation V for calculating the slope steepness factor (S). For agricultural land, irrigated and rainfed plots were assigned a value of 0.21 for the Cfactor, whereas the forest was assigned a value of 0.03 [3]. The P-factor was based on the slope percentage calculated from the clinometer [42].

Model parameters
The R factor determines the soil erosion loss due to rainfall [46]. The average rainfall erosivity of the Aadhikhola watershed was estimated to be 1071 MJ mm ha − 1 h − 1 yr − 1 (Supplementary data; Fig. 1). The erosion rate is higher than the rainfall erosivity factor in the Tinahukhola watershed and is estimated as 733 MJ mm ha − 1 h − 1 yr − 1 . The K factor measures the vulnerability of soil particles to be eroded by rainfall. It depends upon soil properties like soil texture, organic matter present, soil structure, and permeability [47]. The soil erodibility factor in both watersheds is similar. It is 0.092 tons ha MJ − 1 mm − 1 in Aadhikhola and 0.096 tons ha MJ − 1 mm − 1 in Tinahukhola (Supplementary data; Fig. 2).
The topographic feature of the watershed is taken as the LS factor in the model. An increase in this factor denotes steepness and longer slope length, making the soil more prone to erosion. The long and steep slopes of the Aadhikhola watershed contribute to the high LS factor of 4.36, whereas it is 1.81 in Tinahukhola (Supplementary data; Fig. 3).
The land cover in any area is vital in affecting soil erosion. The site with relatively high vegetation cover has a low rate of soil erosion. In contrast, the area with less vegetation has a high rate of soil erosion. The C-factor is determined based on the land-use  Tables 1 and 2). The classification scheme is very reliable, with an overall accuracy of >80% and a Kappa coefficient of >0.80 [48]. Various C-factor values are assigned to different land-use classes, e.g., barren land exposed to high erosion is assigned a high C-factor value. In contrast, waterbody and built-up areas are assigned a value of 0. The C-factor in the Aadhikhola watershed was 0.11, whereas it was 0.13 for the Tinahukhola watershed (Supplementary data; Fig. 4).
The P-factor or the conservation practice factor refers to techniques practiced reducing soil erosion. For example, contour farming, strip cropping, and terracing are practices that affect soil erosion. The P-factor is similar in both watersheds, i.e., 0.95 in Aadhikhola and 0.93 in Tinahukhola (Supplementary data; Fig. 5). The P-factors were almost identical as both watersheds practice contour farming as protective measures.

Annual rate of soil loss
We combined all the factors incorporated in the RUSLE equation and calculated the soil erosion loss from both watersheds. Based on this, soil loss from the Aadhikhola watershed is 41.4 tons ha − 1 yr − 1 with a maximum value of 1302 tons ha − 1 yr − 1 and a minimum of 0.09 tons ha − 1 yr − 1 . The loss mainly occurred in the mid-western part of the watershed (Fig. 3). Similarly, soil loss from the Tinahukhola watershed is estimated as 24.1 tons ha − 1 yr − 1 with a maximum of 278 tons ha − 1 yr − 1 and a minimum of 0.04 tons ha − 1 yr − 1 . Again, the loss is mainly from the mid-northern part of the watershed (Fig. 3).

Precipitation amounts and soil loss in experimental plots
The relation between the amount of rainfall that occurred and soil loss in different land-use classes is summarized in Table 3. We captured 15 rainfall events (observations) in the Tinahukhola watershed during the 2019 monsoon season. The mean amount of precipitation recorded in three replicates of the irrigated agricultural plot is 867 ± 4.87 mm. The amount is comparable to the rainfed agricultural plot with a value of 898 ± 15.4 mm. However, the plots in forests recorded a lower amount (688 ± 1.74 mm). In the Aadhikhola watershed, 13 rainfall events were recorded. Like the Tinahukhola watershed, the plots in forests record less precipitation (450 ± 26.0 mm) compared to irrigated agricultural plot (458 ± 0.84 mm) and rainfed-agricultural plot (475 ± 8.02 mm).
The topsoil erosion, and the overland flow in experimental plots, show the trend: Irrigated agricultural land > Rainfed -agricultural land > Forests (Table 3)  observed is higher in Tinahukhola for the rain-fed agricultural plot and forests, the total soil loss in the irrigated agricultural plots at Aadhikhola is more than in Tinahukhola. This higher soil loss observed in Aadhikhola than in Tinahukhola for irrigated land might be due to the differences in agricultural practices between the two sites. Intense tilling and the absence of crop cover may influence soil loss in the watershed.
The RUSLE model when applied to the plot level showed a similar result to that of the observed value except for the irrigated agricultural land (Table 5). From our observation, we reported that for the Tinahukhola watershed, the value from the RUSLE model was slightly higher for the rainfed agricultural land compared to the observed value (difference = 2.31 tons ha − 1 ). Similarly, the difference of 4.854 tons ha − 1 was seen in the Aadhikhola watershed. On the other hand, the difference between the model and observed values is even smaller for the forested land use type. For the Tinahukhola watershed, the observed value was slightly higher than the model estimates (1.78 tons ha − 1 ). In comparison, for the Aadhikhola watershed, the observed value was slightly higher than the model estimates (0.344 tons ha − 1 ). A significant difference between the model estimation and the observed values were observed for the irrigated agricultural land use type in both watersheds. For the Tinahukhola watershed, the observed value is 11.12 tons ha − 1

Soil bulk density and EM-50 sensor records in plots
The results for soil temperature, moisture, and electrical conductivity from the EM-50 sensors at the Tinahukhola watershed are summarized in Table 4. Similarly, bulk densities at two different depths (0-15 cm and 15-30 cm) for both watersheds are summarized in Table 3. The temperature in the experimental plots at Tinahukhola is rainfed agricultural plot (25.48 • C) > forests (22.65 • C) > irrigated agricultural plot (20.55 • C). The moisture content is higher in the irrigated agricultural land (0.35 m 3 m − 3 ), followed by the forests (0.27 m 3 m − 3 ) and rainfed agricultural land (0.22 m 3 m − 3 ). While the electrical conductivity is high in the irrigated agricultural land (0.06 mS cm − 1 ), it is comparable to the rainfed agricultural plot (0.01 mS cm − 1 ) and forests (<0.01 mS cm − 1 ). Fig. 4 shows the relation between moisture and temperature. The temperature and electrical conductivity change for different landuse classes with the increase in humidity. The increase in moisture content leads to a decrease in temperature for rainfed agricultural land. Electrical conductivity rises, causing the moisture to increase to 0.25 m 3 m − 3 before decreasing again. For irrigated agricultural land, electrical conductivity increases linearly with the increase in moisture content, while the relation with temperature is unclear. Both temperature and electrical conductivity increase with the moisture content in forests.
Bulk density provides compactness of the soil, which in turn helps to determine soil disturbance and vulnerability. The bulk densities at 0-15 cm depth show a clear trend in both watersheds ( Table 3). The irrigated agricultural land has the highest bulk density, followed by rain-fed plots and forests. In soils at a depth of 15-30 cm, bulk density for all three land-use classes is higher in Aadhikhola

Discussion
Among the different soil erosion models available, RUSLE is more widely used than, e.g., USLE and MUSLE because of its versatility and broad usability. RUSLE has been progressively developed and includes an updated research and experimental module to incorporate new data [38]. Improvement in the diagnostic parameters, along with corrections of various inaccuracies, has been made in this model [31,39], resulting in better assessment. For example, there is a correction of the R-factor in the USLE model where if ponded water exists, the model can accurately assess soil cover and management processes, including the influence of topography. The revised model also accounts for soil erodibility due to changes in freeze/thaw and moisture content. In addition, empirical equations have been developed in this model by integrating the updated values of soil conservation practices [49]. Because of these advantages, the application of RUSLE is suitable not only at the local watershed level but also helpful in assessing soil loss at the regional and national levels [22,50].

Estimation of soil loss from the RUSLE model
In the Aadhikhola watershed, we recorded an average rainfall of 2638 mm annually, whereas, in the Tinahukhola watershed, the rainfall was 2010 mm (for 2018). In both watersheds, the average precipitation is higher than the average rainfall in Nepal (1750 mm yr − 1 [43]). High rainfall contributes to rapid erosion (~25 tons ha − 1 yr − 1 ) in these watersheds. Additionally, both watersheds are in the hilly region of the mid-range Himalayan mountain chain. The hilly areas with steep slopes aid in the detachment of topsoil, making it easier for downhill transport. This further increases the amount of soil loss in the hilly regions.
In one of the studies using the RUSLE model in the Chitral district, soil erosion has been estimated as high as 78 tons ha − 1 yr − 1 , and the topographical factor has been identified as the major reason [51]. The authors predict the rate of soil erosion to be worse in the future with changes in land use and climatic conditions. Under RCP (Representative Concentration Pathway) 8.5, it has been predicted that the rate of soil erosion in the same region may reach up to 550 tons ha − 1 yr − 1 during 2030-40 [52].
In Nepalese watersheds, variable rates of soil erosion have been documented. A spatial assessment done in the hilly watersheds of Western Nepal reported that the average soil loss in the Aringale watershed in Salyan is 17 tons ha − 1 yr − 1 (ranging from 0.03 to 100 tons ha − 1 yr − 1 ). The low value is due to less precipitation in the region, contributing to the low R-factor (365-435 MJ mm ha − 1 h − 1 yr − 1 ) than in Aadhikhola (1071 MJ mm ha − 1 h − 1 yr − 1 ) and Tinahukhola (733 MJ mm ha − 1 h − 1 yr − 1 ) watersheds. Similarly, a recent study estimated soil erosion ranging from 0 to 273 tons ha − 1 yr − 1 , which on average is ~38.4 tons ha − 1 yr − 1 for the mid-hills [22]. Consistent with this, the estimates for soil loss in the Aadhikhola watershed in this study are similar. The Gandaki basin recorded an erosion rate of ~30.7 tons ha − 1 yr − 1 . Aadhikhola watershed falls within the larger Gandaki Basin. In contrast to Aadhikhola, the Chure/Siwalik region covered by the Tinahukhola watershed indicate a low soil erosion rate (~6.9 tons ha − 1 yr − 1 ). The Terai floodplains in the Tinahukhola watershed have the lowest erosion rate (~0.1 tons ha − 1 yr − 1 ). It is suggested that for young mountain chains, the rate of soil erosion up to 25 tons ha − 1 yr − 1 is acceptable [53]. Erosion rates in our watersheds are very close to this number or exceeded the proposed limit (>25 tons ha − 1 yr − 1 ), highlighting the vulnerability to soil loss and the immediate need for protective measures that are needed in this region. In particular, most residents in the mid-hills are farmers whose livelihood depends on agriculture; excessive erosion of topsoil from arable land will directly affect their livelihoods.
Using cosmogenic nuclides, an erosion rate equivalent to 1.6 mm/year (average for 15 years) in the Narayani basin was estimated [54]. The authors indicated that this area has a glacial landscape and limited topsoil thickness and supply. Hence, soil erosion contributes ca. 10% of total erosion. However, this study covers large river systems, e.g., Kaligandaki, Trishuli, Marshyangdi, and Narayani, all of which are glacial-fed. The small fluvial watersheds investigated in this study are, however, different. These small watersheds have significant proportions of cultivable land with variable soil depths and gentle slopes. These watersheds are impacted intensely by anthropogenic (land use changes) and natural changes (e.g., erosion and precipitation). Notably, the intense monsoon showers cause more soil erosion than sheet and rill erosion [54]. In addition, landslides during heavy monsoons could contribute to additional soil loss and discharge into these rivers. Although RUSLE underestimates the erosion rate in watersheds with frequent landslides, estimation from this model provides valuable insights. Nonetheless, additional modules must be incorporated into the model to comprehensively assess landslide impacts, mass wasting, and complex erosional processes in heterogenous landscapes [55].
This study divided the watershed into various severity classes based on the soil erosion values ( Table 6). The different classes followed the established range [56] and included: very slight (<5 tons ha − 1 yr − 1 ), slight (5-15 tons ha − 1 yr − 1 ), moderate (15-30 tons ha − 1 yr − 1 ) , severe (30-50 tons ha − 1 yr − 1 ), and very severe (>50 tons ha − 1 yr − 1 ) soil loss [57]. Table 6 summarizes the severity of the soil erosion in both of the watersheds in this study. In the Tinahukhola watershed, nearly 50% of the land belongs to the low erosion class according to the soil loss classification scheme [57]. These areas are mainly in the southern part of Nepal, covering the flat Terai plains, and gentle slopes account for low soil erosion. In contrast, nearly 20% of the land cover is affected by severe soil erosion. Such areas are located in the sloping hilly regions in the northern part of the watershed. Thus, in the Tinahukhola watershed, with diverse geomorphological characteristics, experiences a wide range of soil erosion intensity. However, we notice that the Aadhikhola watershed is more vulnerable to soil loss than the Tinahukhola watershed. The Aadhikhola watershed has only 7% of its land belonging to the low soil loss category, while nearly 70% of its land belongs to the class severe and very severe category. This clearly shows the higher vulnerability of soil erosion in the Aadhikhola watershed, implying immediate soil conservation needs in the region.

Organic matter and erosion
The organic matter in soil increases infiltration rates of precipitation through soil horizons. The high infiltration rates decrease surface runoff and reduce soil erosion and the K-factor by forming stable soil aggregates with organic matter [58,59]. In general, loss of soil organic matter due to degradation (from land-use practices) can significantly alter its physical properties, namely pore size, bulk density, aggregation, and water-holding capacity in soils [60]. Therefore, we estimated the OC flux for both watersheds [20]. The annual transfer of OC from Tinahukhola is higher than Aadhikhola. However, after normalization concerning their size, Aadhikhola has a higher OC flux (1.67-1.07 ton km − 2 yr − 1 ) than Tinahukhola (1.24-1.05 ton km − 2 yr − 1 ). The erosion rates estimated by RUSLE corroborates trends in suspended sediment load and POC flux [20] in these watersheds. However, the annual suspended sediment load for these watersheds is low compared to the RUSLE values because in those research bedload portion has not been included.
Most of the OC flux is reported from input derived from various terrigenous sources, including old soil. Thus, the K factor in catchments could have variable results based on the parent soil type [58]. They conclude that the K-factor depends mainly on texture in granitic terrains, whereas in limestone-derived soils, OM plays a crucial role. Since both watersheds are carbonate (limestone) dominated sites, OM is most likely the principal factor influencing the K factor when using the RUSLE model for estimating soil erosion.

Change in soil loss and rainfall
In both watersheds, there is clear evidence that rainfall is a significant driver for regulating erosion rates. For example, the Tinahukhola watershed depicts a trend coinciding proportionally between rainfall and erosion rates from 1978 to 2017 (Fig. 5B). A similar trend exists in the Aadhikhola watershed except for a few years when this relationship does not hold (Fig. 5A). For example, the annual precipitation increased from 2468 mm (1978)(1979)(1980)(1981)(1982)  Even though the amount of rainfall and rate of soil erosion seems to increase, this is not significant for both Aadhikhola (F significance = 0.56 > α = 0.05) and Tinahukhola (F significance = 0.98 > α = 0.05). In this study, while the amount of precipitation is accounted for, the weather stations did not report its intensity. In this context, we refer to a recent investigation [56], whereby the authors indicate that increased precipitation intensity coincides with a higher soil erosion rate.

Soil loss from erosion plots
Due to the difference in temporal coverage, we observe higher precipitation in Tinahukhola (~900 mm during 15 observations). In contrast, in Aadhikhola, this is ~500 mm during 13 observations (Table 4). Even though there is a difference in precipitation in both watersheds, the amount of soil loss is very close. Notably, soil compaction affects soil loss [61], and bulk density in the soil correlates positively with soil loss [62]. In this study, bulk density in the Aadhikhola watershed is higher than in the Tinahukhola watershed for all three land-use activities. Hence, it is likely that differences in the physical property of soil affected the variable soil loss observed in the catchments.
Similarly, various land-use practices impact soil loss in both watersheds. The highest amount of soil loss occurs in the irrigated agricultural plot, followed by the rainfed agricultural plot (Table 4). Forests have the least amount of soil loss. The irrigated agricultural land is highly disturbed due to frequent tillage during the crop cycles. Thus, human activities and land use changes affect soil quality and increase the rate of soil erosion, which is also prevalent in other parts of the mid-Himalayan region in Nepal [63].
Since the plots with different land use are in proximity (within 1 km 2 ), the decrease in rainfall recorded in the forests should be due to the interception of precipitation by vegetation cover, specifically its crown. The interception not only reduces the amount of rainfall reaching the ground but also reduces the intensity with which the raindrops hit the ground or, in short, reduces the energy in the raindrop. This, in turn, decreases soil detachment and reduces soil loss.
Even though the RUSLE model is quite reliable at estimating soil loss, it could not depict the local conditions, especially for the irrigated agricultural land in our case. We saw considerable differences between the estimated and observed values. For both watersheds, we found that the observed values were higher than the estimated values for irrigated agricultural land. The main contributing factor to this might be the land use management factor. While the irrigated and rainfed agricultural land has different intensities of land use, both were assigned the exact value of the C-factor, which is 0.21. The irrigated agricultural land faces more intense agricultural practice with 2-3 crops in a year, whereas the rainfed agricultural land has only 1-2 crops in a year. This means frequent the tillage is done in irrigated agricultural land compared to rainfed croplands. Thus, we anticipate that the soil is more fragile in the irrigated agricultural land.

Physical properties of soil and their implication in soil loss
Variable soil characteristics present in different land-use types lead to variations in the rate of soil erosion. It is well known that different land-use practices affect physical conditions in soil [45]. The soils with a higher level of disturbance, i.e., soils from agricultural land, have a low quality (high bulk density). In contrast, soils with minor disturbances have a better quality (low bulk density) [45]. Consistent with this observation, we found similar results in our study. The irrigated agricultural land with a high level of disturbance has the highest bulk density in the topsoil, followed by the rainfed agricultural land. In contrast, forests have the lowest bulk density, signifying less anthropogenic disturbance. Another study documented that physical properties in soils affect soil loss [61]. Soils with high bulk density have high soil erosion. Consistent with this, irrigated agricultural land in our study has the highest rate of soil erosion (Table 4), followed by rainfed agricultural land; the lowest rate of soil erosion occurs in forests. The disturbance and soil bulk density followed the trend: irrigated agricultural land > rainfed agricultural land > and forests.
The irrigated agricultural land has a relatively higher amount of moisture content. The high moisture content in the soil makes it more vulnerable to soil loss. The high moisture content blocks all the pores resulting in water-logged conditions, and preventing the percolation of rainwater. This results in more surface runoff. The high surface runoff, in return, transports more soil from the catchment. In this study, with the increase in moisture content, the temperature declined, whereas the electrical conductivity increased. There has been the practice of using a massive amount of chemical fertilizers to increase the productivity of paddy [64]. This may be the major contributor to the observed high levels of electrical conductivity, i.e., 0.06 mS cm − 1 (Table 3). Temperature decreases with the increase in moisture content, and less vegetation cover increases the evaporation rate [65]. These factors have led to less moisture content in the rainfed agricultural soils (0.22 m 3 m − 3 ).
Even though forests have low average moisture content and temperature than irrigated agricultural land, the value is higher than that of rainfed agricultural land. The forest soil is found to have a higher water-holding capacity, at least up to a depth of 40 cm [, [66,67]]. Although forest soil has a high moisture content, soil loss is expected to be high due to anthropogenic disturbances in the rainfed agricultural plots. In both watersheds, the sloping area covered by forests is more than 60% [20]. Even with most of the area under forest cover in both watersheds, soil loss exceeds the acceptable limits. Moreover, urbanization trends in Nepal indicate that forests are rapidly converted into built-up areas and/or agricultural croplands. Thus, land use further increases the vulnerability to soil loss in these mountainous terrains. As both of our watersheds are vulnerable to extreme rainfall events, steep slopes, and landslides, the changeover in land use further increases topsoil erosion in these watersheds.

Conclusions
In this study, we identified areas with a high vulnerability to soil erosion in the mid-hills of Nepal. We analyzed the soil loss from two small representative watersheds using the RUSLE model and experimental soil erosion plots. The experimental field plots were set up in areas with three different land-use activities, i.e., irrigated agricultural land, rainfed agricultural land, and forests in both watersheds. Soil loss is higher in the Aadhikhola watershed than in the Tinahukhola watershed. The primary reason is the high precipitation in Aadhikhola and the steep mountain slopes. Based on land use, irrigated agricultural land had the highest rate of soil erosion, followed by rainfed agricultural land. The forests with high vegetation cover have the least amount of soil loss. The soil loss correlated positively with compactness in the soil. Soils with higher bulk densities are more prone to soil loss. We tested whether the amount of soil loss has changed due to the change in rainfall patterns over the last 40 years, but the results are statistically insignificant. Finally, while the RUSLE model can estimate soil loss with less bias for forests and rainfed agricultural land, the model tends to underestimate soil loss from irrigated agricultural land. Thus, we recommend revising the model that captures the local conditions and farming practices in regions like Nepal. This would allow relevant stakeholders and policymakers to plan for better soil conservation strategies in erosion-prone areas in the mid-hills of Nepal. For examples, a plantation in a barren land, reduced tillage practice, terrace farming in sloping lands, etc. are recommended as simple yet most effective measures for reducing soil erosion.