The impact of land use and land cover changes on soil erosion in western Iran

Analysis of long-term land use and land cover (LULC) changes requires up-to-date remotely sensed data to assess their effects on erosion. This is a particularly important assessment for regions with semi-arid landscapes where soils tend to be scarce and proper management necessitates matching LULC to local conditions to achieve sustainable land use. This study evaluates the impact of LULC changes on erosion using Landsat satellite imagery and the Revised Universal Soil Loss Equation model on the plains around the Jarahi River and the Shadegan International Wetlands. Supervised-classification and maximum-likelihood methods were applied to pre-processed TM, ETM, and OLI images for 1989, 2003, and 2017 to prepare LULC maps. The areal extents of agricultural land, wetlands, waterbodies, and built-up regions increased by 9.48%, 2.52%, 3.44%, and 0.13% from 1989 to 2003, respectively. During this timeframe bare (or barren) lands and wetland vegetation decreased by 11.44% and 4.13%. Between 2003 and 2017, however, the areas of bare lands, waterbodies, and built-up areas increased 12.77%, 1.52%, and 0.30%, while agricultural lands, wetlands, and wetland vegetation decreased by 9.99%, 1.32%, and 3.27%. According to the results, the areal extent of erosion at a rate > 1.1–5 Mg ha−1 year−1 has been increased by about 45.56%, 50.06% and 52.24 between 1989, 2003 and 2017. LULC changes led to increased soil erosion on agricultural and bare lands. This highlights the need to plan and manage changes to LULC to reduce erosion to and below sustainable levels. Nature-based solutions can be effectively used to reduce erosion.


