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

: The world’s ecosystem is severely a ﬀ ected 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 signiﬁcant cause of sediment transport. This study aims to estimate the e ﬀ ect 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). Di ﬀ erent factors were analyzed to study the e ﬀ ect 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 a ﬀ ected 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 e ﬃ cient 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.


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

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 3 /s was recorded while March 10, 1964 witnessed a minimum value of 46 m 3 /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 Sustainability 2020, 12, 5898 4 of 20 full. Throughout summers and winters, the power generation capacity of the dam has dropped to 64,000-20,000 kW, correspondingly.

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.

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.

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.

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.

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, = . . . .

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, 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 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 = 0.2265.MFI 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 = 2.77.10 −6 .M 1.14 .
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.
where As is an area of un-slope drainage per width of contour, β is slope gradient (radians), m 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.

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). 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 Sustainability 2020, 12, 5898 7 of 20 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.     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 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. 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.  Table 4 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.

Estimation of Sediment Yield and Soil Erosion
The mean soil erosion of the entire watershed was 9. year −1 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.

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. Figures 6 and 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 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). 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 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). 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  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 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.

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

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.

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.

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 2 in the atmosphere. A study can be conducted to estimate the change of CO 2 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.