Effects of Landscape Changes on Soil Erosion in the Built Environment: Application of Geospatial-Based RUSLE Technique

Abstract

The world’s ecosystem is severely affected by the increase in the rate of soil erosion and sediment transport in the built environment and agricultural lands. Land use land cover changes (LULCC) are considered as the most significant cause of sediment transport. This study aims to estimate the effect of LULCC on soil erosion potential in the past 20 years (2000–2020) by using Revised Universal Soil Loss Equation (RUSLE) model based on Geographic Information System (GIS). Different factors were analyzed to study the effect of each factor including R factor, K factor, LS factor, and land cover factor on the erosion process. Maps generated in the study show the changes in the severity of soil loss in the Chitral district of Pakistan. It was found out that 4% of the area was under very high erosion risk in the year 2000 which increased to 8% in the year 2020. An increase in agricultural land (4%) was observed in the last 20 years which shows that human activities largely affected the study area. The outcomes of this study will help the stakeholders and regulatory decision makers to control deforestation and take other necessary actions to minimize the rate of soil erosion. Such an efficient planning will also be helpful to reduce the sedimentation in the reservoir of hydraulic dam(s) constructed on Chitral river, which drains through this watershed.

1. Introduction

World land resources are declining day by day due to soil erosion, so much that it has become the main focus of researchers and engineers [1]. Different studies have been carried out for the sustainability of natural habitat [2,3,4,5]. It is a loss of soil rich with nutrients that affects the productivity and sustainability of the original soil [6]. Approximately 80% of agricultural areas are facing higher rates of soil loss, and transported sediments severely affect the natural and built environment. Erosion of soil depletes the storage capacity of reservoirs and dams which ultimately decreases the power generation capacity. It also disturbs the agricultural output and aquatic life by polluting water in the rivers [7]. According to Chuenchum et al. [7], about 2.5 to 4 billion tons of soil is annually eroded worldwide. There are many influencing factors of soil erosion including slope, elevation, rainfall, plane curvature, drainage density, lithology, and lineaments. However, land use and climatic changes are the two most significant factors which affect the sediment transport to river tributaries. It is predicted that human actions will disturb the climate and land use land cover (LULC) significantly. Therefore, it is very important to evaluate the effect of these changes on soil erosion potential [8]. According to the past research, it has been observed that landscape characteristics are responsible for about 65% to 74% of changes in sediment yield and soil erosion. Changes in streamflow discharge due to land use also increase the intensity of soil erosion. Usually, it is observed that areas with more grassland are less vulnerable to soil erosion, whereas arable lands are more susceptible to soil loss [9].

Soil degradation is considered a serious issue due to human actions. Human beings have largely contributed to each influencing factor of land use [10]. It may be due to natural factors including different geomorphological and climatic conditions [11]. Variations in natural vegetation are observed as the first effect of this impact, irrespective of human contribution in the natural environment [12]. When LULC is changed to agricultural land from natural vegetation especially in mountainous regions, it fallouts with an increase in soil erosion rate due to crop production [13]. Land-use changes usually occur during agricultural development, which lead to the hydrological responses that further increases the rate of soil loss. Increased deforestation and growing agricultural and urban land in the tropical areas influence the intensity of soil erosion. It also plays an important role in the performance of hydraulic structures. An increased soil erosion and sedimentation fills the reservoirs of hydraulic dams, which ultimately decreases the power generation capacity.

Land use land cover change (LULCC) is the major influencing aspect that contributes to the increasing rate of soil erosion. It harms the environment by disturbing the supply of water, the capacity of storage basin, agricultural yield, and availability of freshwater in the area [8]. Land cover such as plantation directly affects soil erosion [14,15]. A generalized analysis of land use may be difficult if there is a lack of consistent time-series data for land use or land cover. Recently, scholars are interested in highlighting the effect of the environment on the erosion process. Their sole focus is on the effects of erosion including soil production and reducing the production of crops and affecting water quality by carrying nutrients, heavy metal impurities, and pesticides to surface water bodies. Sediment transport is also responsible for affecting channel, floodplain morphology, and sedimentation of reservoirs [16]. The above-mentioned impacts on soil erosion due to LULC changes can be minimized by estimating the loss of soil. For the estimation of average yearly soil loss, the RUSLE model is most commonly used by researchers, engineers, and planners. Due to its simplicity, ease of use and integration of different factors affecting soil erosion, RUSLE model is preferred over many other methods used for the estimation of soil loss [17,18]. It consists of six parameters including slope length factor (LS), erosion control practices factor (P), soil erodibility factor (K), cover management factor (C), and rainfall erosivity factor (R), which make it possible to calculate average yearly soil loss. The value of each factor is adapted for statistical and empirical data with the integration of GIS or is taken from the previous literature. The soil erosion can also be calculated by the integration of ArcGIS, MATLAB and SolidWorks software [19]. RUSLE in conjunction with GIS has become an effective tool for assessing soil erosion. Many researchers used the integration of RUSLE model with GIS in their study to estimate the rate of soil erosion [20,21,22]. Soil and water conservation practices (SWCPs) have a positive impact on soil properties. A relevant study has been conducted in Ethiopia [23].

