Soil erosion estimation using Geographic Information System (GIS) and Revised Universal Soil Loss Equation (RUSLE) in the Siwalik Hills of Nawalparasi, Nepal


 Soil erosion is one of the gravest environmental threats to the mountainous ecosystems of Nepal. Here, we combined a Geographic Information System (GIS) with the Revised Universal Soil Loss Equation (RUSLE) to estimate average annual soil loss, map erosion factors, compare soil erosion risks among different land use types, and identify erosion hotspots and recommend land use management in the Girwari river watershed of the Siwalik Hills. The annual soil loss was estimated using RUSLE factors: rainfall erosivity (R), soil erodibility (K), slope length and steepness (LS), cover crops (C), and conservation practices (P), and erosion factors maps were generated using GIS. Results indicate highest total erosion occurring in hill forests (13,374.3 t yr–1) and lowest total erosion occurring in grasslands (2.9 t yr–1). Hill forests showed high to very severe erosion due to steepness of hills, open forest types, and minimal use of conservation practices. Also, erosion hotspots (>15 t ha–1 yr–1) occurred in only 4.2% of the watershed, primarily in steep slopes. Overall, these results provide important guidelines to formulate management plans and informed decisions on soil conservation at local to regional levels. While the study is the first effort to assess soil erosion dynamics in the Girwari river watershed, potential for application in other basins largely exists.

factors affecting soil vulnerability to erosion in Nepal are concentrated rainfall, poor soil structure, undulated topography, deforestation, excessive tillage and grazing, and crop residue removal (Chalise et al. a, b). Despite the increasing threat to soil quality, resilience, and services, erosion modeling in small watersheds of Nepal is still in its infancy, thereby affecting conservation plans and policies. RUSLE is a simple and robust model, and has significant advantages in watershed scale erosion and sediment transport modeling (Merritt et al. ). RUSLE offers several improvements in the USLE factors to represent diverse land and crop management, and slope forms (Renard et al. ). Briefly, R values are minimized under flat slopes experiencing high rainfall intensity. The K factor has been greatly improved to account for seasonal variability (e.g. freezing and thawing, soil moisture, and soil consolidation), and to include repercussions of rock fragments on soil permeability and runoff. RUSLE also includes functions for slope steepness and for soil vulnerability to rill erosions relative to inter-rill, and generates a rather linear relation for slope steepness compared with the USLE model. It utilizes a sub-factor method consisting of historical land use, canopy, ground cover, and within-soil effects to improve estimates of weighted average soil loss. Additionally, P factor values are improved for contouring, terracing, stripcropping, and rangeland conservation and management practices (Renard et al. , ). In a recent study, the RUSLE model was used to study the spatial distribution of soil erosion in Nepal (Koirala et al. ), which estimated a nation-wide mean annual soil loss of 25 t ha À1 yr À1 . The mean soil erosion rate under different land use occurred in the following order: barren land (40 t ha À1 yr À1 ) > agricultural land (29 t ha À1 yr À1 ) > shrubland (25 t ha À1 yr À1 ) > grassland (23 t ha À1 yr À1 ) > forests (22 t ha À1 yr À1 ). Also, Uddin et al. However, erosion mapping in the Girwari river watershed of Nepal is still lacking. To that end, we aimed to estimate the spatial distribution of soil erosion in the Girwari river watershed using remote sensing data, the RUSLE model, and GIS.
The study is the first attempt to assess soil erosion dynamics in the Girwari river watershed using high-resolution remote sensing data and field measurements. The specific objectives were to: (i) estimate the average annual soil loss in the Girwari river watershed of the Siwalik Hills, (ii) map erosion factors, (iii) compare soil erosion risks among different land use types, and (iv) identify erosion hotspots that require immediate attention for soil and land use management.

Study area
Nepal is a mountainous country surrounded by India on the east, west and south, and the Tibetan region of China to the north. It encompasses a total area of 147,516 km 2 .
The country is ecologically divided into three regions: Terai ( , ). As such, this watershed was selected for the study due largely to increasing agricultural activities, diversity in land use and terrain landscapes, and vulnerability to soil erosion.

Data description and sources
Data on rainfall, topography, soil properties, and land use and land cover were obtained from multiple sources and through field measurements. A DEM of 20-meter spatial resolution was generated using contour and spot height data from the Department of Survey, Nepal, using GIS software ArcMap (Version 10.2) ( Figure 2). The DEM was used to calculate the LS factor in the RUSLE model. Five years (2010)(2011)(2012)(2013)(2014) of rainfall data were obtained from the Department of Hydrology and Meteorology, Nepal, to calculate the R factor. Soil physical and textural properties were determined from field sampling and laboratory analysis, and used to calculate the K factor. The Normalized Difference Vegetation Index (NDVI) values were generated and mapped from a LANDSAT8 OLI (2015) to determine the C factor. Similarly, the P factor value was obtained from literatures and through field observation. The conceptual framework for soil erosion is explained in Figure 3.

Rainfall erosivity factor (R)
Rainfall-runoff erosivity is defined as the intrinsic capacity of rain to cause erosion, and is computed for either a single or a series of storm events. The R factor is also defined as a product of total rainfall energy (E) and the maximum 30 min rainfall intensity (I 30 ). The R factor and detachment of particles and splash largely depends on the amount, intensity and terminal velocity of rain, raindrop size and distribution, and soil conditions (texture, aggregates and surface roughness). It may also vary with the slope steepness factor and surface roughness. For example, areas with a lower slope and higher vegetation and litter layers tend to have lower erosivity (Farhan & Nawaiseh ). In this study, the rainfall data of the nearest weather stations to the watershed were used to determine the linear relationships between average annual rainfall (Table 1) and computed I 30 values using DEM. In the RUSLE model, the yearly precipitation surface was interpolated to determine the value of each cell based on the values of nearby cells (Sharma & Goyal ). This average precipitation data was interpolated in ArcGIS using an ordinary Kriging tool for the clipped watershed of Deurali VDC. The interpolation of average annual precipitation data was applied to obtain a representative rainfall distribution map and subsequently the R factor map. The interrelationship between annual precipitation and the R factor was determined from Equation (1) below (Renard & Freimund ): where P is annual precipitation (mm).

