CMIP5 climate projections and RUSLE-based soil erosion assessment in the central part of Iran

Soil erosion (SE) and climate change are closely related to environmental challenges that influence human wellbeing. However, the potential impacts of both processes in semi-arid areas are difficult to be predicted because of atmospheric variations and non-sustainable land use management. Thus, models can be employed to estimate the potential effects of different climatic scenarios on environmental and human interactions. In this research, we present a novel study where changes in soil erosion by water in the central part of Iran under current and future climate scenarios are analyzed using the Climate Model Intercomparison Project-5 (CMIP5) under three Representative Concentration Pathway-RCP 2.6, 4.5 and 8.5 scenarios. Results showed that the estimated annual rate of SE in the study area in 2005, 2010, 2015 and 2019 averaged approximately 12.8 t ha−1 y−1. The rangeland areas registered the highest soil erosion values, especially in RCP2.6 and RCP8.5 for 2070 with overall values of 4.25 t ha−1 y−1 and 4.1 t ha−1 y−1, respectively. They were followed by agriculture fields with 1.31 t ha−1 y−1 and 1.33 t ha−1 y−1. The lowest results were located in the residential areas with 0.61 t ha−1 y−1 and 0.63 t ha−1 y−1 in RCP2.6 and RCP8.5 for 2070, respectively. In contrast, RCP4.5 showed that the total soil erosion could experience a decrease in rangelands by − 0.24 t ha−1 y−1 (2050), and − 0.18 t ha−1 y−1 (2070) or a slight increase in the other land uses. We conclude that this study provides new insights for policymakers and stakeholders to develop appropriate strategies to achieve sustainable land resources planning in semi-arid areas that could be affected by future and unforeseen climate change scenarios.