A study was conducted by Sharma et al. [8] on the effect of LULC on erosion process in the reservoir from 1989–2004. This study illustrated that a slight increase in mean soil loss was observed. In 1989, it was 12.11 t/ha/y and in 2004 it became 13.2 t/ha/y. Results indicated that deforestation and increased wasteland in higher slopes increased the mean soil loss rate in 15 years. Another research in the Western Polish Carpathians shows that due to the increase in plantation and a decrease in the rate of cultivation, the rate of eroded soil decreased in the last 160 years (1846–2009). The rate of soil erosion in the year of 1846 was 18.13 tons/ha/y and it decreased to 4 tons/ha/y in the year 2009 [14]. Ouyang et al. [9] demonstrated the behavior of land-use changes, which shows that it has a more severe effect on erosion rate than changes in the properties of soil. For the estimation of soil loss in the Kelantan River basin, Abdulkareem et al. [24] used the USLE (universal soil loss equation) model based on GIS. Arable areas have a greater susceptibility to erosion process; research [25] illustrated that arable lands are approximately ten times more vulnerable to soil loss than orchards. In recent years, many studies have been conducted to analyze the impact of LULCC on soil erosion and sediment transport rate [14,15,20,26,27].

Soil loss is estimated by the nature of erosion process and data availability of the factors that contribute towards erosion. In the past, soil erosion estimation was done by using different methods ranging from factor-based approaches to process-based models. Now, RUSLE model is the most widely used considering its simplicity, and its input parameters are easily available. Investigation of impact of LULCC on soil erosion can be done by analyzing the effect of land use land cover changes (LULCC) on soil erosion potential and discharge of sediments through historic satellite images [8]. In this study, the main purpose is to figure out the potential of soil loss in the mountainous region, which changes due to an increase in agriculture production and decrease in the forest cover. For this purpose, a GIS-based RUSLE model is used to estimate the long-term effect of LULC on soil erosion. The results of the study are demonstrated in the form of temporal severity maps, highlighting the areas with high land use changes and their impact on soil erosion intensity. The quantitative estimation of soil loss due to LULCC is important for researchers and planners to take necessary actions to minimize the rate of soil erosion. Such an efficient planning will also be helpful to reduce the sedimentation in the reservoir of hydraulic dam(s) constructed on Chitral river, which drains through this watershed.

2. Methodology

2.1. Study Area

The northernmost part of Pakistan is the Chitral district which is situated in the world’s largest mountains (Figure 1). Its coordinates are 35°53′15′′ N and 71°48′01′′ E. Its neighbors include Afghanistan, China, Central Asian States, and Northern Areas of Gilgit. Chitral district is 322 km away from the provincial capital of Khyber Pakhtunkhwa. The DEM (digital elevation model) (Figure 2) clearly shows that the study area has rugged mountainous terrain with variation in altitude from 1053 m to 7695 m. In the north-east, it is bounded by the Karakoram, in the south by Hindu Raj Range and north-west by Hindu Kush mountains. Moreover, there are 40 peaks within an area of 14,850 km and altitude ranges from 1094 m at Arandu to 7726 m at Tirchmir. The Chitral district is composed of several valleys. The Chitral Mastuj valley is the largest and most important valley, expanding from the Afghan border to Arandu on the southern tip. Other important valleys include Terich, Shishi, Owir, Mulkhow, Lotkoh, Laspur, Torkhow, and Ashrat. The main Chitral valley is 354 km long and its width varies from 180 m to 4800 m. The whole study area is drained into Chitral river through several tributary channels. From the origin of the glacier Chianter, the river enters Afghanistan at Arandu. The study area is a mountainous tract that receives about 10–25 mm rain per month. The minimum and maximum temperature remain between 21 °C and 37 °C, respectively. Since high mountains prevent much of the monsoonal wind and moisture from reaching the study area, the summer season usually receives less amount of precipitation. Therefore, the region is not considered favorable for vegetation growth. On 16 July 1973, a maximum discharge of 1586 m³/s was recorded while March 10, 1964 witnessed a minimum value of 46 m³/s. The river flows undeviatingly the whole year because of the melting of glaciers and snow. During the season of monsoon, i.e., July–September, some extra runoff can be recorded. The river siphons a wide area of a steep slope and widespread area covered with snow, which contributes half of the discharge to the Kabul River. For irrigation and hydropower generation, a dam was constructed on the Chitral River having a maximum output of 250,000 kW. The location of the dam is about 30 km west of Peshawar. Due to the extensive silting, the reservoir is almost full. Throughout summers and winters, the power generation capacity of the dam has dropped to 64,000–20,000 kW, correspondingly.