Introduction
Changing land uses and land covers (LULC) are the primary environmental changes responsible for global change (Guan et al. 2011;Savari et al. 2020). These changes are principally due to deforestation, urbanization, intensification of agriculture, and overgrazing and they lead to land degradation. Natural fluxes like long-term temperature and rainfall changes can also drive LULC changes (Lambin 1997). Intensification of agriculture and excessive livestock grazing are significant triggers for desertification and land degradation in arid regions (Daliakopoulos et al. 2016). In sensitive and fragile areas, degradation and desertification reduce productive capacities of the land (Eskandari et al. 2016); there is, therefore, a need for better approaches to management. Controlling erosion is crucial to achieving sustainable development. Clark and Dickson (2003) argued for the need to develop sustainability science to facilitate sustainability of societies. Sachs and McArthur (2005) advanced the millennium project to build millennium development goals. Anthropogenic LULC changes can destroy natural resources and affect food supplies to the extent that they may cause severe social, economic, and political consequences (Khosravi et al. 2017;Eskandari-Dameneh et al. 2020;Turner et al. 2007). Studies have focused primarily on LULC changes. Fewer have examined the relationships between land-use change and environmental impacts like soil erosion (Xiao et al. 2006;Rawat and Kumar 2015;Hegazy et al. 2015). In general, inappropriate land use changes can accelerate erosion (Chen 2008). Recently, several studies have analyzed the environmental effects of erosion, reviewing its long-term impacts on soil fertility, crop yields, and soil quality (De Wit and Behrendt 1999;Verstraeten et al. 2002). Land-use change due to anthropogenic influences is one of the main factors driving soil erosion (Houben et al. 2006). It is crucial to examine the potential impacts of LULC-change induced erosion locally and regionally. LULC-change is one of the most important processes diminishing soil fertility, but it is also tied to other problems like flooding, salinization, and water contamination (Rickson 2014;Xiubin and Juren 2000).
Erosion by water is a process that separates, moves, and damages a soil's content and structure. This process can occur naturally or can be initiated and intensified by human actions. Erosion rates can change dramatically depending on soil properties and environmental and climate conditions. Controlling soil erosion and its consequences (soil loss, slope instability, and reduced fertility) depends on good land use management (Bini et al. 2006). Studies have shown a strong correlation between land-use change and soil erosion (Mutua et al. 2006;Sharma et al. 2011). Models have been used to calculate soil erosion (average long-term soil losses). One important model is the Revised Universal Soil Loss Equation (RUSLE) (Renard et al. 1997). The RUSLE model calculates the maximum amount of allowable soil loss from a specific soil type before the maximum sustainable yield of an area is diminished (Ranzi et al. 2012). It can be used to determine appropriate land-use systems and can guide erosion prevention and mitigation efforts (Ranzi et al. 2012;Zare et al. 2017a, b).
Studies have examined the impacts of land-use changes on annual erosion rates using the USLE model-in the Emilia-Romagna highlands of north-central Italy (Brath et al. 2002)-and the RUSLE with GIS (Geographic Information System) in the Wadi Karaka watershed of Jordan (Farhan and Nawaiseh 2015). Isaaca and Aqelel Ashraf (2017) reviewed the relationships of erosion and land degradation to water quality and concluded that erosion adversely impacts water quality. Sharma et al. (2011) examined the impact of LULC changes on erosion potential in agricultural lands (Sharma et al. 2011), andTadesse et al. (2017) assessed LULC-change induced erosion impacts in northeastern Ethiopia (Tadesse et al. 2017). Zare et al. (2017a; used RUSLE and the Conversion of Land Use and Its Effects at Small Regional Extent (CLUE-s) model to examine scenarios of land-use change and their impacts on soil erosion in the Cazilan watershed, Iran.
Considering the importance of agriculture in the Shadegan wetland of Iran, the susceptibility of the ecosystem to environmental change, and the emergence of airborne dust in the region, this study examines the relationships of LULC change to soil erosion in this region using GIS techniques and remote sensing. A motivating goal of this work is to produce useful information for managers and decision-makers who are seeking management solutions and conservation practices to combat erosion in the region.

The study area
The study area is in Khuzestan Province, Iran. It encompasses the plains around the Jarahi River in Ramshir and Shadegan counties and the Shadegan International Wetland. The region with an area of 299,000 hectare is located between 30° 19′ and 30° 51′ N and 48° 41′ and 49° 43′ E, (Fig. 1). As the 11th longest river in Iran (438 km long), the Jarahi River flows through Kohgiluyeh and Boyer-Ahmad Province and Khuzestan Province. The main area under cultivation of agricultural products in winter is wheat and rapeseed, and beans, vegetables and alfalfa are also cultivated very little. The climate of this region is dry and the average rainfall is 218 mm per year and the temperature is 26 °C. The Shadegan Fig. 1 Location of the study area in Iran International Wetland is located at the delta of the Jarahi River on the Khuzestan Plain. Shadegan is the home of permanent and seasonal hydrophilic and mesophilic vegetation. Halophilic plants dominate in the swamp area, but there are also palm orchards and agricultural lands. The four main species of these plants are Cyanus longus (covering about 70% of the wetland area), Typha minima (about 15-20%), Salsola sp (about 10%), and Phragmites sp (about 5%) (Rahimi Blouchi et al. 2013). Elevations range from 3 to 59 m and slope. angles range from 0 to 60%.

Data
Satellite data (path and row 165 and 39) were compiled from Landsat 5, the TM(Thematic Mapper) sensor (on June 3, 1989), Landsat 7, the ETM(Enhanced Thematic Mapper) sensor (on May 29, 2003), Landsat 8, and the Operational Land Imager (OLI) sensor (May 23, 2017). A digital elevation model (DEM) was also acquired from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) data at a 30 × 30 m resolution. These data were accessed from the United States Geological Survey (USGS) website (www. earth explo rer. usgs. gov). The digital elevation model (DEM) is necessary for providing the digital Earth. Nowadays, there are two sets of near-global DEMs produced using remotely sensed data. One is the elevation dataset produced using the C-band singlepass Interferometry Synthetic Aperture Radar (InSAR) data obtained by the Shuttle Radar Topography Mission (SRTM) covering 56° S to 60° N. The other is the Global Digital Elevation Model produced by the stereo processing of the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) covering Earth's land surface between 83° N and 83° S (Ni et al. 2015). A 1:250,000scale topographic map, soil texture map, and climate data  were acquired from Iran's National Cartographic Center, the Agricultural Research Center of Khuzestan Province, and the Khuzestan Meteorological Organization, respectively. The layers were created and merged using ERDAS IMAGINE and ArcGIS10.2 software.