Material and methods
Study area. This study was conducted in an area covering 5833 km 2 in Alborz Province located in central Iran. The area lies between the latitude 35° 31′-36° 21′ N and longitude 50° 10′-51° 30′ E (Fig. 1). The climate of this area is classified as semi-arid bordering to arid 58 . Mean annual rainfall reaches 251 mm and the mean monthly temperature 14.1 °C. During the year, the temperature typically varies from − 2 to 35.2 °C and the precipitation ranges from 1 or 2 mm to 78 mm in the rainy month (https:// weath erspa rk. com). The study area is characterized by a range of land use and land cover categories including rangeland, agricultural land, saline and bare lands.
Data collection and pre-processing. The soil database was elaborated for several years by the Soil Science Department of the University of Tehran. To predict the spatial distribution of soil texture and soil organic carbon of the study area; Decision Tree (DT) and Artificial Neural Network (ANN) models were generated using 70% of data obtained from laboratory analysis of soil samples. Model testing was based on 15% of the data while 15% was used for model validation. As the performance of the DT is better than ANN (Appendix 1 and 2), the output of the DT model was adopted as the main input for calculating the K factor 59 .
Digital elevation models (DEM) and cloud-free Landsat images of the study area were obtained from the United State Geological Survey (USGS) website at (https:// earth explo rer. usgs. gov). The images obtained include Downscaled CMIP5 monthly precipitation parameter was acquired from the WorldClim Data Portal at 1 km resolution (https:// www. world clim. org/). The CMIP5 data from Global Climate Models (GCMs) are available for four representative concentration pathways (RCPs) as released by the twenty-first century in its 5th Assessment Report (IPCC, 2014). The GCM outputs have been downscaled and calibrated, (i.e. bias-corrected using World-Clim v.1.4 as current baseline climate) 60,61 . Previous research has shown that different CC scenarios produced different result in assessments of GCMs performance in Iranian territory 62,63 . In this research, the HadGEM2-ES model was used for the calculation of rainfall erosivity R-factor. For an investigation of the impacts of future CC on SE by water, the output layers of the current climate and projected CC according to the HadGEM2-ES model and several RCP (RCP2. 6, 4.5, and 8.5) were used in our calculation of the R factor. Soil erosion estimation. In this study, RUSLE 31,64 was used for estimating and predict SE. This is one of the universal pioneer methods for SE estimation and modelling 65 . It is recognized as an empirical model limited to calculating rill and inter-rill erosion, without considering gully erosion 21 . SE estimation using RUSLE is based on the following Eq. (1): www.nature.com/scientificreports/ where ϑ represents the annual soil loss (metric tons per hectare per year); R is the rainfall erosivity (megajoule millimeters per hectare per hour per year); K means soil erodibility (metric ton hours per megajoules per millimeter); LS corresponds to the topographic factor (length and steepness-unitless); C is the land cover/ land use factor (unitless); and, P characterizes support/conservation practice (unitless).
Rainfall erosivity R. R factor refers to the kinetic energy of raindrops which could significantly affect the stability of soil aggregates and enhance soil loss [66][67][68][69][70] . In this study, the R factor was calculated using a monthly database approach as following 31,[71][72][73][74] in Eq. (2): where R is a rainfall erosivity factor (MJ mm ha −1 h −1 per year); p i represents monthly rainfall (mm); and, P corresponds to the annual rainfall (mm). The R factor was calculated for two different periods to account for the past, current and future values. The initial value of the R factor was obtained by computing the average values from 1990-2019 (R average (1990-2019) ) and was considered as a representative average result for the past and current time interval. For future climate projection, two different average values were selected for the R factor. The first one was from 2040-2060 (R average (2040-2060) ) and the second one from 2060-2080 (R average (2060-2080) ).
Afterwards, each pixel was assigned a K value in the GIS environment.
Topographic factor LS. The Slope length (L) and steepness (S) play vital roles in SE and reflect the potential contribution of topography in runoff and SE 65 . The LS factor was computed using the following equation 79 Cover management factor (C). The C factor plays a vital role against SE by protecting the soil surface from the direct effect of raindrops, where erosion is significantly correlated with vegetation coverage 72,81,82 . In this study, the C value was generated by applying Eq. 5, which is based on the Normalized Difference Vegetation Index (NDVI) as follows 83 : where ∀ = 2 and γ = 1. Although there are other three approaches for determining the C factor, the remote sensing approach based on NDVI has been widely used 13,84,85 . The average C factor (C x ) of C 2005 ; C 2010 ; C 2015 and C 2019 was calculated and used as a constant input for addressing the impact of climate projection on SE. It is worth to mention here that satellite images were collected in November and April for each target year, then the average NDVI was calculated for both images to get a representative image for each target year. This approach was undertaken to overcome the fact that NDVI varies widely throughout the year because it is affected by vegetation growth dynamics.
Support practice factor P. The P factor refers to soil loss from up and downslope tillage under specific supporting practices. For instance, contouring agriculture, strip-cropping and terracing affect the direction of surface runoff and modify flow pattern 86,87 . In this study, the P-factor map was derived from DEM and the appropriate value was assigned to each category of the slope following Morgan 88 (Table 1). The spatial pattern of SE was derived by multiplying all the factor together (pixel-by-pixel) to generate a current and future SE map of the study area. In terms of future SE, LS, K, P factors were considered as constant (similar to the current situation), C factor was calculated as an average of (C 2005 , C 2010 , C 2015 , C 2019 ), while R factor was estimated as an average for two different periods the 2050s and 2070s. Hence, the future erosion model could be express ass follow: The methodology was represented in a flowchart in Fig. 2, and Table 2 shows a summary of the data sources used in this study.

Statistical analysis.
We estimated mean values, standard deviations and mean errors for SE factors and total erosion rates using the Extract Values by Points' tool of ArcGIS 10.5 (https:// www. esri. com/ en-us/ about/ about-esri/ overv iew ). Finally, the correlation between total SE and each respective factor was determined using a correlation matrix.

Results
Factors influencing soil erosion. In this study, the soil erodibility (K) factor widely varied. It ranges from 0.2 t·ha·h·ha −1 ·MJ −1 ·mm −1 to 0.4 t·ha·h·ha −1 ·MJ −1 ·mm −1 with a mean value of 0.25 t·ha·h·ha −1 ·MJ −1 ·mm −1 . This variation appears to be influenced by land use and soil type. For instance, the K factor values range from 0.11 t·ha·h·ha −1 ·MJ −1 ·mm −1 (less resistant to SE) in the northern, eastern and southern parts of the study area, to 0.45 t·ha·h·ha −1 ·MJ −1 ·mm −1 (most resistant to SE) the central, northwest and western parts. Erodibility is particularly high in the cultivated area (0.3-0.44 t ha h ha −1 MJ −1 mm −1 ) but lower in areas with high relief (0.25-0.01 t ha h ha −1 MJ −1 mm −1 ) (Fig. 3a). Figure 3b shows that the mean value of the slope length (LS) factor is 4. In the study area, the LS factor ranges from 0.5 to 8.6. The LS value is higher in the high and dissected escarpments in northern, northwest, and northeastern parts of the study area than in the south and southwest ones that are characterized by gentle slopes and low runoff potential. Although the average P factor value is 2.5, more than 60% of the study area registered a low P factor (Fig. 3c). The high values of the P factor coincide with the physiography and severity of slopes in the study area. www.nature.com/scientificreports/ Rainfall erosivity factor ranged between 84 MJ mm ha −1 h −1 per year in lower-lying terrain areas but increased rapidly to 164 MJ mm ha −1 h −1 per year at higher terrain (Fig. 3d). The annual mean of R-value in the study area was 112 MJ mm ha −1 h −1 per year. Figure 4 shows the normalized difference vegetation index (NDVI) and the land cover management factor in the study area. Variation in NDVI values (Fig. 4a) during the period of this study were low in marked contrasts to the values of cover management distribution factor. Figure    www.nature.com/scientificreports/ prone to soil erosion as this part of the study area accounts for more than 20 t ha −1 y −1 of erosion in contrast to the southern parts > (1 t ha −1 y −1 ). Table 3 indicates the area affected per SE categories (%) under both the current and projected climate change scenarios. The annual soil loss in most parts of the region range between 0.1 and 5 t ha −1 y −1 . Table 4 shows the matrix of statistical correlation between RUSLE criteria and the values of SE in the study area. The table indicates that although slope length and management practices are correlated with SE in the study area, slope length has a greater influence on SE than management practices. This is in marked contrast with R, K and C factors that seem to be more fragile in relation to semi-dry physiographic features in the study area.
Projected soil erosion. Three scenarios of projected R factor for 2040-2060 and 2060-2080 were determined from the fifth phase of the Coupled Model Intercomparison Project (CMIP5) models. A comparison of the baseline projected R factor values calculated from monthly rainfall rates of 40 years (i.e. 1960-2000) with the future values of projected R factor derived from three Representative Concentration Pathways (RCPs) is shown in Fig. 6. The highest values of the R factor (< 150 MJ mm ha −1 h −1 ) are mainly concentrated in the west (a) and northwestern (b and c) part of the study area. A similar pattern is observed in the projected values ( Fig. 6d-f). Figure 7 shows the changes between baseline and projected SE in the study area under three CC scenarios of RCPs. This confirms that the regions with higher values of the R factor are located in the eastern and northeastern regions.
Projected future SE values indicate that there will be a high soil loss (> 5 t ha −1 y −1 ) in the north, northwest, and far southern parts of the study area in three scenarios of RCPs (Fig. 8). These future changes show that the spatial distribution of SE is similar to the baseline values (Fig. 5). This future simulation indicates that those same areas would be subject to accelerate SE if adequate soil conservation strategies are not developed and implemented. Notably, areas of high erosion values (> 5 t ha −1 y −1 ) reach up to 19.7% in the RCP2.6 (2060-2080) and 19.1% in the RCP8.5 (2060-2080) ( Table 3 and Figs. 8 and 9). Spatial differences between RUSLE under different RCPs scenarios and RUSLE-2019 show an accelerated erosion trend in most areas of the study area.
The highest values of SE were mainly located in the northwestern parts for the three RCPs (re-coloured), and especially under the RCP8.5 scenario.
Detecting the land uses prone to SE. Table 5 shows the results of baseline and future SE quantities according to the land cover types in the study area. There is an upward trend in the quantities of SE. Rangeland areas accounted for the highest amount of SE especially in RCP2.6 (2070) and RCP8.5 (2070) with an overall amount of 4.25 t ha −1 y −1 and 4.1 t ha −1 y −1 , respectively. Then, they are followed by agricultural areas with 1.31 t ha −1 y −1 (RCP2.6 in 2070) and 1.33 t ha −1 y −1 (RCP8.5 in 2070). Also, bare land areas were predicted to register up to 0.34 t ha −1 y −1 and 0.35 t ha −1 y −1 in RCP2.6 (2070) and RCP8.5 (2070), respectively. The lowest amount of  Table 6 summarizes the changes in statistical parameters from the baseline for each land cover. SE was projected to increase in the 2070s under both RCP2.6 and RCP8.5. In contrast, predicted SE in RCP 4.5 is expected to decline in rangelands by − 0.24 t ha −1 y −1 (the 2050s), and − 0.18 t ha −1 y −1 (the 2070s) but it would slightly increase for the other land-use types.

Discussion
In the semi-arid regions of Iran, SE by water is one of the most complex environmental problems threatening agricultural fields and, subsequently, human well-being. SE in these areas has been evaluated by several studies dealing with water erosion. However, there is a limited number of investigations that have rarely approached the topic of SE rates prediction according to climate covariates across remote areas 89 . In the current study, the impact of future CC on SE was investigated in the semi-arid central part of Iran featured by fragile and motivating properties for land and biodiversity degradation. We did not consider wind erosion in this study, but it is relevant to highlight that future approaches, should combine both water and wind types 90 . The five RUSLE factors (R, K, C, LS, and P) were extracted utilizing information from field survey and remote sensing sources. Then, these www.nature.com/scientificreports/ thematic raster layers were modelled and merged in the GIS environment to calculate the annual rates of SE and considers the spatial-temporal dimensions of SE in the Central Part of Iran. In this regard, given the integration of improved methods in effectively calculating erosion with recent data sources, the provided SE values by water could indicate an elevated accuracy and objectivity considering others obtained from prior studies ( Table 7). The baseline and downscaled R factor in this study was computed and mapped based on monthly precipitation data obtained by WorldClim v.1.4 and v2.1data Portal at 1 km. resolution. Accurate mapping of baseline R factor values led to improved results in estimating SE, especially in an area with low annual precipitation rates and highly governed by climatic conditions. These results could give new insights, for example, to foresee especially the occurrence of rills and gullies among other SE processes because they are very sensitive to changes in rainfall patterns and human impacts 99,100 .        www.nature.com/scientificreports/ Based on the validated modelling process (DT and ANN) fed by the analysis results of 362 soil samples, the spatial distribution of the K factor was mapped. Moreover, K factor values are improved because of testing two reliable models in calculating soil properties based on data derived from remote sensing and extended field survey. In this context of statistical calibration, the DT model provided strong correlations in calculating the soil characteristic with regression of R 2 above 70% in all measurements. However, there is still a further way to improve this model if we consider recent investigations. For example, in China, Wang et al. 101 confirmed that based on the nonlinear best fitting techniques, K factor prediction by combining Geometric Mean    www.nature.com/scientificreports/ 0.36 t ha h ha −1 MJ −1 mm −168 , in western Iran, was between 0.20 and 0.59 0.22 to 0.36 t ha h ha −1 MJ −1 mm −1106 , while the average K value was 0.13 t·ha·h·ha −1 ·MJ −1 ·mm −1 in northern Turkey 107 . The intense spatial-temporal variation of NDVI values has greatly affected the annual C factor values 108 . In this sense, C factor values in 2005, 2019 remarkably different from those in 2010, 2015, which could mainly attribute to megadrought events that hit the central part of Iran in that years 109 . However, C Factor, as the R factor, is largely sensitive to the areas characterized by a semi-arid environment and human impacts 110 . Consequently, these two factors led to the complicated and accelerated dimensions of SE as the current study showed. Observing our results, also non-agricultural must be considered when this factor changes such as residential areas or bare lands. Karpilo et al. 111 stated that there is little consensus in the erosion-science community about the correct values of the C factor for the effects of various slope-protection materials. Therefore, we recommend that policymakers and stakeholders pay attention to that areas especially where both factors show drastically changes from nowadays to the simulated scenario to avoid irreparable loss of fertility (bare soils) or floods and extreme sediment discharges (urban areas).
LS and P factors were mapped based on the reclassification of the slope raster layer. These two were found to be the most influencing factors for erosion acceleration. These results agree with Panagos et al. 112 who highlighted that this factor is quite obviated. They estimated that the P factor could reduce the risk of SE by 3%, with vegetation cover and stone walls obtaining the largest positive impact. However, these results can vary at different scales. Paying attention to research conducted at the hillslope scale, Rodrigo-Comino et al. 113 estimated in a Mediterranean old clementine plantation for the LS factor using two pre-established algorithms and ISUM (Improved Stock Unearthing Method) that the micro-topographical changes can show frequent irregularities in SE results. The authors observed high differences among the areas predicted at the moment of furrow construction and the moment of data survey with soil mobilization rates of about 56.9 m 3 (8.3 mm yr −1 ) in 19 years for 360 m 2 . Comparing LS and P factor maps with the final map of the RUSLE model explained that with rising length and percentage of the slope of the area, intensity and rate of soil erosion also is increased which is along with the result of Mohammadi and et al. 54 .
Multi-digital SE mapping enables calculation of the annual rate of SE which was reached to more than 20 t ha −1 y −1 . The spatiotemporal variation of the resulting SE indicates that there is a spatial concentration of erosion in the northern, northeast, and northwestern regions. The given results indicate that the northern, northeastern and northwestern regions were the most affected in 2005, 2010, 2015, and 2019, respectively. These areas are characterized by mountainous terrain and steep slopes. Our results visibly confirmed that soil erosion could be easily affected by a different kind of land cover. The land use map (Appendix 3) showed that rangelands are dominated in the northern, northeastern and northwestern regions which has the highest soil erosion values, where LS factor ranges from 4-85 (Fig. 3b). On the contrary, bare land dominates in the gentle slope area (LS factor = 0.5-1), which minimize erosion processes. In this regard, Table 4 showed that the highest correlation (r = 0.77, p < 0.05) between topography (LS) and soil erosion, which confirm the eminent role of topography in developing the soil erosion process in the study area. In light of this rangeland and following agricultural regions showing the highest value of SE explosibility, demanded higher protection and management. Thus, that area should be considered in any future land conservation plan as a high priority considering topographical changes as key drivers of weather changes and SE intensity [114][115][116] . This finding completely confirmed the result of the research for Borrelli et al. 117 , in which they discussed in the high slope areas with rare vegetation the risk of SE is high. However, bare land has shown the least values of SE in the current situation, which located on a gentle slope (central part) in compression with other land use.
To assess the effect of future CC on SE susceptibility, data derived from CMIP5 by three scenarios of CMIP5-RCP were used in calculating future values of the projected R factor. However, the utilized approaches in the present study are consistent with those presented by Yigini and Panagos 118 which assumed that future changes in precipitation values will inevitably lead to a change in SE rates globally. Future values of SE were predicted in the study area according to the regional CC concerning three scenarios of RCPs. These future values indicate that the northern, northwest and northeastern regions are the most sensitive and vulnerable to CC, especially under RCP8.5 which consistent with the pathway that involves huge amounts of greenhouse gas emissions 119 . High SE rates are located mainly along mountainous terrain; hence, it will be the most affected by changes in precipitation patterns for climatic characteristics in semi-arid areas 120 . Changes include an increase in extreme precipitation events across the study area, thus a greater precipitation intensity with increase SE potentials by runoff in the steep slope regions. Within this context, prediction of future SE is highly important to provide policymakers with appropriate tools for developing action plans against different possible soil erosion scenarios. Alewell et al. 121 , stressed the importance of modeling SE on a different scale for soil conservation planning and policy governance. Meanwhile, Borrelli et al. 117 emphasized the importance of adaptation of conservation strategies based on RCP2.6 and RCP8.5 scenarios.
Future projections of CC in this study provide a spatial interpretation of the future SE in light of different scenarios. Despite the accuracy and quality of available results and the possibility of using them in the management of soil erosion, some inputs lead to uncertainty in present simulation outcomes. For example, the R factor values were calculated based on the monthly and yearly averages of precipitation (Eq. 3), based on 3 scenarios of RCPs, are still not certain, which could explain the low correlation between projected future R factor and SE in Table 4. However, the adopted calculation method is a suitable alternative in light of the scarcity of data required to calculate the values of the R factor according to the kinetic raindrop energy approach (basically, rainfall records at 15-30 min time interval). Besides this, there is a great difficulty in predicting future C factor with complete accuracy, because the C factor is complex and highly sensitive to environmental changes as was above-discussed. The current spatial outputs are of sufficient reliability for the K, LS and P factors. Consequently, this study presented serious and reliable spatial scenarios about the future of SE in the study area. www.nature.com/scientificreports/ Considering all the factors, it is obvious that SE in the northern and northeastern parts which are dominated by rangeland, higher precipitation, and mountains (high relief topography) is suffering from severe erosion and also have the highest potential for SE. After the assessment of the research results and justifying them according to the geological factors in the study area, we concluded that geology also plays an important role in soil erosion activation at large scales but is mainly reflected in the form of the susceptibility of soil erosion (K-factor). Igneous formations that cover the north part of the study area are correlated with the minimum susceptibility of K-factor while basaltic and tuff formations, due to the lightweight and high porosity, can contribute more to soil formation processes and therefore to the soil erosion factor. In the foothill which sedimentation is the main characteristic, soil erosion has experienced a low rate too. Two parts of the study area revealed a high susceptibility, one is in the Taleghan area that the signs of massive erosion are obvious, and another one located at the Eshtehard with the evidence of gully erosion and dissolution erosion types.
Remarkably, this investigation verified the findings of other researchers in this part of Iran as well as many other regions of Iran 94,122 . The output of this research can be used to take measures of sustainable agriculture in an arid study area environment and to work on identifying priorities for spatial conservation. Also, SE could be mitigated by maintaining vegetation cover, using cover crops, reducing soil disturbance by tillage among other measures [123][124][125][126] .

Conclusions
Soil erosion is one of the most pressing environmental issues in light of the accelerating impacts of global climate change. Multisource GIS provides an objective and advanced platform in soil erosion modeling with accurate and reliable results. Spatial correlation between climate change, soil erosion and land cover change using global models, such as RUSLE, can effectively assist in the spatial management process, especially in arid environments. In the present study, the spatial-temporal distribution of potential soil erosion in the central part of Iran was determined based on future climate change scenarios. The experimental RUSLE model was chosen based on the data specificity of the study area. This model provides the possibility to investigate the current and future spatial distributions of soil erosion rate relying on predictive data. The key findings are summarized, as follows: 1. The average R factor was 112 MJ mm ha −1 h −1 per year, P factor was 2.5, and the K value was 0.25 t·ha·h·ha −1 ·MJ −1 ·mm −1 . The C factor was ranged between 0 and 1, while LS was 0.5-8.6. 2. The estimated annual rate of SE is approximately 12.8 t ha −1 y −1 in Central Iran. 3. Projected future SE values indicate that there will be a high soil loss (> 5 t ha −1 y −1 ) in the north, northwest, and far southern parts of the study area in three scenarios of RCPs 4. Rangeland areas registered the highest amount of SE especially in RCP2.6 (2070) and RCP8.5 (2070), followed by agricultural areas. Also, bare land areas were predicted to considerable SE rates 5. The lowest amounts of SE were estimated for the residential areas in RCP2.6 (2070) and RCP8.5 (2070). 6. SE will be increased in the 2070s under both RCP2.6 and RCP8.5. On the contrary, RCP 4.5 showed that the total SE was predicted to be decreased in rangelands and increased slightly under other land use.