2.2. LULC

For the assessment of LULC, Landsat satellite data were obtained. Landsat is the largest program run by NASA/USGS for earth’s imagery. Landsat 7 was used to develop the land cover maps of 2000 and 2010 while Landsat 8 was used for 2020. All three land covers were developed from satellite imageries having zero cloud cover, taken during the period of least snow cover in the area. Along with this data, practical samples were also collected from the field (study area). Based on Landsat satellite data and samples, supervised classification for the years 2000, 2010, and 2020 is made. Accuracy assessment is based on the following two parameters.

2.3. Precipitation

Precipitation data is required for the calculation of rain erosivity factor (R). This factor justifies the extent and intensity of every single rainfall throughout a year. This study utilized Satellite Rainfall Data of GPM (Global Precipitation Measurement). This satellite data is the product of NASA (National Aeronautics and Space Administration). It is widely used because of its enhanced accuracy. The annual rainfall data for the two meteorological stations situated in the study area (mentioned in Figure 1) was obtained from the Pakistan Meteorological Department. It was integrated with the satellite data to validate the results.

2.4. RUSLE Model

The soil loss due to water is most commonly estimated by using the RUSLE model because it is simple, with easy availability of input data, and wide-ranging of its applicability. In this study, RUSLE model is used in combination with GIS to calculate the LULCC. Modeling of soil erosion is an important part of current techniques used in the study of geomorphology along with other practices like photointerpretation, rainfall simulation, and GIS [28,29]. This is described by the following equation,

$$E=R.K.LS.C.P$$

(1)

where E is the amount of average annual soil loss (tons/ha/y), R is rainfall erosivity factor (MJ mm/ha/h/y), K is soil erodibility factor (tons/ha)/(MJ mm/ha/h), LS is a topographic factor (dimensionless), C is vegetation cover and management factor (dimensionless), and P is soil and water conservation factor (dimensionless).

Rainfall erosivity factor is the most important constraint of the erosion process [30]. Analytical calculation of the R factor is impossible because detailed long-term precipitation data is not available. So, R factor is estimated with the help of other methods by using rainfall data available for it [31]. Modified Fournier Index (MFI), the most commonly used empirical method, was used in this research

$$MFI={\displaystyle \sum }_{j=1}^{12}\frac{p_i{}^2}{P}$$

(2)

where, p_i is monthly i-month precipitation (mm) and P is yearly precipitation (mm).

Modified Fournier index and an R factor of RUSLE model are strongly linearly correlated [20]. Licznar [32] determined the correlation between the R factor of RUSLE and MFI for Poland and the power-law equation is defined as

$$R=\mathrm{0.2265.}MF{I}^{1.2876}$$

(3)

This equation has been verified, validated [8] and applied in many studies targeting different parts of Pakistan [4,33,34]. Therefore, it can be assumed that this equation will provide accurate results for the area selected in this study.

There is a minor change in monthly and annual rainfall of two different periods, so they have the same R factor. The R factors of both from equation 3 and erosivity index are compared. Latocha et al. [28] in the Polish Mountains and Sudety and Drzewiecki et al. [35] in Polish Carpathians calculated the quantity of soil loss by using RUSLE model. The factor K denotes the erodibility of soil associated with vulnerability to runoff rate and soil loss. The following equation is used to find K factor,

$$K={\mathrm{2.77.10}}^{-6}.M^{1.14}.\left(12-OM\right)+\mathrm{0.043.}\left(S-2\right)+0.033\left(P-3\right)$$

(4)

where M is particle size parameter (0.002–0.1 mm and 0.002–2.0 mm particle size), OM is organic matter percentage content, S is class of soil structure and P is class of soil permeability. Validation of Equation (4) was not possible due to the unavailability of data. Therefore, soil texture map of the area was used and the values were assigned according to [8,33,36].

LS is the long slope which shows the surface topography. In this study, the equation proposed by Mitasova et al. [37] is used to find out the LS factor. Mancino et al. [38,39] have already used this equation in their research. DEM is used for calculation of flow accumulation and steepness of the slope.