Slope length and steepness factor (LS)
Ganasri & Ramesh () defined L as the distance from the point of origin of overland flow to the point along the slope where either deposition starts or where runoff becomes concentrated in a defined channel. S represents the slope gradient. The LS factor determines total sediment yield from a site; yield being more sensitive to steepness than slope length. In general, an increased slope length and steepness will translate into increased soil loss due to higher runoff velocity and erosivity. Rill erosion generally increases in a downslope direction, whereas inter-rill erosion is reported as a rather uniform process along a slope. The L factor also connotes the ratio between these rill and interrill erosions (Farhan & Nawaiseh  Briefly, (a) fill was calculated from clipped watershed DEM layer using a hydrology tool, (b) flow direction was calculated from a flow direction tool using fill data as the input raster, (c) flow accumulation was calculated from a flow accumulation tool using flow direction data as the input raster, and (d) slope of watershed was calculated in degrees from a slope tool using a clipped watershed DEM as the input layer.

Crop management factor (C)
Surface covers such as vegetation, residues, and mulches are used for crop management to intercept raindrop impacts, lower erosion risks, and delay time to runoff start. The C factor of 0 represents a well-protected soil while the C factor of 1 is used to define the baseline condition with a where NIR (near infrared) is the reflection of the near infra- where A ¼ average annual soil loss rate (t ha -1 yr -1 ); R ¼ rainfall factor (MJ mm ha -1 hr -1 yr -1 ); K ¼ soil erodibility (t hr ha ha -1 MJ -1 mm -1 ); L ¼ slope length factor; S ¼ slope steepness factor; C ¼ crop management factor; P ¼ support practice factor.

RESULTS AND DISCUSSION
This study used the RUSLE modeling approach to estimate the spatial distribution of soil erosion in the Girwari river watershed. Indeed, the amalgamation of RUSLE model with GIS presents an effective tool to predict potential soil

Rainfall erosivity factor (R)
The highest precipitation (2,074 mm) was obtained along the southeast part while the lowest value (2,058.52 mm) was recorded in the northwest part of the study area

Soil erodibility factor (K )
The K value varied between 0.066 and 0.087 t ha. hr ha -1 MJ -1 mm -1 with the highest being in the hill and the lowest in the plain areas of the study site ( Figure 6).

Crop management factor (C )
The lowest C value (0.051) was obtained in areas adjoining the river and in the settlement areas of the plain land ( Figure 9). The C factor was predicted to be higher in the bare areas of hill and Terai region. Even within the forests, areas with degraded vegetation showed higher C values, indicating low crop management compared with the dense forest lands (Figure 10).

Conservation practice factor (P)
The P value varied from 0.5 to 1.0 ( Figure 11). The highest P  contouring, and check-dams were adopted; for example, in the cultivated lands of Terai. Alternatively, a higher P value was observed in areas with no use of conservation measures to regulate run-off flow.
These five factors of the RUSLE model were multiplied in a raster calculator tool to estimate soil erosion rate. The soil erosion rate varied from less than 2 to 100-190 t ha -1 yr -1 as discussed separately below ( Figure 12).
Soil erosion was also determined for six land use and land cover types within the watershed. The total erosion from the total hill forest area (2,147.25 ha) was 13,374.3 t yr -1 , from the total hill cultivated land (672 ha) was 2,446.74 t yr -1 , from the total Terai cultivated land (681.20 ha) was 89.132 t yr -1 , from the total Terai forest area (502.44 ha) was 37.539 t yr -1 , and from the total area of grassland (24.44 ha) of Terai was 2.852 t yr -1 (Table 4).
In simple words, nearly 85% erosion occurred in hill  The results identified nearly 13% of the total land as being highly vulnerable to soil erosion, and a higher landslide frequency ratio in less dense forest area, cropping areas, streams, and urban lands. differ owing to soil types, depth, parent material, land-use, and/or climatic conditions. For example, T values ranging from 2.5 to 12.5 t ha À1 yr À1 were reported in India, due largely to variable soil depths and quality (Mandal & Sharda ). In mountainous ecosystems such as in young mountains of Nepal, soil loss of up to 25 t ha À1 yr À1 is deemed tolerable (Morgan ; Shrestha ; Koirala et al. ).
A default T value of 15 t ha À1 yr À1 is, however, being used in Nepal. Therefore, areas exceeding 15 t ha -1 of soil loss annually were identified and classified as the erosion hotspots in this study watershed to determine soil conservation priorities. Soil erosion hotspots, thus, identified included 176.52 ha (4.22%) of land whereas a 4,008.6 ha (95.78%) area of the Girwari watershed was below the threshold limits (Table 5)     In summary, the use of the RUSLE model and GIS provides a promising tool to map the spatial distribution of soil erosion risks, and results serve a valuable reference guide to other studies and policy makers for effective soil and water resource planning and management. However, RUSLE excludes the effect of mega-rill, gully, bank and channel erosion, and landslides, and deposition of eroded soil, and as such, priority should be given to field measurements of rainfall, soil properties, slopes, and NDVI for inclusion in erosion factors to improve model prediction and accuracy.