Pre-processing of satellite imagery
A geometric correction was applied to the satellite data to facilitate ground accuracy. The TM and ETM sensor data were georeferenced to the OLI in plural using the imageto-image technique with an RMSE of less than 0.5 pixels. Since changes in illuminance affect the radiation intercepted by a pixel, atmospheric corrections were also applied. The ATCORE extension was used in ERDAS IMAGINE 2014 software and the metadata file associated with the satellite imagery for the atmospheric corrections.

Preparation of land-cover maps
Google Earth images, aerial photographs for periods 1989 and 2003, and Global Positioning System (GPS) point locations captured in the field for period 2017 were used to select training samples to carry out a supervised classification. The samples were randomly chosen from the specified area to recognize the type of land cover, using the Region of Interest (ROI) tool provided by ERDAS IMAGINE 2014 software, Google Earth, and field data. Half of the sample pixels were randomly selected as training samples, and the remaining half was used for classification accuracy assessment. The number of the sample pixels used for the classification accuracy estimation was 80 pixels for agricultural lands, 200 pixels for bare lands, 50 pixels for wetlands, 40 pixels for water bodies, 236 pixels for the wetland vegetation cover, 40 pixels for built-up areas.
The maximum likelihood algorithm (Ozesmi and Bauer 2002;Gumel et al. 2020) was used for classification in ERDAS IMAGINE 2014. This method is based on the probability that a pixel belongs to a particular class. The theory assumes that these probabilities are equal for all categories and that the input bands have normal distributions. However, this method requires a long period of time to complete computations, it relies heavily on a standard distribution of the data in each input band, and it tends to over-classify signatures with relatively large values in the covariance matrix. An appropriate band combination was identified for classification by selecting Evaluate in the signature editor menu based on the best average. This classification method is known as a powerful yet simple method due to the fact that training points are carefully selected and introduced to the software (Yousefi et al. 2011;Biro et al. 2013; Pardahan and Soleimian 2009). Band combinations were used to classify data obtained from the TM, ETM, and OLI sensors, respectively )Eskandari Damaneh et al. 2020a, b). Six classes of LULC were identified: agricultural lands, bare lands, wetlands, water bodies, wetland vegetation cover, and built-up areas (Table 1).

Estimation of the accuracy of the classified land-cover maps
The accuracy of the classified image was evaluated using the data not used for classification. The accuracy was assessed with an error matrix by calculating the user's accuracy, producer's precision, overall accuracy, and Kappa coefficients. The user's accuracy (U a ), producer's accuracy (P a ) overall accuracy(A o ) were calculated as the sum of elements of the main diagonal line of the error matrix (Eqs. 1, 2, and 3) (Jupp 1989;da Cunha et al. 2021;Das et al. 2021): where OA is overall accuracy, n is the number of experimental pixels, and ΣPij is the sum of elements of the main diagonal of the error matrix. Pij = number of correctly classified pixels in row i and j. Ai and Cj = total number of pixels in row i and j. The Kappa coefficient regards incorrectly classified pixels and calculates the classification accuracy compared to a wholly random classification (Pal et al. 2017;Mitsova et al. 2011). The Kappa coefficient was calculated according to Eq. 4: K Kappa Coefficient, N Total number of pixels, n Number of classes.