$$LS=\left(M+1\right){\left(As|22.13\right)}^m{\left(sin\beta |0.0896\right)}^n$$

(5)

where As is an area of un-slope drainage per width of contour, β is slope gradient (radians), m (0.4–0.6), and n (1.0–1.4) are exponential factors and they depend on the type of dominant erosion. ArcMap GIS software 10.2 is used to calculate the LS factor.

$$LS=POWER\left(\left[flow accumulation\right].cell size|22.13,0.6\right).POWER\left(sin\left(\left[slope\right]\mathrm{.0.01745}\right)|0.0896,1.3\right)$$

(6)

C represents the vegetation cover and management factor and it is used for the determination of the effect of management actions on soil erosion intensity. The NDVI indices have been used for the extraction of vegetative area, NDWI indices for water body extraction and NDBI for the built-up area. Accuracy of maps was validated from grouped samples. The erosion control practice factor (P) is generally calculated for areas with erosion control practices. Since no such practices were in place in the study area, a constant raster value is used for soil and water conservation factor P [31,40]. This value ranges from 0 to 1.

3. Results

3.1. Dynamics of Land Use Land Cover

General classification accuracy of 79%, 81% and 84% for the LULC maps of 2000, 2010 and 2020 respectively were obtained by using different image processing techniques in collaboration with hybrid classification techniques (

Table 1. Accuracy of land use classes.

Land Use Classes	2020		2010		2000
	Class Accuracy (%)	Overall Accuracy (%) and Kappa Coefficient	Class Accuracy (%)	Overall Accuracy (%) and Kappa Coefficient	Class Accuracy (%)	Overall Accuracy (%) and Kappa Coefficient
Agriculture in Valley	88	84% and 0.82	86	81% and 0.8	85	79% and 0.77
Agriculture in Sloping Valley	78		76		76
Bare Areas	89		87		84
Natural Herbaceous Shrubs	82		80		79
Natural High Shrubs	80		78		79
Snow and Ice	82		77		73
Natural Trees	79		78		77
Water Bodies	91		87		84

). The kappa coefficient of accuracy also improved throughout these years. The results of the year 2020 are more accurate as compared to the previous years. The classes of agriculture in the valley and water bodies have the highest accuracy.

The study area is classified into eight categories, that is, agriculture in sloping valley, agriculture in the valley, bare areas, natural herbaceous shrubs, natural high shrubs, natural trees, snow and ice, and water bodies. The study area is 1,449,120.7 ha and LULC was estimated for years of 2000, 2010, and 2020. Figure 3 clearly shows that the studied watershed is largely covered by natural herbaceous shrubs. The highly thick natural herbaceous shrubs and snow and ice are the two most plentiful land-use types, while a small percentage of the entire watershed is covered by bare areas, natural high shrubs, and natural trees. A very small area was occupied by agriculture in sloping valley, agriculture in valley and water bodies. Satellite examination helped in deriving the thematic maps of the last two decades (20 years). The composition or total extent of the individual LULC category/class is listed in

Table 2. Land use land cover changes from 2000 to 2020.

Land Use Classes	2020		2010		2000
	Area (Ha)	Area (%)	Area (Ha)	Area (%)	Area (Ha)	Area (%)
Agriculture in Valley	39,234.5	3%	28,982.4	2%	14,491.2	1%
Agriculture in Sloping Valley	51,107.9	4%	14,491.2	1%	14,491.2	1%
Bare Areas	161,053.9	11%	130,420.9	9%	115,929.7	8%
Natural Herbaceous Shrubs	515,637.5	36%	550,665.9	38%	565,157.2	39%
Natural High Shrubs	135,184.6	9%	86,947.2	6%	72,456.1	5%
Snow and Ice	465,076.1	32%	507,192.2	35%	521,683.6	36%
Natural Trees	80,471.7	5%	101,438.4	7%	115,929.7	8%
Water Bodies	1354.4	1%	28,982.4	2%	28,982.4	2%
Total	1,449,120.7	100%	1,449,120.7	100%	1,449,121	100%

.

It is obvious from Table 2 and Figure 4 that natural herbaceous shrubs and snow and ice remained the two most dominant land use classes throughout the study period. Together, these two classes occupy 75% of the area.

Table 3. Statistical analysis of land use classes.

Land Use Classes	Soil Erosion Year	Statistics (All Units in tons/ha/year)
		Mean	Min	Max	STD
Agriculture in Valley	2020	144.09	5.97	1578.242	203.87
	2010	122.66	4.21	1464.702	173.43
	2000	111.14	3.57	1370.052	143.54
Agriculture in Sloping Valley	2020	198.9	6.45	5681.93	562.02
	2010	155.25	4.91	5242.39	513.59
	2000	123.57	4.48	4865.9	468.21
