Land use change affects water erosion in the Nepal Himalayas

Soil erosion is a global environmental threat, and Land Use Land Cover Changes (LUCC) have significant impacts on it. Nepal, being a mountainous country, has significant soil erosion issues. To examine the effects of LUCC on water erosion, we studied the LUCC in Sarada, Rapti and Thuli Bheri river basins of Nepal during the 1995–2015 period using the Remote Sensing. We calculated the average annual soil loss using the Revised Universal Soil Loss Equation and Geographical Information System. Our results suggest that an increase in the agricultural lands at the expense of bare lands and forests escalated the soil erosion through the years; rates being 5.35, 5.47 and 6.03 t/ha/year in 1995, 2007 and 2015, respectively. Of the different land uses, agricultural land experienced the most erosion, whereas the forests experienced the least erosion. Agricultural lands, particularly those on the steeper slopes, were severely degraded and needed urgent soil and water conservation measures. Our study confirms that the long term LUCC has considerable impacts on soil loss, and these results can be implemented in similar river basins in other parts of the country.


Introduction
Soil erosion is a severe ecological issue that humanity is facing [1] as it washes away the fertile topsoil, deteriorates soil quality and increases the soil sediments in stream channels [2]. Extensive use of available land for agriculture increases the soil loss at a global scale [3]. It is one of the major environmental issues of hill and mountain ecosystems [4]. Soil erosion is a significant feature in Nepalese terrain, given the hilly topography and rugged mountains, much-concentrated rainfall events in the monsoon (June-September) and increased human influence in the removal of natural vegetation [5,6]. Numerous studies have been performed to assess the soil erosion in Nepal, mostly in the Middle Mountain region [7] but only a few have addressed the erosion in the High Himalayas [4,8,9]. Variation in topography, slope, land use patterns and population pressure across different physiographic regions produces different rates of soil erosion in Nepal, ranging from zero in the lowland areas to 420 t/ha/year in the shrublands [10]. Soil erosion rates of 11.17 and 10.74 t/ha/year have been reported in the Aringale Khola Watershed and Sarada river basin, respectively in the Siwalik Hills [11,12]. Erosion due to natural causes is the a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 and 81˚49'17'' and 83˚16'49'' E longitudes and has an area of 12,182 km 2 (Table 1). It encompasses five distinct physiographic regions of Nepal: Terai Plains (516-700 m), Siwalik Hills (700-1,500 m), Middle Mountains (1,500-2,700 m), High Mountains (2,700-4,000 m) and High Himalayas (4,000-7,264 m) and covers the area of Jajarkot, Rukum and Rolpa districts and some parts of Surkhet, Dang, Pyuthan, Salyan, Baglung, Myagdi, Dolpa, Dailekh, Kalikot, Jumla and Mugu districts. Of the three river basins, Thuli Bheri is the most diverse as it encompasses all the physiographic regions available, with elevation ranging from 523 m in the South to 7,264 m in the North. It is also the largest river basin amongst the three, covering an area of 9,338 km 2 followed by Rapti and Sarada river basins with 1,972 and 872 km 2 , respectively.
The data used in the analysis are the rainfall measurements (1980-2016) from 53 rainfall stations (Please see S1 File) obtained from the Department of Hydrology and Meteorology,  For all the images (TM and OLI), radiometric calibration and Fast Line-of-sight Atmospheric Analysis of Hypercubes atmospheric correction were performed to remove the atmospheric influence, and dark object subtraction method was used to remove the effects of smoke, dust and haze [29]. The LULC maps were developed using the Maximum Likelihood Classification (MLC) because of its strong theoretical base and simplicity [29]. MLC is one of the most preferred parametric techniques as it calculates the likelihood of an unknown vector based on the highest probability of fit based on Bayesian equation [29,30]. Accuracy assessment of all the classified images was carried out using the points collected from the satellite images. High-resolution satellite images were available for 2007 and 2015 (QuickBird and WorldView) dataset; however, the 1995 dataset did not have high-resolution images. Therefore, the Google Earth image information was supplemented through the author's knowledge of the area as well as expert input. The overall producer and user accuracies were obtained for each classified image, and the Kappa coefficient [31] was calculated to check the accuracy of classified images.

