Assessment of Land Cover Change and Its Impact on Changes in Soil Erosion Risk in Nepal

Abstract

Land cover change is a critical driver for enhancing the soil erosion risk in Nepal. Loss of the topsoil has a direct and indirect effect on human life and livelihoods. The present study provides an assessment of the decadal land use and land cover (LULC) change and consequent changes in the distribution of soil erosion risk for the years, 1990, 2000, and 2010, for the entire country of Nepal. The study attempted to understand how different land cover types change over the three decades and how it has changed the distribution of soil erosion risks in Nepal that would help in the development of soil conservation priority. The land cover maps were produced using geographic object-based image analysis (GEOBIA) using Landsat images. Soil erosion patterns were assessed using the revised universal soil loss equation (RUSLE) with the land cover as the input. The study shows that the forest cover is the most dominant land cover in Nepal that comprises about 6,200,000 ha forest cover. The estimated annual erosion was 129.30 million tons in 1990 and 110.53 million tons in 2010. The assessment of soil erosion dynamics was presented at the national, provincial, and district level. District wise analysis revealed that Gulmi, Parbat, Syangja, and the Tanahu district require priority for soil conservation.

1. Introduction

Land degradation is identified as one of the major environmental challenges that affects all landforms of the fragile Earth’s skin. Soil erosion is one of the key drivers of land degradation in mountain regions due to the loss of topsoil, thereby reducing soil fertility and consequently decreasing agricultural productivity and influencing habitat destruction [1,2]. Soil erosion causes an alteration in the water retention property of soil. Eroded soils are carried by the flow of water and deposited in the river channel that intercepts the water flow [3,4,5]. The debris flows due to erosion from upper land deposited in the river that reduces the carrying capacity, causing siltation in the irrigation canals, and damage to water control structures and turbines [6]. Erosion has been accelerating due to land use and land cover change and unsustainable management practices of land use and land cover change [2,7]. Countries, like Nepal, which have complex and diverse topography and improper land use management practices and has water-induced soil erosion, which is eventually transported to the Ganges delta, a downstream river basin forming islands in the Bay of Bengal [8]. Different case studies have revealed that land use land and land cover changes (LULCC) influence soil erosion at various spatiotemporal scales [2,7,9]. The area with high erodible soil is also prone to a higher risk of landslides [10,11]. Therefore, it is crucial to understand LULCC and the potential impacts it has on soil erosion dynamics [12]. Recognizing the trend and distribution of LULCC and soil erosion risk would be useful for better management of land cover and land use and reduce the risk of disaster from landslides.

Several models have been developed to assess soil erosion driven by water and associated sediment yield, which are categorized into empirical, regression, and conceptual models. Empirical or regression models include the universal soil loss equation (USLE) [13], and the revised universal soil loss equation (RUSLE) [14]. Among the various models, the RUSLE model is widely used to evaluate and identify the soil erosion due to its modest data requirement and simple, comprehensible model structure. The model is widely used for countries, where a lack of sufficient input data limits the implementation of complex models [15]. The model is useful to identify the spatial distribution of soil loss for a larger region [3] as it can estimate erosion risk on a regular grid [16]. Assessing soil erosion through the integration of soil erosion models, field data, and data obtained from remote sensing platforms using geographic information systems (GIS) methods has shown satisfactory results [3,7,17]. The integration has made the assessment and prediction easier and more efficient to address the issue of soil erosion [9].

Few studies on the estimation of soil erosion at a small scale have been conducted for different parts of Nepal. Optical images of high and medium resolutions with the GIS technique were used to examine soil erosion-prone areas in the Mustang watershed [18]. The RUSLE model in GIS platform has been used to study soil erosion in the Trijuga watershed [19], the entire Koshi river basin [17], and the Karnali river basin [20] of Nepal. Another model, such as the Morgan, Morgan, and Finney model [21], in a GIS environment, has been applied to assess the Likhu Khola valley [22]. Both these models were used to map soil erosion in the Kalchi Khola watershed of Nepal [23]. However, the assessment of soil erosion risk and distribution at the national scale is still lacking [24,25,26,27]. The overall purpose of this study was to identify the spatial patterns and temporal trend of soil erosion in Nepal and to determine how changes in land cover types affect the distribution of soil erosion risks. We have used land cover maps of 1990, 2000, and 2010 to estimate soil erosion risk represented by soil erodibilty for those years to understand how different land cover types affects the rate of soil erosion, and developed an erosion risk map for Nepal that would be useful for soil conservation prioritisation.