Bare Areas	2020	177.04	5.03	5467.25	461.49
	2010	164.5	4.57	4963.82	414.54
	2000	152.15	4.12	4531.28	398.43
Natural Herbaceous Shrubs	2020	593.53	5.44	158,647.45	7067.94
	2010	605.96	5.87	159,302.21	6746.44
	2000	611.4	6.12	159,623.75	6783.98
Natural High Shrubs	2020	116.32	5.64	8421.65	403.38
	2010	104.89	4.32	7988.97	367.45
	2000	103.57	4.22	7701.54	313.54
Snow and Ice	2020	1276.08	5.4	335,719.74	15,947.62
	2010	1152.43	4.99	334,365.09	15,562.19
	2000	1033.03	4.53	333,126.55	14,907.43
Natural Trees	2020	124.84	6.173	4767.86	286.48
	2010	115.44	5.54	4511.21	234.65
	2000	112.01	5.12	4186.56	212.47
Water Bodies	2020	70.23	5.64	377.412	85.31
	2010	63.8	4.65	361.982	73.54
	2000	59.17	4.13	340.552	63.56

shows the changes in soil erosion of every individual class for all the studied years.

The area occupied by natural herbaceous shrubs is reduced to 515,637.5 ha (36%) in 2020 as compared to 550,665.9 ha (38%) in 2010 and 565,157.2 ha (39%) in 2000. This reduction is approximately 8.76% in comparison to the original area in 2000. The snow and ice is the second major LULC class during the period of study. It also shows decrement. The area covered by snow and ice is reduced to 465,076.1 ha (32%) in 2020 as compared to 507,192.2 ha (35%) in 2010 and 521,683.6 ha (36%) in 2000.

There is a gain in the bare area during the study period. Content of minerals and organic matter in soil is suitable for soil quality because it affects the soil functioning positively [8]. The bare areas are mostly rich in organic matter and mineral content as compared to the agricultural land. This gain in the area is 28% in 2020 as compared to the total bare area in 2000. Natural high shrubs increased to 135,184.5 ha (9%) as compared to 86,947.2 ha (6%) in 2010 and 72,456.1 ha (5%) in 2000. It shows an increment of 33.33% as compared to the original area in 2000. Water bodies of this region are reduced in these last 20 years. The water bodies decreased to 1354.4 ha (1%) in 2020 as compared to 28,982.4 ha (2%) in 2010 and 2000. Natural trees of the study area also showed a decrement of 30.5% as compared to the original area in 2000. These natural trees decreased to 80,471.7 ha (5%) in 2020 as compared to 101,438.4 ha (7%) in 2010 and 115,929.7 ha (8%) in 2000.

In contrast to all these land use classes, agriculture shows an increment in both valley and sloping valley. Agriculture in the valley increased to 39,234.5 ha (3%) in 2020 as compared to 29,892.4 ha (2%) in 2010 and 14,491.2 ha (1%) in 2000. The increase in agriculture in the valley is gradual. The increase in agriculture in the sloping valley is not gradual. In the first decade, there was no change in agriculture in the sloping valley, but it rapidly increased in the second decade. It increased to 51,107.9 ha (4%) in 2020 as compared to 1491.2 ha (1%) in 2010 and 2000.

Table 4. Land use class-wise soil erosion severity.

Land Use Classes	Soil Erosion Year	Soil Loss Class (%)
		Very Low	Low	Moderate	High	Very High
Agriculture in Valley	2020	41	16	17	19	7
	2010	38	23	29	6	4
	2000	32	26	31	7	4
Agriculture in Sloping Valley	2020	51	21	17	9	2
	2010	47	20	20	11	2
	2000	41	19	24	13	3
Bare Areas	2020	7	18	26	35	14
	2010	18	22	24	27	9
	2000	21	23	26	23	7
Natural Herbaceous Shrubs	2020	39	18	21	11	11
	2010	41	23	30	4	2
	2000	40	31	20	5	4
Natural High Shrubs	2020	50	22	17	9	2
	2010	53	20	15	10	2
	2000	51	19	17	10	3
Snow and Ice	2020	61	26	7	5	1
	2010	57	21	12	6	4
	2000	60	25	8	5	2
Natural Trees	2020	39	18	21	11	11
	2010	41	23	30	4	2
	2000	40	31	20	5	4
Water Bodies	2020	61	26	7	5	1
	2010	57	31	7	4	1
	2000	60	25	8	5	2

and Figure 5 show the land use class-wise soil erosion of every individual class in detail. It provides the ease to see and understand the changes throughout the study period because of its sub-divisions.