Evaluation of soil erosion
The RUSLE model was used to estimate the average annual soil erosion. This model contains six parameters: soil erodibility (K), rainfall erosivity (R), vegetation cover (C), the length (L) and steepness (S) of slopes, and conservation practices (P) in place. Sensitivity to erosion depends on soil characteristics. Changes in soil characteristics are related to LULC and topography (Pradhan et al. 2012). Soil erosion is calculated using Eq. 5 (Wischmeier and Smith 1978): where A is average soil erosion (Mg ha −1 year −1 ), R is rainfall erosivity (MJ mm ha −1 year −1 h −1 ), K is soil erodibility (Mg h MJ −1 mm −1 ), LS is the topographic parameter, C is vegetation cover, P is the set of conservation practices in use. Each of these measures is described below.

Rainfall erosivity (R)
The rainfall erosivity factor was proposed by Wischmier and Smith to include the effects of weather on soil erosion and is defined as the potential for rainfall to cause erosion. It depends upon the physical properties of raindrops and is associated with direct energy in impact, kinetic energy of raindrops, and maximum 30-min rainfall intensity (Wischmeier and Smith 1978). There are few meteorological stations equipped with rain gages in the study area. Therefore, annual and monthly rainfall-based indices like the Fornier Index (an indicator of rainfall "aggressiveness") were used in the USLE and RUSLE models (Ferro et al. 1991;Renard and Freimund 1994). A modified Fornier Index was calculated for this study using Eq. 6. Inserting this index into Eqs. 7 and 8 is a way to calculate R for areas without detailed rainfall data (Renard and Freimund 1994). The CRU database rainfall networking data were used to evaluate the R value. The CUR model is a dynamic-simulation model that uses ground-based meteorological data recorded worldwide (Harris et al. 2014). In addition to high spatial power, this database covers a longer period of time than other global data (Harris et al. 2020). In this study, monthly rainfall data between 1991 and 2017 were used to determine the value of R.
The Inverse Distance Weighting (IDW) method was used to generalize the entire study area (Renard and Freimund 1994): where Pi is the average monthly precipitation (mm) in a month i, P is the average annual precipitation (mm), R factor is in MJ mm ha −1 h −1 year −1 .

Soil erodibility
Soil erodibility indicates the inherent sensitivity of soil to erosion and the ease of dispersion of soil particles due to raindrops' kinetic energies and transportation of particles by runoff force (Khademalrasoul and Amerikhah., (2021). Characteristics such as soil texture, amount of topsoil, soil structure, and soil permeability directly impact soil erodibility (Gelagay et al. 2016). The more resistant the soil is to erosion, K approaches 0. For soil that is less resistant to erosion, K approaches 1 (Marondedze et al. 2020). Given that in this study soil descriptors like texture and soil organic matter were not available to calculate the K erodibility coefficient, therefore, global soil data from (www. soilg rids. org) were used (Hengl et al. 2017). The following Eq. (9) was used to calculate the soil index.

Topographic factor (LS)
The slope (S) (measured as %) and its length (L) (measured in meters) influence erosion. Flow accumulation is the accumulation of upslope flows for each cell, cell size is the network cell size (30 m), and slope was derived from the DEM. The constant 0.01745 was used to convert slope measures from degrees to radians in GIS.

Vegetation-cover factor (C)
Vegetation-cover (C) is the loss ratio of soil from a region with a specific vegetation cover to a plot in tilled farmland without plant residue (Wischmeier and Smith 1978). C was calculated with Eq. 11 (Lin et al. 2002): where NDVI is the normalized difference vegetation index, calculated by Eq. 12: NIR and RED are the reflectance values of a location in the near-infrared and red bands. This index ranges from − 1 to + 1. Negative and zero NDVI values refer to no vegetation, and higher values indicate dense vegetation (Drisya and Roshni 2018; Eskandari Damaneh 2021).

Protection support practice (P)
Protection support practice (P) is the soil loss ratio produced using a specific tillage practice that may promote or battle erosion. Straight-row cultivation is the worst approach in terms of erosion, particularly if rows are oriented downslope. Conservation practices include contouring, strip farming, or terracing. P is lower as conservation practices are increasingly effective at preventing soil erosion, and less soil loss occurs (Wischmeier and Smith 1978). The P index is dimensionless, which shows the effects of different soil protection methods against other erosion. The value of this index is between 0 and 1. Values close to 0 indicate complete soil protection against erosion, and values up to 1 show an increasing lack of soil protection against erosion and increased erosion (Yang et al. 2003). Also, the possibility of growing soil erosion increases due to increased tillage (Obiahu and Elias 2020). The values of the P factor were determined by reclassifying land cover ( Table 2). As no conservation practices were used in the study area, P was determined by land cover classes.