2. Materials and Methods

2.1. Study Area

The study area covers all of Nepal (Figure 1), which falls between latitudes, 26°22′ N to 30°27′ N, and longitudes, 80°04′ E to 88°12′ E, sharing borders with China to the north and India to the west, south, and east, with a total land area of 147,181 km². It is divided into five physiographic regions: High Mountain (34,475 km²), Mountain (28,895 km²), Hill (43,503 km²), Siwalik (18,886 km²) and Terai (21,422 km²). Administratively, Nepal has seven provinces and 77 districts. Elevations range from 60 m from the southern plains to 8848 m at Mt Everest in the north, which is the highest point on the Earth [28]. Nepal encounters subtropical to alpine climates due to the variation of topography and geographical location. Lowest precipitation occurs (160 mm) in the north-western region and highest precipitation (5500 mm) in the central region of Nepal. More than 80% of the of annual precipitation occurs from June to September [29]. According to the Nepal national population and housing census of 2011 (NPHC 2011), the population of Nepal is 26,494,504 people with a 1.35% annual population growth rate [30]. In Nepal, around 66% of the population is directly engaged in small-scale agriculture farming, which provides about 33% of the gross domestic product. Maintenance of soil resources, particularly soil erosion, is critical due to the agriculture practice of grazing on rangeland beyond the carrying capacity and topography [24].

2.2. Assessment of Land Cover