3.2. Estimation of Sediment Yield and Soil Erosion

The mean soil erosion of the entire watershed was 9.21 tons/ha/y in 2000, 12.43 tons/ha/y in 2010, and 15.63 tons/ha/y in 2020. Figure 6 and Figure 7 show the spatial variation and distribution of different RUSLE factors and also the resultant map of soil erosion for all three years. The Warsak Dam which is the reservoir part of the watershed does not contribute to the net soil erosion, but it behaves as a sink for eroded particles of soil which was excluded from further analysis of soil erosion. The range of R factor was found to be 349.769 MJ mm ha⁻¹ h⁻¹ year⁻¹ as highest and 197.347 MJ mm ha⁻¹ h⁻¹ year⁻¹ as the lowest value, and variation in its spatial view is shown in the Figure 6a. The range of K factor value lies between 0.1 and 0.4. The spatial distribution of the K factor is shown in Figure 6b. The calculated values of the LS factor using SRTM ranges from a minimum value of 0.2 for flat terrain to a maximum value of 6.3 for high-elevation areas, predominantly for hill slope region. The spatial distribution of the LS factor is shown in Figure 6c.

The spatial distribution of p factor is not generated because there are no precautionary or conservation measures being taken against the soil erosion in this region. The values of four parameters of RUSLE out of five remain constant but the value of management (C factor) shows variation both spatially and temporally. The spatial and temporal distributions of the C factor are presented in Figure 7a–c. The values of the C factor range from 0 to 1 for all three years of study, but changes in LULC make their spatial distributions different. Sediment yield of the watershed was obtained by putting all the five parameters into the RUSLE, and the obtained results are illustrated in Figure 8a–c. The figures assist us to extract the information that the central part of the study area is predominantly affected by the sediment yield. The percentage of the affected area kept changing throughout 20 years. In the year 2000, it was 20%, 50% in 2010, and then changed to 80% in 2020 (

Table 5. Soil erosion severity classes showing the area affected from 2000 to 2020.

Soil Erosion Severity Class	Soil Loss (tons/ha/year)	2020		2010		2000
		Total Area (ha)	Area (%)	Total Area (ha)	Area (%)	Total Area (ha)	Area (%)
Very Low	<5	623,122	43%	536,175	37%	521,684	36%
Low	5–10	173,895	12%	275,333	19%	333,298	23%
Moderate	10–20	304,315	21%	362,280	25%	362,280	25%
High	20–50	231,859	16%	202,877	14%	173,895	12%
Very High	>50	115,930	8%	724,56	5%	579,65	4%
Total		1,449,121	100%	1,449,121	100%	1,449,121	100%

). Due to anthropogenic activities, all those areas which fall in the high class in 2000 and 2010 have been transformed into a very high soil erosion hazard zone. Hence, it can be said that man-made activities further increase the potential of soil erosion in areas which are already vulnerable to soil erosion. This increased percentage of the affected area is alarming and needs to be addressed properly.

The soil erosion of the studied watershed showed a variation on the spatial scale to some extent in all three years of this study. Singh et al. [41] suggested that the soil erosion maps of the studied years are to be classified into four soil erosion classes from a management perception such as low, medium, high, and very high. The area affected by the soil erosion was divided into five categories (Figure 9) to have a better and clear understanding. Quantification of the soil hazard is very important to identify the hazard level and hotspot regions. Usually, hazard is defined in terms of very low, low, moderate, high and very high [42]. Ranges for the soil loss tolerance level for the study area are defined by OECD [43]. Multiple researches, carried out in several parts of the region, same as study area, followed the similar soil erosion classification method [44,45,46]. Therefore, similar quantification technique is followed in this study while keeping ranges of this region in mind. It was categorized as very low (less than 5 tons/ha/y), low (5–10 tons/ha/y), moderate (10–20 tons/ha/y), high (20–50 tons/ha/y), and very high (more than 50 tons/ha/y). The results of the difference between the three maps are presented in Table 5. The area affected by soil erosion showed increment. It is obvious from the increment that the land cover use kept changing throughout these years. There is a net change in the total area under the very low, low, moderate, and high categories by 7%, −11%, −4%, and 4%, respectively. The percentage of affected areas under a very high category was 4% in the year 2000, and it changed to 5% in 2010 and 8% in 2020. The area under very high category is increased by 4% in the comparison between 2000 and 2020, which is a major issue of concern. The change in soil erosion and its spatial distribution is depicted in Figure 10.

