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

Often defined as a geomorphological and land degradation process, soil erosion presents a major threat to natural and managed ecosystems (Montgomery 2007; Montanarella et al. 2016; Poesen 2018). Drivers and impacts of soil erosion are, however, complex and multifaceted. Drivers may involve overexploitation of natural resources, land mismanagement, weak governance, population growth, and climate change, among others (Boardman et al. 2003; Montanarella et al. 2016; Chalise et al. 2019a, 2019b; Wynants et al. 2019).

Soil erosion affects global agricultural production and its long-term sustainability. The impacts of soil erosion vary with soil structure, vegetation cover, and land topography (Pimentel 2006). Soils with medium to fine texture, low organic matter, and poor structure have higher erosion rates due to low infiltration, reduced aggregate stability, and higher runoff and interflow. Effects may manifold through soil slaking, compaction, and crusting. Soils with limited contact cover, marginal areas, and steep lands are subjected to greater raindrop impact and shearing by wind and water (Lal 2001; Pimentel 2006). Soil erosion results in preferential loss of nutrients, soil organic matter, water availability, and soil thickness, thereby reducing soil quality, biodiversity, and agricultural productivity (Pimentel 2006).

Soil erosion is one of the gravest environmental threats to the mountainous ecosystems of Nepal causing soil degradation, loss of water quality, and aggravated crop productivity (Al-Kaisi 2000; Le Bas 2007; Nyssen et al. 2009; Chalise et al. 2019a, 2019b). For example, Chalise et al. (2018) estimated a mean annual soil erosion rate of 11.17 t ha^–1 in the Aringale Khola watershed of Nepal resulting in serious ramifications on soil function and services. Ghimire et al. (2013) estimated that 64 t ha^–1 yr^–1 of sediments is eroded within the Khajuri watershed of Siwalik Hills in Nepal. Rivers in Nepal transport nearly 336 million tons of soil and sediments down to India every year; posing a threat to regional soil and water quality (Thapa 2009). Key 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. 2019a, 2019b). 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.

A geographic information system, remote sensing, and modeling tools are increasingly used in environmental studies for land use land cover planning and water management (Nearing et al. 2017; Barakat et al. 2019). For example, Lamqadem et al. (2019) used remote sensing and vegetation field measurements to understand and evaluate land-use change and their environmental effects for the last three decades in the Feija Basin, Morocco. Several studies have utilized different empirical, conceptual, and physically based models to predict the detachment, transport and deposition of soil particles (Yu et al. 1997; Viney & Sivapalan 1999; Nearing et al. 2017; Phinzi & Ngetar 2018; Wang et al. 2020). Empirical models such as the Universal Soil Loss Equation (USLE), RUSLE, Chemical Runoff and Erosion from Agricultural Management systems (CREAMS), Agricultural Nonpoint Source Model (AGNPS) and Modified Universal Soil Loss Equation (MUSLE) are increasingly being used to predict erosion rates at spatio-temporal scales (Prasannakumar et al. 2012; Uddin et al. 2016, 2018; Chalise et al. 2018; Koirala et al. 2019). The USLE model was originally proposed in 1965 to estimate rill and sheet erosion from cultivated lands in the US (Wischmeier & Smith 1978). However, the model is continuously being revised and used in different countries and ecosystems.