Assessment of Nepal land cover of 1990, 2000, and 2010 was done using the Landsat Surface Reflectance Level-2 data (Figure S1). Thirteen Landsat images covered the whole study area (each scene approximately 185 × 185 km). All Landsat data were downloaded via United States Geological Survey (USGS) data archived portal (https://earthexplorer.usgs.gov). For the estimation of land cover, we adopted an object-based image analysis (OBIA) method, which is widely known as a geographic object-based image analysis (GEOBIA). In recent years, the GEOBIA approach has showed potential to efficiently classify land cover information with a higher accuracy [31,32]. During the image analysis, the non-overlapping image segmentation technique was used to produce meaningful image objects based on the input satellite image information as well as other thematic layers. In the region, GEOBIA was used for the national land cover change assessment of Bhutan [33]. In the case of Nepal, the method was applied for the 2010 land cover study [28]. The same GEOBIA method was later applied to generate land cover for 1990 and 2000. All the land cover data produced through these studies were made available to the public through the International Centre for Integrated Mountain Development (ICIMOD) (http://rds.icimod.org).

The downloaded Landsat images were uploaded in the eCognition Developer 8.1 software for classification (Figure 2). For the segmentation of images, the “multiresolution segmentation” algorithm was used, which consecutively clusters pixels into image objects based on their neighbours, and relative homogeneity criteria [34,35]. During the multiresolution image segmentation in the eCognition developer, image layer weights were used as 1 for the TM and ETM+ band; only for band 6, the weight was used as 0, scale parameter used as 16, Shape as 0.1, the composition of homogeneity creation shape used as 0.1, and compactness used as 0.5. Multiresolution segmentations of image objects assumed that similar features would have a similar kind of spectral responses and that the spectral response of an element would be unique concerning all other aspects of interest. In the multiresolution segmentation process, the homogeneous areas used to combine into larger image objects and the heterogeneous regions into the smaller object or segment [36].

During the classification, information of image layers’ spectral values generated image indices, like the normalized difference snow index (NDSI), normalized difference vegetation index (NDVI), normalized difference water index (NDWI), land water mask (LWM), and bare soil index (BSI), which were created through band ratioing. Image indices were critical to the object based image classification. Image ratioing is a “synthetic image layer” created from the existing bands of a multispectral image [37,38]. The next step is to set rules for these image objects according to their attributes, such as the NDSI, NDVI, NDWI, LWM, and BSI layer value. In this process, the image objects that represent patterns are determined with the help of ground truth information and comparison with high resolution images from Google Earth. Rule sets were developed to generate a land cover map, exploring image band and indices’ mean value by 2d features plot in the ground truth segment. After computing the class separation using the segmentation approach, classification was performed to obtain a land cover map for each scene. Due to the data analysis performance and satellite image resolution, the prepared land use and land cover map and its change might differ significantly from the ground reality. To ensure the confidence and reliability, statistically robust and transparent approaches are critical to ensuring the integrity of land change information [39,40]. After the preparation of land cover maps, the results were validated using available field data and Google Earth images of 2000 and 2010 over randomly selected points. Once the image classification for each scene was finalized, all the scenes were mosaicked to obtain the land cover map of Nepal. After compiling a land cover map for the entire country, the land cover changes were assessed and the individual land cover maps were used for deriving the support practice factors for the erosion risk estimation.

2.3. Estimation of Soil Erosion Areas

The models applied to estimate the soil erosion driven by water and associated sediment yield are categorized into empirical, regression, or and conceptual models. RUSLE is useful to identify the spatial distribution of soil loss for a larger region [3]. Assessing soil erosion through the integration of soil erosion models, field data, and data obtained from remote sensing platforms using geographic information systems (GIS) methods have shown satisfactory results [3,7,17]. The well-known RUSLE equation is presented as [13]:

$$\mathrm{A}=\mathrm{R}\times \mathrm{K}\times \mathrm{L}\times \mathrm{S}\times \mathrm{C}\times \mathrm{P}$$

where, A is the soil loss in tons per hectare, R is the rainfall-erosivity index, K represents the soil erodibility index, L represents the slope length, S is the slope steepness factor, C represents the land cover management factor, and P represents the supporting practices factor. The description of input data, data preprocessing, and calculation of each factor are discussed in greater detail below. These factors (RKLSCP) were combined using the RUSLE framework in the ArcGIS model builder environment for the estimation of average soil loss. The overall processing flow is shown in Figure 3.

Development of Specific RUSLE Data Layers

The equations used for the calculation of individual factors are given in

Table 1. Input data sources, and published equations used in RUSLE.

Soil Erosion Factor	Input Data	Data Source	Equation or Published Factor
Rainfall erosivity factor (R)	Precipitation (ESRI grids, 10 arc-minutes)	World climate precipitation data [41]	R = 0.0483 × P1.610 where P = annual precipitation (mm)
Soil erodibility factor (K)	Soil maps from Nepal		Literature review
Slope length factor (L)	SRTM 30 m digital elevation data	[42]	L = (λ/22.13) m where λ is the field slope length (m), and m assumes a value from 0.2 to 0.5 [43]
Slope steepness factor (S)	SRTM 30m digital elevation data	[42]	S = (0.43 + 0.30 s + 0.043 s2)/6.613 [13]
Land cover management factor (C)	NDVI from Landsat	[44]	C = 0.431 − 0.805 x NDVI [45]
Support practice factor (P)	Land cover map	ICIMOD ([28]	Literature review

. The rainfall erosibity factor, slope length factors, and soil erodibility were estimated using several site-specific equations.

Rainfall Erosivity Factor (R)

The rainfall erosivity factor (R) used to be an important layer for soil loss assessment under future land use and climate change [17,46]. The rainfall erosivity factor shows the ability of rainfall to cause erosion on the earth. Normally, raindrops with higher intensity control the released amount of soil and runoff as sediment. There were not enough records of rainfall intensity in Nepal, and the records of monthly rainfall data were used for the determination of the R-factor using the average annual value. The R-factor was evaluated from world climate precipitation data [47] using the following equation [48]:

$$\mathrm{R}=0.0483\times {\mathrm{P}}^{1.610}$$

where, P = annual precipitation (mm).

Soil Erodibility Factor (K)

One of the mandatory layers for soil erosion study is the soil erodibility factor (K), which is associated with soil types that provide information on the ability of soil release or soil loss in different environmental situations, topographies, and geographic locations. The K factor also measures the susceptibility of a soil surface to erosion due to precipitation and runoff. The K factor is related to crucial soil factors, such as the soil texture, structure, organic matter, and permeability triggered soil loss [49]. The soil erodibility factor (K) is a lumped parameter and quantitative description that represents an integrated annual value of a particular soil. For a particular soil association, the soil erodibility factor is the rate of erosion per unit of the erosion index from a standard plot. In this study, K values are obtained from the published work conducted in and around of Nepal [50,51,52,53,54]. If there was two different researchers who had produced two different K values for the same soil type, we used the mean value of the two values (Figure 4).

Slope-Length Factor (L)

The slope-length factor (L) is defined as the distance from the point of origin of overland flow to another point on the site erosion potential. The slope-length factor is the ratio of field soil loss to the corresponding soil loss of 22.13 m slope length and is estimated as

$$\mathrm{L}={\left(\mathsf{\lambda }/22.13\right)}^{\mathrm{m}}$$

where $\lambda$ is the slope length, and m assumes the value of 0.2 to 0.5. [13] have identified varying values of exponent, m, for different slopes and these are as follows [43,53]. For this purpose, the m map was created, using the slope map as the input (

Table 2. Slope length exponent.

Slope Gradient	Value of m
1%	0.2
1–3%	0.3
3–4.5%	0.4
4.5% or more	0.5

).

The L value was estimated by using the field slope length calculated from SRTM DEM [42] and is given as follows:

$$\mathrm{L}={\left[90/22.13\right]}^{\mathrm{m}}$$

Slope-Steepness Factor (S)

Steep slopes and undulating topography provides a critical platform to springing precipitation water into a lower place. The slope-steepness factor (S) represents how fast water can flow on a particular surface interacting with the angle of the ground and affecting the magnitude of soil erosion. The slope factors differ according to the shape and size of the ground; at the same time, soil loss increases more with slope steepness than it does with slope length. It is the ratio of soil loss from the field gradient to that from a 9% slope under otherwise identical conditions. The relation of soil loss to the gradient is influenced by the health and density of vegetative cover and soil particle size. For estimating the S factor [13], the following equation was used:

$$\mathrm{S}=\left(0.43+0.30 \mathrm{s}+0.043 {\mathrm{s}}^2\right)/6.613$$

where s is the slope in percentage.

Cover Management Factor (C)

The vegetation cover is the second most crucial factor, and is close to topography that controls soil erosion risk and impacts on erosion [55]. When rainfall drops on the vegetated areas, the vegetation canopy intercepts the rainfall, increases the infiltration, and reduces the rainfall energy. The C factor is usually assigned based on the type and density of vegetation cover. Such an approach results in transforming land cover into discrete weighted data.

Alternatively, De Jong [56] developed the following relationship between field calibrated C factor values with a Landsat satellite-based normalized difference vegetation index (NDVI) to produce a continuous C factor surface:

$$\mathrm{C}=0.431-0.805\times \mathrm{NDVI}$$

where $\mathrm{NDVI}=\frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red}}$.