The impact of LULC changes on soil erosion
To determine the effect of LULC change on soil erosion, a map of land cover for each year compared to the map of soil erosion during the same year. Using the natural breaks (Jenks) classification method (Jenks et al. 1971;do Carvalhal Monteiro et al. 2018), the soil erosion intensity map was categorized into five soil erosion classes, i.e., very low (< 0.5), low (0.51-1), medium (5-15), high (2.1-5) and very high (> 5).

3
The erosion rates and classes of erosion (Table 3) were determined for each land cover class.

Evaluating classification accuracy
The assessment of LULC class efficiencies (Table 4) reveals that the highest coefficients of overall accuracy and kappa were 93% and 0.91 for the classification for 2017, and the lowest were 88% and 0.86 for 1989.

Rainfall erosivity (R factor)
R was calculated for the CRU precipitation data period 1991 to 2017 (Fig. 4). It ranged from 40 to 89 MJmm ha −1 year −1 h −1 across the study area.

Topography (LS factor)
LS was mapped based on the DEM and Eq. 7 (Fig. 6). LS ranges from 0 to 41 in the study region.

Vegetation cover (C factor)
C was determined by combining NDVI with Eq. 10 and it was mapped for 1989, 2003, and 2017 (Fig. 7). C ranged from 0.09 to 0.85 for 1989, 0.13 to 0.76 for 2003, and 0.15 to 0.82 for 2017. The highest and lowest values correlated to areas lacking vegetation (i.e., bare land) and densely vegetated areas, respectively.

Protection support practice (P factor)
P was mapped for each year by reclassifying land cover classes and assigning the respective numbers from Table 3 (Fig. 8). The values of P ranged from 0.12 to 1 in all three years.

Annual soil erosion
The five factors were layered in ArcGIS and were combined to determine annual soil erosion rates and annual soil loss. The values were grouped into five "risk" classes and these were mapped for 1989, 2003, and 2017 (Fig. 9). The total area of each soil erosion class was calculated (Fig. 10). The trends among the soil erosion classes varied and changed direction from 1989 to 2017. During the first period (1989 to 2003), areas classified as very low, low, moderate, and high increased by 1.38%, 1.42%, 3.27%, and 1.2%, respectively, while the very high soil erosion regions decreased in size by 7.29%. From 2003 to 2017, the areas of very low, low, and high erosion increased by 1.92%, 5.53%, and 4.49%, respectively, areas of moderate and very high erosion decreased by 2.31% and 9.63%, respectively. During the third period, which encompassed the entire study period 1989 to 2017, areas classified as very low, low, moderate, and high increased by 3.29%, 6.94%, 0.96, and 5.71%, respectively, but areas of very high erosion decreased by 16.92%.