RUSLE is a simple and robust model, and has significant advantages in watershed scale erosion and sediment transport modeling (Merritt et al. 2003). RUSLE offers several improvements in the USLE factors to represent diverse land and crop management, and slope forms (Renard et al. 1991). 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, strip-cropping, and rangeland conservation and management practices (Renard et al. 1991, 1994). In a recent study, the RUSLE model was used to study the spatial distribution of soil erosion in Nepal (Koirala et al. 2019), which estimated a nation-wide mean annual soil loss of 25 t ha⁻¹ yr⁻¹. The mean soil erosion rate under different land use occurred in the following order: barren land (40 t ha⁻¹ yr⁻¹) > agricultural land (29 t ha⁻¹ yr⁻¹) > shrubland (25 t ha⁻¹ yr⁻¹) > grassland (23 t ha⁻¹ yr⁻¹) > forests (22 t ha⁻¹ yr⁻¹). Also, Uddin et al. (2018) used the RUSLE-GIS interface to estimate erodibility and to map erosion risks in Nepal under land cover change over the past two decades. RUSLE has been applied to a few small watersheds in Nepal to predict soil erosion and to develop soil management plans. To cite examples, Jha & Paudel (2010) used RUSLE and a revised Morgan, Morgan and Finney model to predict spatial erosion pattern and rates in the Kalchi Khola watershed of Nepal. Chalise et al. (2018) used GIS and RUSLE to model soil erosion rates in the Aringale Khola watershed of the middle hills in Nepal. 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.

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². The country is ecologically divided into three regions: Terai (lowland), hills, and mountains. The physiographic characteristics of different ecoregions are available from Paudel et al. (2014). The Girwari watershed boundary of Deurali VDC in eastern Nawalparasi lies between 27°44′04″N and 27°38′28″N latitude and 84°00′27″E and 84°04′36″E longitude with a total area of 41.9 km² (4,189.15 ha) (Figure 1). Under recent administrative restructuring of Nepal, the Deurali VDC falls under the Hupsekot Rural Municipality. The Girwari river watershed is located in the Siwalik Hills, and has a subtropical climate with a maximum temperature of 43 °C, minimum of 5 °C and average precipitation of 2,493 mm (Department of Hydrology and Meteorology). The elevation ranges from 195 to 1,005.46 m above the sea level.

The Girwari river is the main river in the Deurali VDC, which flows from the uppermost part of the VDC and meets the Narayani river. The watershed encompasses different land use and land cover: (i) hill forest, (ii) Terai forest, (iii) hill cultivated land, (iv) Terai cultivated land, (v) grassland, and (vi) river. The majority of the area in the watershed is hilly with more forest trees and less agricultural land. The base of the hill is relatively flat/plain and is dominated by agricultural land, forest and residential areas. Within the Terai land, rice, maize and wheat are the principal agricultural crops. In hilly parts of the watershed, rice/millet maize cropping system is practiced. While agriculture remains a mainstay of the watershed dwellers, intensive farming and land mismanagement is increasing within the watershed. In the past, Siwalik Hills have observed several soil erosions due to both increased agricultural and non-agricultural activities (Tamrakar et al. 2002; Ghimire et al. 2006, 2013). 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 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–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 literature and through field observation. The conceptual framework for soil erosion is explained in Figure 3.

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.

The P factor is the ratio of soil loss with a specific cropland conservation practice to the corresponding loss with slope-parallel tillage (Wischmeier & Smith 1978). The P factor varies between 0 and 1. A lower P factor indicates effective and anthropic erosion control practices and facilities whereas a P factor of 1 indicates non-anthropic facilities (Koirala et al. 2019). Control practices generally reduce erosion risk by influencing drainage patterns, and runoff characteristics such as concentration, velocity, and hydraulic force against soil (Ganasri & Ramesh 2016). Values of the P factor were assigned according to the land use in vector format and converted to raster format to use with other factor layers. The P values were assigned as per Table 2.

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 erosion rates for the target watershed. Five different factors including rainfall erosivity, soil erodibility, slope length and steepness, land cover management, and soil conservation practices were calculated individually under our settings and used to estimate average annual soil loss per unit area. While the model has been widely used in different parts of the world, its use is still limited to very few watersheds in Nepal. To the best of our knowledge, this study is the first attempt to model and map erosion risks in the Girwari watershed in Nepal. Below, we present results for each factor separately.

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 (Figure 4). The R factor ranged between 10,569.41 and 10,441.65 MJ mm ha^–1 hr^–1 year^–1 with the highest value being in the southeastern and the lowest value in the northwestern parts of the watershed (Figure 5). Erosivity increased towards the eastern part of the Girwari watershed due to the relatively higher monsoon (rain). The monsoonal effect is higher in eastern Nepal relative to western Nepal; and this trend exists for our watershed as well. Morin et al. (2018) observed a strong correlation between rainfall events and sediment export in the Narayani basin, Nepal; suggesting a dominant effect/role of monsoon on the erosional process and sediment fluxes. The relative contribution of rainfall, however, varies with site and land management. For instance, a heavy storm event coupled with low plant cover and steep slope could exacerbate the erosion process and sediment transport.

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). Erodibility has been found to be higher in most of the hilly area in comparison to the Terai counterparts of the watershed, which may possibly be due to the dominance of silt and fine sand, among other factors. Sand and silt particles of the hilly area are highly susceptible to losses through surface run off; thereby generating higher K values.