Annual composite NDVI of 1990, 2000, and 2010 were created from Landsat data in a Google Earth Engine (GEE) environment. Spatial distribution of the cover management factor for 1990, 2000, and 2010 can be seen from Figure 5.

Support Practice Factor (P)

The erosion control practice factor (P) reflects the impact of support practices on the average annual erosion rate, and it also describes the ratio between soil losses on a particular field where the erosion control practice is determined. It is measured as the ratio of soil erosion with a specific support practice to the corresponding loss with ploughing up and the downslope [45,57]. Briefly, the P factor map was prepared utilizing the 1990, 2000, and 2010 land cover map of Nepal (Figure 6), where the erodibility factor for the support practice was assigned based on the mean value of secondary published data for a particular land cover class in the region [48,50,58,59,60]. Within the study area, support practice factor values were between 0.003 and 0.906.

2.4. Determination of Soil Priority Areas

Soil erosion maps for 1990, 2000, and 2010 and trends in the soil erosion maps were generated to support the determination of soil erosion priority areas. During this analysis, we combined and calculated the cross-tabulated soil erosion rate between 1990 to 2000 and 2000 to 2010. The estimated soil erosion rate, and changes in the soil erosion rate were calculated to identify the priority areas for soil erosion management. For the identification of priority areas for controlling erosion, the map of the estimated erosion was divided into eight erosion classes: very low, low, low to medium, medium, medium to high, high, very high, and extremely high erosion severity.

