Analysis of the Impact of Land Use Changes on Soil Erosion Intensity and Sediment Yield Using the IntErO Model in the Talar Watershed of Iran

: Land use change is known as one of the main inﬂuencing factors on soil erosion and sediment production processes. The objective of the article is to study on how land use change impacts on soil erosion by using Intensity of Erosion and Outﬂow (IntErO) as a process-oriented soil erosion model. The study has been conducted under land use changes within the period of 1991–2014 in the Talar watershed located in northern Iran. The GIS environment was used to prepare the required maps including Digital Elevation Model (DEM), geology, land use, soil, and drainage network. The climatology data including average annual precipitation and air temperature as well as the volume of torrential rain were extracted from the data of meteorological stations located inside and around the study watershed. The results indicates that, within the period of 1991–2014, the forest area decreased by 12,478.04 ha (6%), while the other land uses including rainfed agriculture, rangeland, irrigated agriculture, and residential area increased by 7248.25, 4481.05, 476.00, and 273.95 ha, respectively. The estimated outﬂow with 100 year return interval was 432.14 m 3 s − 1 in 1991, which increased to 446.91 m 3 s − 1 in 2014. It can be concluded that the probability of larger and/or more frequent ﬂoods waves in the Talar River is expected to increase. In addition, the amount of production of erosion material (gross erosion) in the watershed increased from 1,918,186 to 2,183,558 m 3 yr − 1 , and the real soil losses per year (sediment yield) of the watershed increased from 440,482.4 to 501,421.3 m 3 yr − 1 . The results clearly emphasized how the lack of appropriate land management and planning leads to increase the maximum ﬂow discharge and sediment yield of the watershed.


Introduction
Soil erosion is one of the most significant causes of land degradation and an important environmental hazard throughout the world, especially in developing countries [1]. Sediment yield and soil erosion are two main constraints on sustainable management of water resources and soil [2]. The quantification of these processes is crucial to design any scientifically based soil and water conservation plan and integrated land management [3,4]. The acceleration of soil erosion due to human activities on a global scale has led to an increased sediment flow in many parts of the world [5]. Unwanted complementary effects of soil erosion, such as loss of soil fertility, reduced water quality, alteration of the hydrological systems, and environmental contaminations, have been identified as a serious problem for human sustainability [6,7]. Land cover and land use changes are key factor in controlling the hydrological response of a watershed. Many studies have shown that there is a significant relationship between land use change and soil erosion [8,9]. Land use change may result in an increase of sediment and nutrient supply to rivers and may affect the water balance in the watershed and its variability, which must be assessed on a local scale [10]. Direct measurement of soil erosion in watersheds and water-sediment sampling is very time consuming and costly [11]. Therefore, the use of soil erosion and sediment yield models at watershed scale is globally raising the interest of specialists. The quantitative understanding of hydrological process at watershed scale also needs the modelling of microscale processes, such as infiltration, permeability and even water and particles transport processes in porous soils [12,13]. Many models, such as the Water Erosion Prediction Project (WEEP), Chemicals, Runoff, and Erosion from Agricultural Management Systems (CREAMS), European Soil Erosion Model (EuroSEM), and Soil and Water Assessment Tool (SWAT), have been developed with varying degrees of complexity in order to fulfil the growing request for a reliable and easy to manage tool to predict erosion and sediment yield. Examples of comparison between the Erosion potential method (EPM) and the Revised Universal Soil Loss Equation (RUSLE) is well presented in the research of Efthimiou et al. [14], and storm event interaction with sediment transport is presented in the research of Panagoulia et al. [15].
The main problem of the process-based models is the large number of input parameters and the lack of data to validate the model predictions [16]. Therefore, empirical models for soil erosion assessment play an important role in soil conservation planning [17]. The EPM is empirical model originally developed for Yugoslavia [18,19] and used in many studies [20][21][22][23][24][25], especially to investigate the effect of land use on soil erosion and sediment yield [26,27]. Applicability of the EPM method in analyzing erosion potential using spatial data manipulation techniques (GIS environment) was first tested in the research of Globevnik et al. [28].
The Intensity of Erosion and Outflow (IntErO) of Spalevic [29] is a program package with the EPM integrated into the algorithm for Windows Operating System [11,30,31].
In some scenarios, climate projections indicate an increase in temperature of as much as 2 • C and a 20% decrease in daily precipitation. These predictions urge experts, especially in the Middle Eastern countries, to assess future changes in water distribution capacities in watersheds and to reconsider the management of reservoirs and dams. It is necessary to work in a broader context to develop adaptation strategies for the optimal use of the available water [11].
The main purpose of this study is the application of the IntErO model to evaluate the effects of land use changes for two periods (1991-2014) on soil erosion and sediment yield in Talar watershed in Iran. With this study we try to create one of the sustainable forms of modelling that would be calibrated and validated in the close region of the studied Talar catchment and afterwards used to evaluate how water resources can meet requirements in water for local community, eco sectors, and agricultural production.