Soil erosion estimation
We used the RUSLE model [32] to calculate the average soil loss in the study area using the following equation: Where E is the estimated average soil erosion (t/ha/year), R is the rainfall factor (MJ mm/ha/h/ year), K is the soil erodibility factor (t ha/MJ/mm), LS is the combined slope length and slope steepness factor (dimensionless), C is the cover management factor (dimensionless), and P is the support practice factor (dimensionless).
The R factor represents the potential erosivity of soil [33] and relies on the intensity and volume of rainfall occurring in a place over a specified period [34]. The R factor was estimated using the following equation [35]: Where r is annual rainfall in mm The K factor represents the soil susceptibility to detachment caused by the beating action of rainfall and runoff water [32,36]. We acquired soil textural maps of the study area from the National Land Use Project, Nepal, and prepared the K factor map assigning values to soil texture as proposed by [37] ( Table 2).
The slope length (L) and slope steepness (S) are the topographic factors, and they account for the effects of slope length and steepness on soil loss [36]. We calculated the LS factor using the following equation given by Wischmeier and Smith [38] and Šurda, Šimonides [39]. Where Cell size = Grid cell size (20 m for this study), m = 0.2 to 0.5 (0.2 for slopes less than 1%, 0.3 for 1-3%, 0.4 for 3-4.5% and 0.5 for slopes exceeding 4.5%) [40] and 's' is the slope in percentage.
The C factor represents the relationship between erosion on bare land and erosion on cultivated land [41] whereas the P factor represents the consequences of soil protection measures such as contouring and terracing [12]. The C and P values were assigned as per the land uses in the vector format [42][43][44][45] (Table 3)  We have also established corn and bare fields as soil erosion experimental plots in the study area to facilitate the comparison of the soil erosion rates computed from the RUSLE modelling. Replicated four times each, five treatments: no-tillage + mulch, no-tillage + no mulch, tillage + mulch, tillage + no mulch and bare fallow each with 3 × 3.75 m plot size, were established in an unbalanced complete randomized block design (Fig 2). Manakamana-3, a popular corn variety in the Nepalese hills, was sown with a spacing of 75 cm × 25 cm, and trench of 50 cm was dug at the lower end of experimental plots to collect the soil runoff. Corn was grown in all the treatments except in the bare fallow for two years. The primary objective of this experiment was to assess the impacts of tillage and mulch on soil erosion and corn yield.

Spatial interpolation
Kriging, a linear geostatistical interpolation technique, was used to prepare the digital map layers for R and K factors of the RUSLE. The acquired data were fitted to spherical semivariogram and then interpolated with the ordinary kriging technique in ArcMap. Kriging is a vital tool to estimate values at non-sampled places based on the sampled data [46,47]. It is an interpolation method which approximates the value of a function at a given point as a weighted amount of the function values at the adjacent points [48].

LUCC dynamics
Our classification accuracies were found to be 85% for 1995, 86% for 2007 and 84% for 2015 with Kappa Coefficient of 0.81, 0.82 and 0.79 for the year 1995, 2007 and 2015, respectively. The relative proportion of different LULC classes during the study periods and percentage change in LULC through the years are shown in Table 4 and Fig 3. The land use change matrix of the study area is presented in