3. Results

Classification results of the preprocessed Landsat images for 1990, 2000, and 2010 are presented in Figure 6 and the statistics for the three different years of land cover with estimated soil erosion rates are shown in

Table 3. Land cover and estimated erosion rates in Nepal in 1990, 2000, and 2010.

Land Cover	Land Cover Area (ha)			Annual Soil Loss (t)			Mean Erosion Rate (t/ha/yr)
Year	1990	2000	2010	1990	2000	2010	1990	2000	2010
Forest	6,668,336	6,148,401	6,202,809	9,557,443	4,968,762	5,488,194	1.43	0.81	0.88
Shrubland	328,142	346,930	342,986	1,223,757	849,352	1,449,670	3.73	2.45	4.23
Grassland	1,728,561	1,379,485	1,264,552	25,198,391	13,941,821	13,982,884	14.58	10.11	11.06
Agriculture area	3,753,933	4,096,968	4,039,820	69,054,803	36,115,163	60,071,672	18.40	8.82	14.87
Barren area	1,006,831	1,702,002	1,535,851	21,299,992	38,469,816	23,266,676	21.16	22.60	15.15
Water body	81,052	73,051	72,685	2,029,087	1,608,671	1,597,417	25.03	22.02	21.98
Snow/glacier	1,168,741	974,176	1,255,347	919,372	682,081	4,537,558	0.79	0.70	3.61
Built-up area	32,916	47,499	54,462	33,952	39,974	152,439	1.03	0.84	2.80
Total	14,768,512	14,768,512	14,768,512	129,316,796	96,675,641	110,546,509	8.76	6.55	7.49

. The results show that forest was the most significant land cover in 1990 with 6,668,336 ha, 6,148,401 ha in 2000, and, in 2010, 6,202,809 ha of the total area of Nepal, followed by agriculture area in 1990 with 3,753,933 ha, 4,096,968 ha in 2000, and 4,039,820 ha 2010. Estimated shrubland cover in 1990 was about 328,142 ha, about 346,930 ha in 2000, and about 342,986 ha in 2010. The results show that there are variations in aspect and slope controlling the vegetation growth in the mountain zone of Nepal [37]. However, altitude and its influence on climatic conditions have much influence in the land cover distribution. Forest cover is mostly distributed in the Terai, Siwalik, and Mahabharata mountain ranges. Agriculture area is primarily distributed in the Terai area, and in and around the Kathmandu Valley, Pokhara Valley, Gandaki River Basin, Koshi River Basin, and Carl River Basin. Shrubland and grassland are mainly found in the high Himalayan region.

From Table 3, it can be concluded that land cover change was more noticeable between 1990 and 2010. The study shows that between 1990 and 2000, forest area declined, but between 2000 and 2010, forest area increased. Considering the total area of Nepal, annual net forest loss was 0.035% between 1990 and 2000, whereas between 2000 and 2010, annual net forest gain was 0.037% (Tables S4–S6). Built-up area has less significant land cover, but steadily increased with 32,916 ha in 1990, 47,499 ha in 2000, and 54,462 ha in 2010. Between 1990 and 2000, the built-up area increased to 14,583 ha and between 2000 and 2010 it increased by 6963 ha.

The analysis by Nepal physiographic regions revealed that forest area dominates in the Hill region, followed by the Middle Mountain region (Tables S4–S6). Out of the total forested area of Nepal, in 1990, the Hill, Middle Mountain, Siwalik, Tarai, and High Mountain physiographic regions contained 38.57%, 30.08%, 21.71%, 6.75%, and 2.89% of the forest area, respectively. Whereas, in 2010, the Hill, Middle Mountain, Siwalik, Tarai, and High Mountain physiographic regions contained 37.88%, 30.43%, 22.21%, 6.65%, and 2.83% of the forest area, respectively. Similarly, out of the total agriculture area, the Hill physiographic region covered the highest proportion (40.79% in 1990 and 42.35% in 2010), whereas the Eastern High Mountain region covers the lowest (0.11% in 1990 and 0.07% in 2010). Physiographic region-wise analysis (Figure S2) showed that the forest area had decreased in the High Mountain region by about 0.10% between 1990 and 2000. On the otherhand, between 2000 and 2010, the forest area increased by about 0.04%. Similarly, in the Hill region, the agriculture area decreased by about 2.72% of the total forest area between 1990 and 2000.