The slope degree value ranged from 0 to 62°; very strong and steep slopes (≥17°) were observed in the hilly areas of the watershed (Figure 7). The LS value was lowest (0) in the plain areas and highest (1.2) in the more sloped areas across the watershed (Figure 8). Higher slope and steepness of the hilly area could possibly accentuate flow accumulation and the magnitude and potential of soil erosion in hilly areas.

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).

The P value varied from 0.5 to 1.0 (Figure 11). The highest P value (1.0) was obtained in areas covered by forest and dense vegetation in the hilly parts as well as in the Terai lands while the lowest (0.5) was obtained in grassland as well as Terai cultivated lands. The intermediate P value (0.8) was obtained in cultivated lands of hill. The P value was low where soil conservation measures such as terracing, 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).

Based on six soil erosion classes, the erosion rate was found to be <2 t ha^–1 yr^–1 in 51.62%, 2–5 t ha^–1 yr^–1 in 18.29%, 5–10 t ha^–1 yr^–1 in 19.65%, 10–50 t ha^–1 yr^–1 in 10.24%, 50–100 t ha^–1 yr^–1 in 0.14%, and 100–190 t ha^–1 yr^–1 in 0.032% of the total Girwari river watershed area (Table 3).

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 forest, 14.26% erosion in hilly cultivated land, 0.56% erosion in Terai cultivated land, 0.23% erosion Terai forest, and 0.017% erosion in grasslands within the watershed. Furthermore, the results illustrate the vulnerability of areas with very strong and steep slopes (≥17°), lower plant cover, and intensive cultivation practices to soil erosion. Similar to our findings, Koirala et al. (2019) also observed the highest erosion rate under cultivated lands compared with native grasslands and forest. In the Sabaragamuwa Province of Sri Lanka, RUSLE was used to develop erosion hazard maps using land-use change and a landslide frequency ratio and to identify sustainable land management options for erosion and sediment control (Senanayake et al. 2020). 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.

The soil erosion class was defined based on the amount of soil loss from an area over time (Morgan et al. 2004). Soil losses amounting to <2 to 10–50 t ha^–1 yr^–1 occurred in 99.8% of the watershed area. Such areas were classified as very slight to high erosion risk zones. About 50–190 t ha^–1 yr^–1 of erosion occurred in just 0.2% of the area, which were classified as severe to very severe erosion risk zones of the Girwari watershed. Less than 1% of the watershed area had a higher degree of slope and poor vegetation cover, which accounted for higher erosion rates. Similar to our study, Uddin et al. (2016) divided erosion risks in Koshi basin, Nepal, into eight classes from very low (<0.5 t ha⁻¹ yr⁻¹) to extremely high (>50 t ha⁻¹ yr⁻¹) based on erosion rates and provided recommendations for land use and land cover management at high risk areas. Recommendations included implementing conservation practices in the priority areas through a voluntary approach and multi-disciplinary partnerships among governmental and non-governmental agencies.

Soil erosion hotspots were identified as the area exceeding maximum permissible soil loss tolerance rate (T) within the watershed. They are often defined as critical sediment source areas. While the T value of 11.2 t ha⁻¹ yr⁻¹ (Wischmeier & Smith 1978) is widely documented, T values may 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⁻¹ yr⁻¹ were reported in India, due largely to variable soil depths and quality (Mandal & Sharda 2011). In mountainous ecosystems such as in young mountains of Nepal, soil loss of up to 25 t ha⁻¹ yr⁻¹ is deemed tolerable (Morgan 1986; Shrestha 1997; Koirala et al. 2019). A default T value of 15 t ha⁻¹ yr⁻¹ 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). Erosion hotspots mainly occurred in the hill forest and cultivated hilly areas where the slope was higher (≥17°). On the contrary, Terai cultivated lands and Terai forest areas showed limited to no erosion hotspots (Figure 13). The present results are also consistent with those obtained from similar ecoregions and climate. For example, Mandal (2017) observed higher erosion in steeply to very steeply sloping hilly and mountainous terrain lands in the Khadokhola watershed of Nepal, a site very similar to our study. Koirala et al. (2019) mapped the spatial distribution of soil erosion in Nepal and observed the highest mean soil erosion rate for steep slopes (34 t ha⁻¹ yr⁻¹) and the lowest for gentle slopes (3 t ha⁻¹ yr⁻¹). Areas at lower elevation with flat terrains have low erosion risks (Uddin et al. 2016). Uddin et al. (2018) observed spatial patterns in soil erosion between different physiographic regions of Nepal, and reported potentially higher erosion risks in the central part and lower risks in the lowlands of Nepal due to steep slopes and flat terrains, respectively. Indeed, sloped areas with higher slopes are highly vulnerable to erosion with hotspots likely to occur. Effects may multiply with human activities (e.g. tillage) and land mismanagement, and therefore sustainable land management and conservation interventions must be targeted realizing site-specific controlling factors. Senanayake et al. (2020) laid forth an increased emphasis on soft engineering conservation strategies, such as ecological farming, to minimize soil loss hazards. The use of conservation tillage, hedgerows, grass strips, terracing, mulch, intercrops, and cover crops will therefore improve soil cover and quality, and provide sustainable protection against erosion hazards (Acharya et al. 2019; Chalise et al. 2019a, 2019b).