RUSLE factors and soil erosion maps
The spatial distribution of R, K, LS, C and P factors are given in Figs 6 and 7, whereas the erosion maps for the years 1995, 2007 and 2015 are presented in Fig 8. As expected, soil erosion had increased through the years; mean erosion rate of 5.35 t/ha/year in 1995 had increased to 5.47 and 6.03 t/ha/year in 2007 and 2015, respectively. Soil erosion rates in the study area were regrouped into eight severity classes, and the respective areal coverage was presented according to guidelines developed by Uddin, Abdul Matin (9) (  (Table 7), and the erosion rates were low in the High Himalayas with values of 1.10, 0.07 and 0.86 t/ha/year in1995, 2007 and 2015, respectively. Increasing soil erosion rates through all the physiographic zones, except the High Himalayas, were experienced over 21 years.

Discussion
The observed changes in the different LULC classes throughout 1995, 2007 and 2015 can be explained with one of the following reasons. First, the increase in population by 43% in 2011 as compared to 1991 [49] had significantly increased the built-up areas during the 1995-2015 period. Increase in the agricultural area (52.87%) at the expense of bare lands and forests through the years can also be attributed to the population rise. Second, increasing awareness of community-based forestry proved its significance in protecting forests throughout the country [50]. Even though an overall decrease of 11.17% in the forests was observed during the study period, there was a rise in the forest area by 18.55% during the 2007-2015 period (Table 4);   2000 [50] have become a substantial success in protecting forests throughout the country, and the study area is no exception regarding this. Overall, agricultural lands, built-up area and snow cover had increased, resulting in the subsequent decline in the bare lands, forests and water bodies. The increasing population had increased the built-up areas, but the migration of people abroad or to urban centres in search of better job opportunities have left some agricultural areas uncultivated. Land Use Policy (2013, 2015) has banned the illegal conversion of one land use to another [51]; however, partitioning and illegitimate conversion of agricultural lands for commercial purposes are becoming a major issue for the society. This practice also exacerbates soil loss. An inspection of the change matrix (Table 5) shows an increase in the agricultural areas at the expense of bare lands and forests during the 1995-2007 period. This change has led to an increase in the soil erosion in the area for the described period as the erosion rates are higher in the agricultural lands in comparison to the bare lands and forests. Similarly, a decrease in the bare lands with a corresponding increase in the forests (15%) and agricultural lands (8%) was reported for the same period. Since the erosion rates in the forests were meagre (0.3 t/ha/ year), increase in the forests contributed very little to the decrease in soil loss, whereas the agricultural lands having higher erosion rates significantly increased soil loss during the 1995-2007 period. Again, an increase in the agricultural lands with a corresponding decrease of bare lands (38%) and forests (8%) was observed during 2007-2015. Snow cover had reduced by 2015, which also had a marked increase in soil erosion during the period. An overall increase in the agriculture, built-up areas and snow cover and decrease in the bare lands, forests and water bodies were observed during 1995-2015. Soil erosion rates in the forest, water bodies and snow were negligible. Although with higher erosion rates, built-up areas occupied less area, so the contribution of the built-up area in determining the rate of soil loss was also negligible. Thus, it can be concluded that the increase in the agricultural lands at the expense of bare lands and forests escalated soil loss in the area.
Not only the erosion rates have increased through the years, but also an increase in severely eroded area by 20% and a decrease in the less severe area by 4% were reported. Soil erosion was the highest in the Siwalik Hills, ranging from 8.77 to 9.84 t/ha/year, and the lowest in the High Himalayas (0.07 to 1.10 t/ha/year). Presence of erodible silt particles and occurrence of steep slopes have made the Siwalik Hills the most eroded zone across the study area. The second most eroded zone was the Middle Mountains followed by the High Mountains, Terai Plains and High Himalayas. Gently sloping flatlands in Terai Plains [52] and presence of snow cover in the High Himalayas make these physiographic regions less susceptible to soil erosion. Agriculture was the most eroded land use in the study area. Lack of vegetation cover is the main reason for higher soil losses in the agricultural lands [53]. Erosion rates through agricultural fields had reduced through the years, but with the upsurge in agricultural lands, the contribution of agriculture to gross soil erosion ultimately increased. Surprisingly, bare lands lost soil at lower rates than the agricultural fields; this may be due to the continuous and excessive tillage in the crop fields [14] while soil may have stabilized or compacted to resist the soil erosion to some extent in the bare lands. Community-based forestry showed positive impacts in safeguarding the forest resources after 2007.
The RUSLE has normally been utilized to estimate soil erosion in gentle slopes across the world. However, it has been increasingly used to assess soil loss across all the slopes ranging from gentle to steep slopes [54,55]. The RUSLE has been extensively used in Nepal as well to calculate the erosion rates in areas with steep slopes too [9,40,56].
Although the results come with several future possibilities, this method has some limitations too. Some improvements can be made in the calculation of the RUSLE factors. In this research, annual rainfall data were used to calculate the R factor, whereas more intense rainfall events which are more likely to have a marked impact on soil loss [40] were not considered. Erosion potential of a given rainstorm is calculated by multiplying kinetic energy of the rain with maximum 30-minute rainfall intensity. Since the meteorological stations of Nepal do not have adequate laboratory facilities to measure 30-minute rainfall data, they were not considered in calculating the R factor. However, there is a possibility of improving the R factor in the future. Similarly, better estimates of the K factor can be obtained if soil texture (sand, silt and clay %) and soil organic matter data are available instead of just the soil textural classes. It also lacks the automation of the procedure in calculating the RUSLE factors.
It would be much useful to verify the soil erosion estimates with real soil loss observations from long term soil erosion plots. However, the availability of only a few soil erosion location data restrained us from validating the estimates and, thus, analyzing the errors further. We compared the soil erosion rates computed from the RUSLE with plot-based erosion measurements from another experiment conducted in the study area. The two-year data from the field plots produced annual soil loss rates in the range of 9-10 t/ha/year (approx.) which is higher than the RUSLE derived soil loss. We confirmed the synergistic interaction of mulch with tillage to lower the losses of soil organic matter and total nitrogen and the effectiveness of no-till to reduce the soil losses in the study area as well. Erosion plots include both the corn and bare fields; where all the plots except the bare ones received tillage/no-tillage and/or mulching/no mulching practice. The mean soil loss in the study area will be lower than what is reported in the erosion plots as a large part of the study area are forests, and soil loss will be much lower in the forests.
This study used a modelling approach to develop spatial distribution maps of water erosion for 1995, 2007 and 2015 using the RUSLE, GIS and Remote Sensing. Spatial distribution of soil erosion maps, such as produced here, can be a vital tool for planning, particularly where the local economy is mainly agriculture-based and rapid urban development has taken over the agricultural lands [57], as in this instance.

Conclusions
This study provided information on the changes in LUCC in Sarada, Rapti and Thuli Bheri river basins of Nepal over the 1995-2015 period and analysed the effects of these changes on the rates of soil loss in the study area. The RUSLE was used in conjunction with GIS and Remote Sensing taking three temporal Landsat imageries from 1995, 2007 and 2015. The results suggest that the long term LUCC has a significant impact on soil loss in the study area as soil loss rates have increased over time. Increase in the agricultural lands at the expense of forests and bare lands have increased soil erosion. Of the different land uses, erosion was the highest in the agricultural lands. Forest cover seemed much more effective in reducing soil loss, even at the steeper slopes too. Soil erosion rates were the highest in the Siwalik Hills and Middle Mountains as compared to other physiographic regions. Proper agricultural management is needed to reduce the soil loss in the agricultural lands, focusing more on the steeper slopes. The use of the RUSLE, GIS and Remote Sensing to estimate the soil erosion dynamics over other soil erosion estimation models has many advantages as it saves a substantial amount of time [12] and also avoids tedious measurements of soil sediments over time. The analysed land use maps and estimated soil erosion severity maps may be helpful for policymakers and planners in developing better soil and water conservation programs across the country. Since LUCC is inevitable in the future, sustainable land use plans should be formulated to keep soil erosion rates under control so that it does not pose added threats to the ecological sustainability of the area.