The overall accuracy of the 2000 land cover map, assessed using the 450 reference points, was 86.67%, with a kappa value of 0.82, standard error kappa of 0.0211, 95% confidence interval of 0.78–0.86, and a 0.93 maximum possible un-weighted kappa, given the observed marginal frequencies (Table S1). Similarly, the overall accuracy of the 2010 land cover map, assessed using the 750 reference points, was 85.87%, with a kappa value of 0.80, standard error kappa of 0.018, 95% confidence interval of 0.77–0.84, and a 0.92 maximum possible un-weighted kappa, given the observed marginal frequencies (Table S2).

The years of 1990, 2000, and 2010’s soil erosion risk maps of Nepal were used to identify priority areas for soil erosion control and can be seen in Figure 7. The mean annual soil erosion for Nepal was estimated at 8.76, 6.55, and 7.49 ton/ha/yr accordingly for the year of 1990, 2000, and 2010, respectively. The identified annual erosion was 129.32 million tons in 1990, and 96.68 million tons in 2000 and 110.55 tons in 2010 million tons. The southern and northern parts of Nepal were found to be less soil erodible areas, while the middle part of the study area was found to be highly erodible. The annual soil erosion risk map of Nepal shows the spatial pattern between the different elevation zones.

The very low erosion class was found mainly in the Terai area of Nepal due to the relatively flat terrain. The differences in erosion levels in the middle elevations and the northern and southern parts of the study area are due to the topography. The maximum per hectare average soil loss was estimated in the 4000 m to 5000 m elevation zone, and the minimum was estimated below the 100 m elevation zone. Table S3 shows soil loss according to elevation zones. Table 3 presents the distribution of erosion according to land cover and land use, revealing that agriculture contributes around 50% of the total erosion, followed by barren areas with 20% of the study area. In 1990 and 2010, agriculture land contributed 69,054,803 and 60,071,672 tons to soil loss, respectively.

The erosion class under high to medium, high, very high, and extremely high represent erosion levels of 5–50 ton/ha/yr.

Table 4. Change matrix for soil erosion classes from 1990 to 2000 (%).

Soil Erosion Risk Rank (t/ha/yr)	Very Low (<0.5)	Low (0.5–1)	Low Medium (1–2)	Medium (2–5)	High Medium (5–10)	High (10–20)	Very High (20–50)	Extremely High (>50)	Total 2000
Very low (<0.5)	29.56	5.62	2.59	0.63	0.80	0.60	0.39	0.08	40.27
Low (0.5–1)	1.44	2.68	3.29	0.37	0.19	0.11	0.08	0.02	8.18
Low medium (1–2)	1.50	0.92	6.01	4.11	1.73	0.55	0.32	0.06	15.20
Medium (2–5)	1.03	0.21	0.80	1.28	2.17	1.04	0.42	0.06	7.00
High medium (5–10)	1.38	0.40	0.75	0.50	2.57	3.02	1.79	0.19	10.61
High (10–20)	0.90	0.35	0.42	0.39	0.74	2.64	3.54	0.64	9.63
Very high (20–50)	0.59	0.09	0.32	0.08	0.42	0.95	2.91	1.81	7.16
Extremely high (>50)	0.19	0.01	0.05	0.03	0.02	0.15	0.52	0.98	1.95
Total 1990	36.57	10.27	14.24	7.40	8.66	9.06	9.97	3.84	100

and

Table 5. Change matrix for soil erosion classes from 2000 to 2010 (%).