Effects of LULC changes on the trends of soil erosion rates
Among the land uses examined, soil erosion has increased on agricultural and bare lands. It did not increase in wetlands, built-up areas, waterbodies, or areas covered by wetland vegetation. Classification of bare-land erosion rates during the first period (1989 to 2003) showed that there was an 11.36% increase of lands classified as very low, 18.66% increase of land classified as low, 12.10% increase of lands classified as moderate, and 0.34% increase of lands classified as high (Fig. 10). The area of lands classified as very high decreased by 42.47%. During the second period (2003 to 2017), areas classified as very low diminished by 11.92%, low diminished by 8.79%, and moderate diminished by 9.38%. The high and very high classes increased in size by 9.75% and 19.44%, respectively. During the entire period (1989 to 2017), the areas classified as very low increased by 0.35%, low increased by 9.87%, moderate increased by 2.72%, and high increased by 10.9%. The area of the very high class decreased by 23.02%. Fig. 10 The percentages of the five erosion classes across all LULC types in each year (brown, yellow, and green) and the percentages in each erosion class for bare lands (a) and agricultural lands (b) Erosion patterns associated with agricultural lands changed over the course of the period from 1989 to 2017. During the period from 1989 to 2003, LULCs classified as very low and low erosion increased in area by 12.38% and 0.44%, respectively, while moderate, high, and very high erosion areas decreased by 2.99%, 6.66%, and 3.17%, respectively. During the second period from 2003 to 2017, very low and low areas decreased by 13.56% and 4.98%, while moderate, high, and very high areas increased by 2.03%, 9.98%, and 6.52%. During the entire period, very low, low, and moderate zones decreased by 1.18%, 4.54%, 0.96%, and 10.9%, respectively. The areas classified are high and very high increased by 3.32% and 3.36%. Areas of very low, low, and moderate erosion decreased by 3.25%, 8.38%, and 5.49%, respectively. High and very high zones increased by 0.88% and 16.19% over the entire period.