Study Area
The Talar River Basin is located in northern Iran, in the Mazadaran Province, south of the Caspian Sea. Talar is a mountainous watershed with an area of 2055.75 km 2 and subjected to a Mediterranean rainfall regime. The study area has an annual minimum and maximum temperature of 7.7 and 21.1 • C and an average annual rainfall of 552.7 mm, respectively [32]. The location of the Talar catchment is presented in Figure 1. In terms of lithology, the watershed is mostly covered by igneous, metamorphic, and sedimentary rocks [32]. The greatest amount of the geological units (about 61%) consists of Mesozoic rocks. In the study area, 58 types of soil with different physical and chemical properties are present. A, B, and C hydrological soil groups cover 60, 21, 19% of the watershed, respectively, whereas the hydrological group D is not present. The average discharge measured at the Shirgah flow gauge in the 1971-1998 interval is 7.95 m 3 s −1 , whereas the highest discharge ever recorded is 93.46 m 3 s −1 . The most important land uses of the study area are forests, irrigated lands, rainfed agriculture, rangelands, and residential areas.

IntErO Model Application
In Iran, the Intensity of Erosion and Outflow-IntErO program package of Spalevic [29], based on Erosion Potential Method-EPM of Gavrilovic [18,19], proved to return good predictions of runoff and soil erosion intensity, specifically in the Region of Shirindareh River basin, north-eastern Iran [21,[33][34][35][36][37] as well as Khamsan Representative Watershed, western Iran [33]; Chamgardalan Watershed of Ilam Province, western Iran [38]. In this study the IntErO was used to obtain estimates of soil erosion intensity in the Talar river basin, Mazandaran Province, Iran.
The method has clearly defined procedures and subjective evaluations are reduced to a minimum [11,31].
The IntErO model, an upgrading of the Surface and Distance Measuring [39] and River Basins [39,40] programs and can be used for handling a large number of data with the processing of 27 inputs, returning, after the calculations, 22 final result parameters [11]. The IntErO flowchart is presented in Figure 2.
The input data for the IntErO model processing are presented in Table 1 and Figure 3.   The annual volume of soil detached due to soil erosion (W year ) was calculated in m 3 km −2 y −1 by the following Equation (1) [11,19,29,41]: P is the annual average rainfall in mm, T is the temperature coefficient, calculated from the following equation [11,29,41]: t represents the mean annual temperature in • C which is extracted from the meteorological data collected from the Meteorological Organization of Iran.
The Z coefficient describes the intensity of the erosive process. It is calculated using the following equation [11,19,29,41]: X represents the coefficient of soil protection (Table 2), it is a dimensionless parameter based on the vegetal cover and catchment's land use, and it varies from 0.05 to 1; values close to 0 indicate low soil protection while values close to 1 mean high soil protection [42].
Y is soil erodibility coefficient (Table 2), which depends on the pedological and lithological characteristics of the watershed. It indicates the resistance of lithological units and soils to erosion.
The values of this dimensionless factor, which is usually between 0.25 and 2, can be determined either by field measurements or by laboratory experiments [40,42]. Lower values indicate low erodibility, whereas higher values represent strongly erodible formations [42,43]. In the present study, the pedological and the geological map and data of the study area were used to evaluate the Y factor. Then, the coefficient values were assigned to each soil type based on the EPM guidelines and previous studies [42,[44][45][46][47]. Φ coefficient depends on the active erosion and on the degree of extension of linear erosion and mass movements ( Table 2). It is also a dimensionless factor ranged from 0.1 to 1 [48,49]. Each type of erosion form was accordingly related to a value of Φ based on the guidelines of the EPM method [42,44,47]. Finally, the weighted average value was calculated based on the area of all polygons with various values [31].
The total volume of sediment produced in the different areas of the watershed does not fully reach the main river stem. A portion of it is redeposited at the slope foot or deposited in upstream tributaries alluvial plain; therefore, it is essential to calculate the actual specific sediment production (Gsp) in m 3 km −2 y −1 by the following equation [11,19]: This is possible by multiplying the average annual production of sediments Wyear by a delivery ratio coefficient Ru calculated as Lv represents the length of the main river in km, O is the perimeter of the catchment area in km, and D is the average elevation (in m) of the basin with respect to the closing section.
The Maximum flow discharge, Qmax, was calculated using the following formula: A represents the basin shape coefficient and is computed by the following formula: S1 is the coefficient of water permeability of the area calculated from the following equation: Coefficients fp, fpp, and fo are the parts of the river basin that consist rocks of high, medium, and low permeability, respectively. S2 is vegetation cover coefficient computed from the following equation: Coefficients fs, ft, and fg are the parts of the river basin under forest (fs), grass, meadows, pastureland, and orchards (ft), and bare land, plough land, and soils without grass vegetation (fg).
W is analytical expression of inflowing water retention and is presented by the following equation: hb is torrential rain volume in meters, g is acceleration due to gravity (m s −2 ), D is mean height difference of the basin in meters. F is basin area in km 2 .

Land Use Changes (1991-2014)
Landsat satellite images for 1991 and 2014 and the maximum likelihood method were used to prepare the land use map for the Talar watershed ( Figure 4). The study area was classified into five classes including forest, rainfed agriculture, irrigated agriculture, residential area, and rangeland. The comparison of land use change between these two periods indicates that the forest decreased by 12,478.04 ha (6%), whereas rainfed agriculture and rangeland increased by 3.49 and 2.18%, respectively (Table 3). Irrigated agriculture and residential remained almost the same (Table 3).

Climatic Characteristics, Geological Structure, and Soil Characteristics of the Studied Area
The average annual temperature and precipitation were calculated 17 • C and 729 mm based on the Gharakheil meteorological station data. The geology information was obtained from the geological map of Iran [50]. The description of type and area (Ha) of the different formations are reported in Table 4. The Dark grey shale and sandstone, called Shemshak formation, underlie the largest area with 936.3 ha. There are three types of permeability class of formation considered in the IntErO model such as poor water permeability rocks (f0) (Class 1), medium permeable rocks (fpp) (Class 2) and very permeable products from rocks (fp) (Class 3) which occur in 0.3%, 88.2%, and 11.4%, respectively, of the watershed area. The geology and soil maps are shown in Figure 5.

Vegetation and Land Use
The land use map shows that in 2014 forest, irrigated land, rainfed agriculture, rangeland and residential areas cover 34.8%, 0.7%, 14.8%, 48.8%, and 0.9%, respectively, of the Talar watershed area. Most of the basin is covered by rangeland and forestland. The coefficient of the river basin planning (Xa) and the coefficient of the vegetation cover (S 2 ) and Numeral equivalents of visible and clearly exposed erosion process (Φ) are calculated and shown in Table 5.

Soil Erosion and Runoff Characteristics
The data presented in Table 6 describe the result of calculation of the IntErO model for the Talar watershed. The coefficients of the river basin form (A), average river basin width (B) and watershed development (m) were calculated to be 0.2, 29.88 km, and 0.63, respectively. The value of peak discharge, with a return interval of 100 years (Q 100 ) and for a land use setup of 2014 resulted 446.91 m 3 s −1 . The (A) symmetry of the river basin was calculated as 0.31, which indicates that there is a possibility for large flood waves to appear in the studied river basin. In other word, the probability of large flood waves is about once per three year. Q 100 , in fact, is about five times larger than the maximum discharge measured in the 1971-1998 interval for which a maximum return interval of 28 years can be optimistically assumed. In the recent decades human activities such as urban development and land use changes have increased flood hazard and watershed vulnerability to rainfalls and rain storms and an increase of peak flows is to be expected [51,52].  Table 6. The outputs of the IntErO for Talar watershed.

Output Variables Amount and Unit (1991) Amount and Unit (2014)
Coefficient of the river basin form 0. The drainage density of the study river basin (G) we calculated as 0.77 km km −2 . The G index indicates that there is a medium density of the hydrographic network in the Talar watershed. The drainage density is an important factor affecting the flood hydrograph and erosion process. This index depends on the amount of flow through the channel and on soil type [52][53][54]. The index of average river basin decline we calculated to be 43.4%. The value of this index shown that in the studied watershed very steep slopes prevail.
The Height of the river basin local erosion base and the coefficient of basin erosion energy were calculated to be 3790 m and 179.16 m, respectively.
According to the IntErO report Z coefficient is calculated to be 0.492. That indicates that the study area belongs to the 3rd destruction category (out of five); with the medium strength of the erosion process, where the surface erosion is predominant.
The production of erosion material in the Talar watershed was calculated to be 2,183,558.33 m 3 year −1 and the coefficient of the deposit retention (sediment delivery ration) resulted 0.23. It means that 23% of the total eroded material reaches the outlet point whereas the remaining 77% is deposited on the basin slopes and within the hydrological drainage system. The detailed report for Talar watershed hydro morphological parameters is shown in Table 6.
Most of the outputs presented are based on the multiplication of the model parameters, except the Coefficient of the Drainage density (D), Average annual temperature (T 0 ), and the Average slope of the basin (I). According to Dragičević et al. [51], most of these parameters are categorized as high-or medium-sensitivity, whereas those in the multiplication form are classified as very high-sensitivity parameters.
The soil erosion and the actual soil losses based on the land use in 1991 and in 2014 are presented at the Figure 6. Taking into consideration the land use changes in the study catchment reported in the Table 3, the most relevant changes in the Talar are a decrease by 6% of forest lands and an increase of about 4% in agricultural land. The maximum outflow (Q 100 ) increase by 14 m 3 s −1 (3.41%) with land use of 2014. In addition, the amount of sediment production (eroded material) in the river basin and real soil loss increases by 265,372 m 3 yr −1 (6%) and 60,938.9 m 3 yr −1 (12%), respectively. These findings are in line with other studies findings that the degradation and reduction of forest lands and the development of agricultural lands is commonly associated with an increase of flow and sediment yield [55][56][57][58][59][60][61][62][63][64].
The best measure for protection of soil erosion processes and torrential floods is prevention. Establishment of this structured IntErO database for the studied Talar river basin will support in the future planning of permanent control of erosion processes, which will be achieved by an integral river basin management system [65].
Findings from this research are in line with Rodrigo-Comino et al. [66] who clearly demonstrate that the soil erosion intensity is highly dependent on the agriculture use. Furthermore, various strategies and research reports have shown the importance of grown cultures rotation to control the soil losses [30]. The key management aspect to control the soil losses is agriculture and this is clear in many crops: vineyards [66], olive [67], and the other cultures, what leads to the decrease in soil erosion intensity because of the vegetation recovery [68].

Conclusions
Mazandaran Province is one of the northern provinces of Iran, which experienced an appreciable land use change throughout the last three decades. Over more than two decades, in the Talar watershed, a large part of the forest area decrease of about 30% due to road and urban development, mine exploration, and construction of factories. In addition, a large amount of rangelands was turned into agricultural as part of rural areas development programs. All these changes have accelerated soil erosion process. In this study, the IntErO model was used to predict the effect of land use change on soil erosion and sediment yield in the Talar watershed. According to the model results, we concluded that an increase in flood intensity has to be expected following the recent land use change. The study showed that the Talar basin belongs to the third destruction category (out of five), i.e., erosion process is medium and erosion type is predominantly surface erosion. The degradation of forest lands, the development of agricultural land and the conversion of rangelands into rainfed farming is leading to an increase of peak flow and sediment yield in the Talar watershed. The IntErO model proved to be a suitable tool for assessing the watershed response to land use change in terms of peak flow, erosion, and sediment yield, especially in countries with data limitation. The projected climate changes in rainfall amounts and intensity, predicted for the Middle East, may exacerbate the effects of the land use change and remarkably increasing the current hazard level.

Acknowledgments:
The Authors are grateful for support of the organizers of the Green room Sessions 2018 International GEA (Geo Eco-Eco Agro) Conference (www.greenrooms.ucg.ac.me, accessed on 22 March 2021) for funding participation of some of the authors to provide the main findings of this research as Oral presentation, with publishing the Abstract of this research in the Conference Book of Abstracts [69].

Conflicts of Interest:
The authors declare no conflict of interest.