Soil Erosion Risk Rank (t/ha/yr)	Very Low (<0.5)	Low (0.5–1)	Low Medium (1–2)	Medium (2–5)	High Medium (5–10)	High (10–20)	Very High (20–50)	Extremely High (>50)	Total 2010
Very low (<0.5)	32.44	2.23	1.64	1.08	1.55	1.02	0.63	0.19	40.79
Low (0.5–1)	3.87	3.16	1.19	0.16	0.31	0.37	0.16	0.01	9.23
Low medium (1–2)	1.49	2.03	5.11	0.61	0.68	0.33	0.26	0.06	10.57
Medium (2–5)	0.68	0.34	4.34	1.30	0.51	0.46	0.08	0.02	7.72
High medium (5–10)	0.83	0.22	2.08	2.60	2.97	0.91	0.56	0.02	10.20
High (10–20)	0.57	0.11	0.48	0.89	3.30	3.23	1.16	0.22	9.94
Very high (20–50)	0.32	0.07	0.31	0.32	1.17	3.02	3.24	0.55	9.00
Extremely high (>50)	0.06	0.01	0.05	0.05	0.12	0.30	1.09	0.87	2.55
Total 2010	40.27	8.18	15.20	7.00	10.61	9.63	7.16	1.95	100

show the erosion change matrix between 1990, 2000, and 2010. To control erosion in these high, very high, and extremely high erosion zones, high priority is required to be given to soil erosion conservation. Identified priority areas of Nepal for soil erosion control can be found in Figure 8. Based on the risk, the highest priority districts for soil conversation are identified as Gulmi, Parbat, Syangja, and Tanahu and the second highest priority districts are identified as Nuwakot and Ramechhap. The third highest priority districts are Dolpa and Mustang (Table S7). Conservation measures, such as creating plant communities that are environmentally favourable, bioengineering during construction of infrastructure, cultivation on steep slopes, and controlling overgrazing, could help to decrease the erosion severity in this zone.

4. Discussion

Considering land cover as a critical factor for the environmental and erosion study of Nepal, the present study provides findings and its analysis of land cover changes from 1990 to 2000 to 2010 using Landsat images of a 30m spatial resolution with a consistent methodology and standardized legend. The main worth of Nepal land cover change lies in its ability to provide complete, coherent, and harmonized national land cover maps. Land cover maps serve as a resource for regional to national scale applications. Although, national level land cover assessments were done during 1963 and 1988 as part of a forest resources’ survey using aerial photography and land resource mapping project using satellite data [61,62,63]. In the past, most of the study has engaged in the manual digitization process. There are various merits and demerits for thematic maps’ preparation using digitization. Most of the cases of manual digitization introduce error into the database, and the amount of error is often left unmeasured [64,65]. In this study, we applied state-of-the-art image segmentation, an emerging paradigm of land cover mapping, using object-based image analysis to produce high-quality land cover data for Nepal [66,67]. Recently, Uddin, Shrestha, Murthy, Bajracharya, Shrestha, Gilani, Pradhan, and Dangol [28] have developed the 2010 national land cover database with 12 land cover classes for Nepal using Landsat images. Similar year land cover maps for 2010/11 were prepared using RapidEye satellite images for forest resources’ assessment [63,68]. Since then, there have been no national-level decadal land cover change assessments carried out for Nepal. Although national land cover data is very much essential for the national reporting to the United Nations Framework Convention on Climate Change (UNFCCC) for forest reference level (FRL), assessing forest carbon fluxes as per the good practice guidelines of IPCC, developing national biodiversity conservation strategies, and evolving national level REDD mechanisms and sustainable land management [28,69] is also important.