Discussion
Soil erosion and sediment transportation are natural processes that have been accelerated by forest fires (Di Prima et al. 2018), grazing (Antoneli et al. 2018, agriculture , and road and railway construction (Hazbavi et al. 2018). Models are vital tools to discern soil erosion processes, and they inform us about the temporal and spatial changes taking place within regions. Soil erosion is also a significant driver of land degradation globally and throughout Iran (Ahmadi 2006;Rahman et al. 2009). This study investigated the impacts of LULC changes on soil erosion. The overall classification accuracy and kappa coefficient for these 28 years are 88% and 0.86, respectively, indicating that the classification by the maximum-likelihood algorithm of land-use change determination was highly accurate (Lu et al. 2019;Eskandari Damaneh et al. 2020a, b). The assessment of LULC changes from 1989 to 2003 revealed that bare land and wetland vegetation cover were diminishing and agricultural land,wetland, Built-up areas and waterbodies were increasing in areal extent. Agricultural lands, wetland vegetation cover and wetlands decreased from 2003 to 2017, while bare lands and Built-up areas increased in the area during that time period.
Areas with soil erosion rates 1-5 Mg ha −1 year −1 have diminished in size and areas with rates 1 < and > 5 Mg ha −1 year −1 have increased in extent. These changes are a consequence of land abandonment which eventually contributes to increasing vegetation cover and reduced soil erosion and runoff . Studies in eastern Spain under specific climatic conditions (300 and 500 mm of precipitation year −1 ) showed that semiarid Mediterranean landscapes would respond to abandonment with low vegetation recovery rates and high erosion. In contrast, vegetation in wetter parts of the Mediterranean recovers more quickly after abandonment and erosion rates remain low . Vegetation is the critical factor influencing soil erosion. This has been shown in many ecosystems, but primarily in semi-arid landscapes and agriculture systems where water and sediment delivery are very active due to high connectivity (Cerdà et al. 2018a;. LULC change has increased soil erosion in the study area. The region has been affected by drought, improper water management in upstream watersheds, dam construction, inter-basin water transfers, and a lack of assigned water rights to downstream lands. LULC changes, however, immediately increase erosion and can trigger desertification. Research has shown this in other regions of the world (Ahmadaali et al. 2021;Ghorbani et al. 2021;Fu et al. 2011). Studying the effects of land-use change on erosion trends, a similar study undertaken in the Yezat Basin of northwestern Ethiopia showed that vegetation cover decreased 91% between 2001 and 2010 and then increased 88% between 2010 and 2015 due to the implementation of a comprehensive water management program (Tadesse et al. 2017).
Accelerated erosion destroys agricultural soils, degrades soils' productivity, and contaminates water bodies with sediment. In this study, Landsat imagery was used to assess LULC change and model RUSLE to calculate soil erosion. Erosion seems inevitable and natural, but it should not exceed acceptable levels. Erosion is usually deemed acceptable when it does not exceed the rate of soil formation. This principle is a feature of the United Nations Goals for Sustainability (reviewed by Keesstra et al. 2016).
Based on the available information and considering all factors, the mean rate of soil formation in the study area is approximately 1.1-2 Mg ha −1 year −1 . This should be the limit of the rate of acceptable erosion, so approximately 18.39%, 21.66%, and 19.35% of the area experienced unacceptable levels of soil erosion in 1989, 2003, and 2017. According to the results, the areal extent of erosion at a rate > 1.1-5 Mg ha −1 year −1 has been increased by about 45.56%, 50. 06% and 52.24between 1989, 2003and 2017. This can be explained by climate factors, short-term weather phenomena (like drought), anthropogenic disruption of natural systems, and improper management (Savari et al. 2021). Climate factors are usually natural phenomena that occur over time, but because they reduce vegetation cover, droughts also reduce rainfall erosivity. Disruptions include dam construction, inter-basin water transfers, and land-use changes. The impacts of human activities can be minimized with proper management, consideration of environmental issues, adequate allocation of water rights, comprehensive conservation planning, and sustainable development. This study highlights the need to address these issues by presenting solutions that support management and conservation goals.
One of the main issues this study clarified is that LULC changes are crucial to explaining sediment production in a catchment. LULC changes can lead to land degradation when exceeding sustainable rates of soil losses. Erosion must be controlled over the next decade to maintain soil quality and to prevent land degradation. New policies are needed. A new approach that could be supported is one that improves vegetative cover both in fields under production and on the bare lands. Bare lands should be covered with straw to decrease water and soil losses and to support increased vegetative growth. The government should encourage using the catch crops, mulches, weeds, pruned and chipped branches, and even geotextiles on agricultural lands. Kirchhoff et al. (2017) found that organic farming contributes to soil and plant recovery and controls soil erosion. Cerdà et al. (2018a; demonstrated that erosion could be reduced if the soil surfaces below the canopies on citrus plantations' were covered with chipped or pruned branches. Cerdà et al. (2018a; also used catch crops and weeds in similar settings of the Canyoles River watershed in Spain to achieve the same soil protection. The results show that the trend of medium and high erosion classes is increasing, which is different in different land use. Other studies including Eskandari Damaneh et al. (2021), Senatore et al. (2019) have shown that the southern regions of Iran, especially the south of Khuzestan are affected by climate change such as reduced rainfall and increased temperature. On the other hand by performing human activities such as dam construction and transfer of surface water resources from this area to other area, the trend of land use changes is increased, including reducing the area of waterbodied, which leads to an increase in barren lands and a decrease in agricultural lands.

3 5 Conclusions
Changes to LULCs were determined to be important soil erosion promoters. Rates of change of LULC and change of soil erosion rates due to the soil and vegetation properties were connected to specific land use types. Maps of LULCs were created with supervised classification and maximum likelihood methods using pre-processing of TM, ETM, and OLI images from 1989. RUSLE was used to calculate soil erosion rates. One of its advantages is that it can be used to in regions that lack sufficient, comprehensive data. With the development and advancement of remote sensing science and GIS, sufficient up-to-date and complete databases can be developed to calculate erosion. This study estimated erosion under the changes observed in satellite imagery between 1989 and 2017. The results showed that erosion at rates 2.1-5 Mg ha −1 year −1 increased between 1989 and 2017 due to increased areal extents of both barren and agricultural lands. Based on the moderate and high erosion rate, approximately 45.56% of the area experienced unacceptable levels of soil erosion in 1989. In 2003, 50.06% was eroding faster than the threshold acceptable value. And by 2017, this area grew to 50.24%. According to the results, the areal extent of soil erosion has been increased. And LULC change has increased soil erosion in the Case Study. This study highlights the need to plan and manage LULC changes to diminish erosion to levels that are sustainable. Nature-based solutions should be applied throughout the study region to reduce soil losses.