Numerous studies suggest the RUSLE model as a cost-effective and efficient tool to estimate erosion at watershed scale. In Narmada river basin, India, Mondal et al. (2018) compared soil erosion estimates from RUSLE, USLE, and Morgan-Morgan-Finney models against sediment load data observed in 2009. The deviation from observed data was –4.80, 9.60, and 39.45%, respectively; thereby suggesting RUSLE as the best-fit model under their site. Also, in a Chinese Loess Plateau, Li et al. (2020) compared RUSLE predictions for soil erosion with the Pan-European Soil Erosion Risk Assessment (PESERA) model. They observed that soil loss rates from RUSLE were relatively closer to check-dam sediment yield measurements than from PESERA. While the RUSLE model has shown significant potential to explore spatial difference in erosion risks over an entire catchment, and to develop soil management and conservation plans in such several studies, many challenges and uncertainties still exist. For example, Chalise et al. (2019a, 2019b) compared RUSLE estimates against the Intensity of Erosion and Outflow (IntErO) model and erosion plots in Sarada river basin in the western hills of Nepal. The erosion estimates between the IntErO model and erosion plots were relatively closer, but the erosion rate was significantly underestimated by RUSLE. Also, Wang et al. (2019) compared the RUSLE results and the actual situation of soil erosion within the Nanling National Nature Reserve in South China and observed significant differences in the spatial distribution of soil loss, possibly due to complex topography and interactions between several factors. RUSLE excludes the effect of several geomorphological processes in its modelling such as mega-rill, gully, bank and channel erosion, and landslides, and deposition of eroded soil and sediments under transport (Chalise et al. 2019a, 2019b). Indeed, estimations of erosion risks are affected by the complex interaction of multiple factors, and as such, results from RUSLE warrants ground-truthing. However, the RUSLE model offers a huge advantage in estimating and mapping soil erosion, when extensive field measurements are challenging due to rough terrains and resource availability, such as in this study. Based on literature, we suggest that field validation of erosion estimates are needed, and thus recommend more field data accrual to increase the predictive capacity of model and reliability of outputs.

Soil loss due to erosion is a major problem worldwide. The Girwari river watershed of the Siwalik Hills has observed certain erosion risks, and utilizing GIS and RUSLE modeling could identify erosion hotspots and prioritize soil conservation strategies. This study shows the watershed areas at different severity levels: very severe (0.032%), severe (0.14%), high (10.24%), moderate (19.65%), slightly (18.29%) and very slightly (51.6%). Less than 50 t ha^–1 yr^–1 of soil loss occurred in 99.8% of the study area, suggesting no major conservation needs. However, 0.2% of the land area had 50–190 t ha^–1 yr^–1 of soil erosion, which warrants immediate attention and the use of soil conservation practices. Areas with steep slopes, low plant cover, and poor conservation measures reported elevated soil loss. Lower erosion rates were, however, reported in areas with sustainable land management practices, gentle slopes, and stable soil structure, such as Terai forest and grasslands. Overall, erosion estimates epitomize low erosion hazard and risks in the Girwari river watershed.

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.

The authors are responsible for the views expressed in this paper and do not necessarily represent or reflect the views and policies of the Oklahoma Department of Mines.

The authors extend their sincere thanks to the National Agricultural Research and Development Fund (NARDF) for providing financial support.

D. B. Tiruwa, B. R. Khanal and S. Lamichhane conducted the research and mapping, and D. B. Tiruwa, B. R. Khanal and B. S. Acharya wrote the paper.

Data cannot be made publicly available; readers should contact the corresponding author for details.