With the limitation of the availability of field-based measurements, we have attempted to adopt the RUSLE model to bring out national level erosion estimates of Nepal [70]. The estimated erosion rate for different land cover classes was evaluated with published field based erosion measurements and model-based results mostly in the mid and high hills of Nepal (Tables S8–S10). Due to deep shadows of high mountain ranges, land use dynamics, and a high degree of spatial variability of rainfall, a lot of earlier studies focused significantly on erosion assessment concerning the middle mountains of Nepal’s Himalayas [50,51,53]. Studies addressing national level erosion assessment and understanding its spatial trend and its relation to land use practices and rainfall regimes are lacking. Also, there is no consistent protocol-based field measurements that exist over the country to provide an understanding of the erosion levels and associated factors. It can be found that the estimates were within the ranges of the published studies [23,71,72]. However, it should be mentioned here that there is heterogeneity in rainfall, topography, soil and cultural practices, causing the high range of variation in erosion levels in the assessment at different sites. Therefore, one to one comparison of estimates over a set of sites is essential. However, it was not possible make a comparison over a set of sites to bring out error and bias of the estimates as very few studies were reported in Nepal. Given this, it could be mentioned that the estimates are semi-quantitative and reliable on relative terms. These spatial estimates would be useful to assess the spatial variability of erosion and prioritizing areas of conservation.

Land cover is the physical materials of the earth surface, which anchors all life on Earth. Many studies have revealed that the soil erosion rate is triggered by change in land-use and land cover [7,73]. Deforestation in natural forested areas, which anchor the soil with their roots, is a direct cause of erosion. Increased demand for agriculture commodities generates incentives to convert forests and grasslands to cropland. The transition to farmland from healthy natural vegetation often reduce the holding capacity of the soil. In our analysis, the mean soil erosion rate was decreased from 8.76 t/ha/yr in 1990 to 7.49 t/ha/yr in 2010. The main reason for the changes in soil erosion between the land cover types is due to the changes in the land cover, which were reflected in the greenness of NDVI, and changing the C factor values, which affected the soil erosion strongly. The results could be interpreted as a change in land use practice due to the internal migration [74,75]. The main reason for this decline is the progressive decrease of the agricultural areas; a large percent of the cultivated area changed to uncultivated areas in the study area as well as the community forestry [76,77].

The rainfall and slope length factor (LS) and soil were the major drivers of erosion in the RUSLE equation. The rainfall erosion potential is mainly controlled through the product of total storm energy and maximum 30-min storm intensity. In the absence of detailed rainfall data at the sub-hourly interval and also the lack of Nepal specific equations on the rainfall erosion factor, we adopted available annual rainfall based equations that are suitable for hilly areas. Concerning quantification of the slope length factor, there are several equations to estimate the slope length factor in the GIS domain using the digital elevation model. Most of these equations in our study found a yielding overestimation of erosion and the best among them, which provides estimates in tune with the published literature, was chosen. Hence, the LS factor is the most important variable, which requires calibration over the study area to address reliable quantitative estimates.

The soil erodibility factor was weighted at the soil order level using published literature [3,17,22,23]. In the case of the availability of soil texture and organic carbon information, the erodability factor can be better estimated. The crop management factor and support practice factors are found to be more site specific and broadly match with the generic properties of the sites. Any improvements to be made on weightage to be given for these parameters require intensive local specific ground observations. RMF and MF models are said to yield better estimates over hilly terrains. These models require extensive ground driven data and calibration. Therefore, a holistic discussion on the level of accuracy requirements of erosion levels should be made to plan for appropriate model and ground measurements.

The study disclosed different estimations of land cover and erosion statistics. In the quest for producing better quality land cover and soil erosion data sets, it is essential that the study focused on putting feet’s utilization of the modern image analysis technique with a sufficient training sample and comparing this with high-resolution Google earth images. In the mountain area, the yawning shadow is a challenge to extracting land cover information, although several tools are available for shadow removal. For a national land cover reporting time, country boundaries, and its projection, raster and vector file overlapping issues might produce little different figures.

5. Conclusions

This study provided baseline information as well as information about the changes in land cover over the past two decades (1990, 2000, and 2010). Some applications and analyses can be executed on these three land cover layers for environmental management in Nepal. The land cover map might support specific national and international reporting, such as greenhouse gas inventories and environmental-economic accounting. In addition, an easily adaptable RUSLE based method was developed for assessing erosion risk in Nepal. Indeed, a qualitative assessment method is the best alternative to field-based, time consuming soil erosion assessments. Assessed spatial erosion zones comprise the necessary information for the planning of successful and sustainable management practices. The analyzed land cover and soil loss dynamics may support decision-makers and planners to take relevant forest and soil conservation priority steps and thereby reduce the soil loss and land degradation issue of Nepal.