The study area is a mountainous region. The figure shows that changes were severe in areas with high topographic potential for soil erosion during the study period. Degraded forestland, predominantly in foothill regions due to increasing human activities, is one of the main reasons behind this severe soil erosion. Invulnerable areas, land use changes in the intensification of agricultural activities and deforestation make the land more prone to soil erosion. From these results, it can be inferred that any land transition to cropland would be harmful, as it was the major source of sedimentation. On the other hand, the forest was the most effective barrier to soil loss.

With the passage of time and advancement in the knowledge of land, the cover has become essential to overcome the issues related to biogeochemical cycles, biodiversity, deterioration of environmental quality, loss of productive ecosystems, and loss of agricultural lands. The main reason behind the LULC changes includes rapid growth in population, rural-to-urban migration, reclassification of rural areas as urban areas, lack of valuation of ecological services, poverty, ignorance of biophysical limitations, and use of ecologically incompatible technologies. The present study area of Chitral is a mountainous and developing town. During the past few decades, the study area had witnessed substantial increase in population, economic growth, industrialization, and transportation activities that harmed environmental health of the region.

4. Discussion

The accuracy of the generated maps has improved throughout the study period. One of the main reasons is that the advancement in technology has enabled us to produce improved GIS, which can encounter more complex problems and issues. The results of the year 2020 are more accurate as compared to the previous years. The study conducted by [8] shows the same trend of improvement in accuracy but another study [47] shows that the accuracy for the maps in 2002 was 81.6% and it decreased to 80.5 percent for the year 2009. Thus, accuracy may depend upon one’s personal skills and availability of data.

Due to the involvement of multiple data sets, we used the latest technologies like remote sensing and GIS to quantify LULC. Interpretation of remote sensing imagery, GIS and existing study area conditions enabled us to classify the study area into eight categories, that is, agriculture in sloping valley, agriculture in the valley, bare areas, natural herbaceous shrubs, natural high shrubs, natural trees, snow and ice, and water bodies. To determine the hazard potential and vulnerable hotspots, such zonation has been made in several studies. Those studies are taken as base line. Four concepts (regenerative economics, nature-based solutions, connectivity, and systems thinking) were introduced by Keesstra et al. [48] to explain the neutrality of land degradation. A study in the past includes strategies to reach UN Sustainable Development Goals more effectively [49]. Visser et al. [50] discussed the process to achieve Sustainable Development Goals (SDGs) by using transitions for soil-water system. Another study was conducted in the past to examine the geography of soil [51]. The classification of same nature has already been made in previous researches [8,47]; some other classes can also be found in the literature [52].

Table 2 shows the area covered by every land use class for the studied duration. It shows reduction in area covered by agriculture in valley, agriculture in sloping valley, bare areas, and natural high shrubs. The area occupied by natural herbaceous shrubs was reduced. The possible reason for the reduction in natural herbaceous shrubs might be the degradation of land because of nutrient depletion by intense soil erosion due to the lack of required soil conservation measures. The major class of snow and ice also shows decrement. This reduction is the consequence of global warming, the greenhouse effect, and the changing climates due to the mentioned reasons [53]. The class of bare area shows an increased percentage. It might be due to the settlement of more people in urban areas which leaves the rural areas bare [54]. Degraded forest, river sand, bare soil, heavily eroded land, and rock outcrop are included in the bare areas described in this study [8]. Natural high shrubs increased because the study area is governed by high altitude, receiving comparatively more rainfall, which excels the growth of natural high shrubs. The area covered by the water bodies is reduced due to lack of conservation practices. The agriculture in the area is increased, which is not very unusual considering the increasing demand of population and food. The highest soil erosion has a strong connection with cultivated land [8]. The strata of the cultivated soil are disturbed when it is ploughed. This disturbance of soil strata makes it more vulnerable to soil erosion. The increase in agricultural production is possible by the use of fertilizers turning uncultivable land to cultivable lands. Our study shows contrast to many other studies in which agriculture area or cultivated land is decreased due to many reasons such as conversion of agricultural land to settlements and to other land use classes such as forest [14].

Different land use types have great impact on soil erosion. The rate of soil erosion may change from area to area. It is not fixed that it will increase every time. It depends on conservation practices and management against soil erosion. Calculation of soil erosion can be divided into two categories, i.e., one for short period and the other for long period. Studies for short-term soil erosion rates are conducted in the Polish Carpathians, [35,55,56,57,58,59,60]. The Szymbark IG&SO PAS Research Station (49°38′04″ N, 21°07′08″ E), located in the Polish Carpathians, where soil erosion measurements covered the last 30 years is an exception [61,62,63,64]. A similar analysis was conducted by Latocha et al. [30] for soil erosion for a long-term period (in the last 150 years). However, their study was in the Sudetes Mountains, where after World War II, immigration actions led to depopulation and then forest expansion. This consequently decreased the soil erosion rates [14]. Our study is for the short period of 20 years and it shows that the rate of soil erosion increased from 2000 to 2020. The increment in the mean soil erosion may be due to slight transitions in land and land use. The main reason for increase in the mean soil erosion of the watershed over the studied period could be allocated to the drop in percentage of forest area, the increment of bare areas, and increased cultivation practice in areas that are more prone and vulnerable to soil erosion. It is important to remember that the change in the mean soil erosion may significantly affect the sedimentation process of the reservoir. LULC dynamics in the watershed result in more sedimentation in the reservoir, which would affect its life negatively. LULC dynamics obtained in this study indicate that the studied watershed is affected by moderate management and some other local human activities. RUSLE model was used to evaluate the subsequent effects of LULC dynamics on soil erosion potential.

This study shows increment in soil erosion rates. Changes in land use and land cover contributed to this increment. The LULC change was rapid in the second decade of the 21st century in the study area. It shows that the area is under many human activities such as a reduction in natural trees, snow, and ice, and increment in agriculture in plain and sloping areas as well. Increase in soil erosion can also be attributed to changes in the land cover, which are reflected by NDVI, and in turn changing the C factor values which strongly affect the soil erosion. Global climate system may be altered due to LULC changes. Nearly two-thirds of the precipitation that falls on the earth’s surface is returned to the atmosphere via evaporation. Precipitation amount and intensity and potential evapo-transpiration is directly affected by climate change, and it indirectly affects plant water-use efficiency through altering plant growth rate and species composition, which result in change of land cover. The figures show that the areas with high soil erosion potential were mostly located along the foothill regions in the southern part and some portions of the central watershed during studied years. This study shows contrast to many previously conducted studies such as [14,47]. In these studies, the rates of soil erosion decreased with changes in land use. It is obvious from all these references that soil erosion rates may change and differ for different land, land uses and land covers depending upon the topography, settlement trends and meteorology of the area. The outcomes of this study will help the stakeholders and regulatory decision makers to control deforesting, controlled agricultural farming and to take other necessary actions to minimize the rate of soil erosion. Such an efficient planning will also be helpful to reduce the sedimentation in the reservoir of hydraulic dam(s) (Warsak Dam) constructed on Chitral river, which drains through this watershed. Apart from its local implications, this study may also provide useful guidelines for the other parts of the globe having similar climate and geographical settings.

5. Conclusions

The study presented in this paper emphasizes on rapid assessment of LULC changes and their consequences on regional soil erosion potential in a mountainous watershed of Chitral, Pakistan. For this purpose, RUSLE model is used in combination with GIS, which is proved to be an effective tool to extract land use land cover changes. The results of the study demonstrate that LULC changes were significant during the period from 2000 to 2020. The risks of soil erosion have been increased by 4%. The mean erosion potential for the whole watershed was increased from 9.2 tons/ha/y in 2000 to 15.63 tons/ha/y in 2020. A significant expansion in bare areas, agriculture in the valley, and agriculture in sloping valley has been noticed. On the other hand, there is a decrease in areas covered by water bodies and natural trees. These results indicate the significant impact of different land cover factors and human activities on LULC change, and consequently, on soil erosion.

Massive deforestation and increased agricultural land could be treated as most important reasons for increase in soil loss during the study period. This demands an immediate attention of the stakeholders, regulatory decision makers, and environmental management groups to propose comprehensive guidelines for possible LULCC. The farmers need to be provided with the capacity development programs to raise their level of awareness. An efficient system to control soil erosion will be equally helpful to reduce the sedimentation and increase the life and efficiency of hydraulic dam(s) constructed on Chitral river, which drains through this watershed. The soil loss hazard maps generated in this study could also be helpful to develop strategies to overcome food and water problems in the coming decades.

6. Limitations and Future Studies

In the future, the effect of LULC changes can be estimated in other hilly and plain areas. Effects of all other influencing factors like slope, elevation, and rainfall should also be estimated by integrating RUSLE and GIS. Global climate system may be altered due to LULC changes. Nearly two-thirds of the precipitation that falls on the surface of earth is returned to the atmosphere via evaporation. The prediction about the changes in the precipitations can be made depending upon LULC changes in the study area. Additionally, LULC alters the amount of CO₂ in the atmosphere. A study can be conducted to estimate the change of CO₂ in the atmosphere due to LULCC. The storage capacity of the Warsak Dam is decreasing day by day due to changes in LULC which directly contributes to soil erosion. Thus, there is a need to propose some control structures in Chitral river by using BIM (building information modeling) to control soil erosion.