Insights for Estimating and Predicting Reservoir Sedimentation Using the RUSLE-SDR Approach: A Case of Darbandikhan Lake Basin, Iraq–Iran

Abstract

Soil loss (SL) and its related sedimentation in mountainous areas affect the lifetime and functionality of dams. Darbandikhan Lake is one example of a dam lake in the Zagros region that was filled in late 1961. Since then, the lake has received a considerable amount of sediments from the upstream area of the basin. Interestingly, a series of dams have been constructed (13 dams), leading to a change in the sedimentation rate arriving at the main reservoir. This motivated us to evaluate a different combination of equations to estimate the Revised Universal Soil Loss Equation (RUSLE), Sediment Delivery Ratio (SDR), and Reservoir Sedimentation (RSed). Sets of Digital Elevation Model (DEM) gathered by the Shuttle Radar Topography Mission (SRTM), Tropical Rainfall Measuring Mission (TRMM), Harmonized World Soil Database (HWSD), AQUA eMODIS NDVI V6 data, in situ surveys by echo-sounding bathymetry, and other ancillary data were employed for this purpose. In this research, to estimate the RSed, five models of the SDR and the two most sensitive factors affecting soil-loss estimation were tested (i.e., rainfall erosivity (R) and cover management factor (C)) to propose a proper RUSLE-SDR model suitable for RSed modeling in mountainous areas. Thereafter, the proper RSed using field measurement of the bathymetric survey in Darbandikhan Lake Basin (DLB) was validated. The results show that six of the ninety scenarios tested have errors <20%. The best scenario out of the ninety is Scenario #18, which has an error of <1%, and its RSed is 0.46458 km³·yr⁻¹. Moreover, this study advises using the Modified Fournier index (MIF) equations to estimate the R factor. Avoiding the combination of the Index of Connectivity (IC) model for calculating SDR and land cover for calculating the C factor to obtain better estimates is highly recommended.

1. Introduction

Erosion, in its two types, i.e., water and wind [1], is one of the major threats to soil worldwide [2]. Water erosion is affected by climate, land-surface topography, lithology, vegetation, and human-induced activities [3]. Pal [4] stated that more than 19.03 million km² (12.78%) of the world’s land suffers from water erosion as a result of human-induced degradation. Soil loss (SL) directly correlates with reservoir sedimentation (RSed). Increasing SL leads to an increase in the RSed and, consequently, essential problems for water resource development, particularly by increasing the siltation and sedimentation of the reservoirs [5,6]. Furthermore, it can contaminate and degrade the river systems, which costs a lot of money and effort to overcome [5]. Globally, ~84% of land degradation is caused by erosion [7]. SL detaches the uppermost fertile topsoil, which has high concentrations of rich organic matter and nutrients. This will negatively affect soil productivity [8].

The phenomenon by which the sediments are of eroding, transporting, and depositing sediments into the reservoirs by streams is called RSed [9]. It is a fundamental concern for dam operation and affects the dam’s lifetime period by decreasing the reservoir storage capacity [10]. Several factors control the RSed, which leads to difficulties in estimating the amount of the RSed deposited in the reservoir [11]. The rapid development in computer applications contributes to performing several models of the SL and the RSed [12]. SL models help to determine the areas that suffer from erosion susceptibility, assess the rate of erosion, and identify the agent reasons, thereby helping to perform land management [13]. These models can be classified into four widely utilized groups, which are empirical, conceptual, physically based, and hybrid models [1,2]. Empirical models are simple, and their data requirements are less than those that are required for the other three groups [1,13]. In general, empirical models are based on statistical observations. While physically based models are based on the conservation of mass concept. The conceptual models are a combination of empirical and physically based models. Moreover, the hybrid models are a mixture of dynamic and empirical soil-erosion evaluation techniques [1].

The Revised Universal Soil Loss Equation (RUSLE) model [14] is the revised version of the Universal Soil Loss Equation (USLE) model [15]. It is one of the widely used empirical models and is recognized by the scientific community [14]. The RUSLE calculates the average annual soil loss worldwide [16] by calculating the result of six factors, which are rainfall erosivity (R), soil erodibility factor (K), slope length factor (L), slope steepness factor (S), cover management factor (C), and support practice parameter (P) [17].

The major deficiency in the RUSLE model is that, if one of the factors is not accurately estimated, the result will show a significant error [18]. Each factor can be calculated in various ways. The R and C factors are the main factors that highly impact exceeding the SL tolerance limits on erosion control measures [16]. Therefore, properly implementing the R and C factors is fundamental not only for the RUSLE but also for its use to estimate the RSed [16], such as the RUSLE-SDR model [12].

The term “RUSLE-SDR” comes from the integration between the RUSLE model and the sediment delivery ratio (SDR) [12]. The SDR is the fraction of gross erosion that is delivered from a specific catchment to a specific outlet in a specific time interval [19]. A huge number of articles used the empirical SDR-area power function to estimate SDR [20], and other studies used a constant number (between 0 and 1) to treat the SDR [20,21]. Borselli et al. [22] suggested a model that depends on the drainage basin’s hydrological and sediment connectivity to calculate the SDR.

Several works have been accomplished in the Zagros region and surrounding areas to estimate SL and RSed in Iraq [23,24,25,26], Turkey [27,28,29,30,31,32,33], and Iran [25,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50]. However, these investigations used different combinations of equations to estimate RUSLE factors, SDR, and RSed. Three of these articles predicted individual factors of RUSLE rather than the estimation of the RUSLE itself [37,47,50]. Almost all of the articles within the Zagros region estimated RUSLE [25,26,27,28,29,30,31,32,33,34,36,40,41,42,43,44,46,48,49]. In comparison, few of them estimated the RSed by using the RUSLE-SDR model [23,24,35,38,39,42,45]. Among all the above articles, only Zare et al. [42] validated their results. They produced one scenario to estimate the RUSLE-SDR model without testing different combinations of equations to estimate the RUSLE-SDR. To this date, choosing “the best” combination of equations to estimate RUSLE factors and RSed constitutes as a major issue, despite the large number of studies conducted worldwide and specifically in the Zagros region.

This study fills this gap by testing different R, C, and SDR models to estimate the RSed, usually neglected by previous studies accounting only for the most suitable areas where natural erosion may occur. Interestingly, few studies accounted for the sediments that are carried out into reservoirs and validated with bathymetry. Therefore, the main aims of this study are (a) to offer a RUSLE-based model proper for modeling SL in rugged lands and mountainous regions, such as Zagros; (b) to propose a suitable SDR model coupled with the SL models to estimate the RSed; and (c) to validate the models result with the estimated RSed from the bathymetric field measured data of Darbandikhan Lake (DL) in the northwest segment of the Zagros range.

This study’s motivation comes in the context of natural hazard mitigation and water management for one of the major dams in Northern Iraq (Kurdistan Region). Estimation and future prediction of sedimentation in DL, which has neither been studied to estimate SL nor to estimate sediment yield, will benefit the performance of the Darbandikhan Dam and reduce the risk of destructive flooding in this high-population region downstream. Outcomes and motivations also have implications for similar environments worldwide

Table 1. List of acronyms and variables mentioned in the main text.

Term	Abbreviations	Term	Abbreviations
C	Cover management	P	Support practice parameter
CRSed	Sedimentation catchment of its reservoir	R	Rainfall erosivity
DL	Darbandikhan Lake	RI	Topographic surface roughness
DLB	Darbandikhan Lake Basin	RSed	Reservoir Sedimentation
DEM	Digital Elevation Model	RUSLE	Revised Universal Soil Loss Equation
HWSD	Harmonized World Soil Database	S	Slope steepness
IC	Index of Connectivity	SD	Standard deviations
IDW	Inverse Distance Weighting	SDR	Sediment Delivery Ratio
K	Soil erodibility	SL	Soil loss
L	Slope length	SRTM	Shuttle Radar Topography Mission
MCM	Million cubic meters	TRMM	Tropical Rainfall Measuring Mission
MIF	Modified Fournier index	USLE	Universal Soil Loss Equation
NDVI	Normalized Difference Vegetation Index	UTM	Universal Transverse Mercator

is provided to explain all acronyms and variables to improve the readability of the manuscript in forthcoming sections.

2. Darbandikhan Basin

The Darbandikhan Lake Basin (DLB) is located in the northeastern part of Iraq/northwestern part of Iran between 45°11′20″E–47°58′43″E longitudes and 34°13′53″N–35°47′20″N latitudes (Figure 1). It covers a total area of 16,463.1 km², and the major part is located in the northwestern part of Iran (13,155.28 km², i.e., 79.91% of the total catchment area), while the minor part (3307.82 km², i.e., 20.09% of the total catchment area) is located in Kurdistan Region, the northeastern part of Iraq (Figure 1).

The basin is located in a mountainous area, the elevation ranges from 450 m to 3351 m, and the major slope is directed towards the southwest of the area. The average annual precipitation for the last 20 years is 680 mm·y⁻¹, with a major occurrence between October and April. The maximum temperature goes up to 45 °C in August, while the minimum temperature drops down to 3 °C in February [51].

The DL was formed after the construction of the Darbandikhan dam was completed, in November 1961 [52]. The DL feeds by the Sirwan (Diyala) River, which is a seventh-order river (according to the classification of [53]) and is located 55 km southeast of Sulaymaniyah city. The Sirwan River flows 270.4 km inside Iran before becoming part of the Iraq–Iran border for 43.1 km; then it flows 30.9 km inside Iraq until reaching the Darbandikhan dam, with a total length of 344.4 km. The maximum, minimum, and average annual discharge of the Sirwan river are 459.27 m³ s⁻¹ (in 1969), 41.86 m³ s⁻¹ (in 2000), and 153.26 m³ s⁻¹ (1931–2004) [54].

From November 1961 to 1978, the dam’s catchment covers 16463.1 km², which is also the sedimentation catchment of its reservoir (CRSed). In 1978, Iran built the first dam (Vahdat dam) within the catchment of the Darbandikhan dam [55]. Vahdat dam leads to a decrease in the CRSed to 15,403.5 km². With the continuation of the dam construction (13 dams) within Iran, the water supply and the CRSed to DL decreased (Figure 2). The present situation shows that the CRSed for DLB is 5965.8 km², representing 36.2% of the original catchment (

Table 2. Variation of the area of the Darbandikhan Dam catchment from 1961 to 2018.

Period	Area of the Sedimentation Catchment for Darbandikhan Dam (km2)	Area of the Catchment %	Event and the Year	Reference of the Event
1961	16,463.1	100	Building Darbandikhan dam	[52]
1978	15,403.5	93.6	Building Vahdat dam	[55]
2004	13,329.8	81.0	Building Gavoshan dam	[56]
2012	12,253.9	74.4	Building Azadi dam	[57]
2013	11,865	72.1	Building Garan and Ziviyeh dam	[57]
2018	5965.8	36.2	Building Hirwa and Daryan dams	[58]

).

3. Materials and Methods

3.1. Materials

Six scenes of the Digital Elevation Model (DEM) gathered by the Shuttle Radar Topography Mission (SRTM) were mosaicked [59]. The DEMs have a 30 m spatial resolution and are used to extract the drainage network, slope gradient, and flow accumulations. Due to the lack of in situ meteorological data, 3B43-V7 of the Tropical Rainfall Measuring Mission (TRMM) data [60] to construct a precipitation map were utilized. The monthly TRMM is with 0.25° × 0.25° spatial resolution [61]. The soil erodibility was determined by using the Harmonized World Soil Database (HWSD), which has a pixel size of ~852 m [62]. The AQUA eMODIS NDVI V6, with a spatial resolution of 250 m from 2003 to 2021, a scaling factor (of 10⁻⁴), and a radiometric resolution of 16-bit, was used to extract the C factor. The average of the C-factor value for the period from 2003 to 2021 has been used to cover the period from 2021 to 2002. The NDVI product can be accessed from the USGS webpage (https://earthexplorer.usgs.gov/ (accessed on 3 October 2022)) and is obtained by applying Equation (1), which is as follows:

$$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=\frac{\mathrm{\rho }\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{\rho }\mathrm{R}}{\mathrm{\rho }\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{\rho }\mathrm{R}}$$

(1)

where $\mathrm{\rho }\mathrm{N}\mathrm{I}\mathrm{R}$ and $\mathrm{\rho }\mathrm{R}$ correspond to the surface reflectance of both the near-infrared (NIR) and red (R) wavebands of the MODIS, respectively.

Each NDVI scene covers 10-day composited datasets [63]; therefore, the approach suggested by Almagro et al. [64] was applied where four scenes per year were selected to cover the four seasons (i.e., 1–10 January, 1–10 April, 1–10 July, and 1–10 October). Such a strategy enabled us to capture of the variability of the natural vegetation.

Yearly global maps of MODIS land cover (MCD12Q1) were used, which was composed of a supervised classification (decision tree) algorithm of eight observation days by MODIS [65]. These data come with ~463 m spatial resolution, Nadir BRDF-Adjusted surface Reflectance [66], and land-surface temperature [67]. The data were downloaded from the main repository (https://lpdaac.usgs.gov/data, accessed on 3 October 2022) and covered the period from 2001 to 2020 with HDF file format, Sinusoidal grid, and multi-classification maps [65]. In this study, the International Geosphere-Biosphere Programme (IGBP) classification map was selected, encompassing 17 major land-cover classes in the scene. However, our study area includes only nine classes of land cover. All the data mentioned above are free of charge, and they were resized to 30 m cell size and reprojected to the Universal Transverse Mercator (UTM) projection/WGS 1984 datum within zone 38N, using the nearest neighbor resampling method to fit with the DEM scenes spatial resolution.

ArcGIS 10.8 software was used [68] to prepare the data, subset, and mosaic and perform the data operations, such as rater calculator, raster conversion, slope gradient, and stream flow accumulation. The connectivity index toolbox [69] was used to estimate the Index of Connectivity (IC). This toolbox works as a plugin within the ArcGIS environment. The drainage network and watershed boundaries were extracted by using TecDEM 2.2, a MATLAB-based toolbox [70]. Finally, r-based scripts were utilized to implement the statistical analysis [71].

3.2. Methods

Among all approaches, RUSLE has been widely used to estimate soil erosion under different conditions because it meets the need better than any other models available [2]. It has a huge number of works in the scientific literature and a large comparability of results, allowing researchers to adapt the model to nearly every type of condition and region of the world [2]. It is being used to predict long-period rates of rill and inter-rill erosion under different management practices around the globe [12,25]. It is a robust tool to estimate water erosion rates [72] and gives estimates on large spatial scales. Therefore, the RUSLE was chosen to estimate the SL in DL. Eighteen scenarios of RUSLE resulted from six R equations, and three models of the C factor, in addition to the LS, K, and P factors, were used to estimate the SL (

Table A1. The eighteen RUSLE scenarios of the R and C factors combinations were examined in Darbandikhan Lake for the period between November 1961 and 2 May 2022.

Scenarios	C-Equation	R-Equation	Scenarios	C-Equation	R-Equation
1	18	3	10	17	6
2	18	4	11	17	7
3	18	5	12	17	8
4	18	6	13	LC	3
5	18	7	14	LC	4
6	18	8	15	LC	5
7	17	3	16	LC	6
8	17	4	17	LC	7
9	17	5	18	LC	8

in Appendix A). The RUSLE was integrated with SDR (RUSLE-SDR) to estimate a 30 m spatial resolution of the RSed map for the DLB. RUSLE is widely employed [73] to assess the RSed, while the SDR is used to estimate the rate of erosion delivered to the lake. The model of Renard et al. [74] was used to calculate the RUSLE (Equation (2)).

$$\mathrm{A}=\mathrm{R}· \mathrm{K} · \mathrm{L}\mathrm{S} · \mathrm{C} · \mathrm{P}$$

(2)

where A is the average annual rate of the soil loss (t·ha−¹·y⁻¹), R is the annual rainfall and runoff erosivity factor (MJ·mm·ha⁻¹·h⁻¹·y⁻¹), K is the soil erodibility factor (t·ha·h·ha⁻¹·MJ⁻¹·mm⁻¹), LS is the slope length and slope steepness factor, C is the cover management factor, and P is the support practice factor. The LS, C, and P factors are dimensionless. The following subsection describes the RUSLE’s factors.

3.2.1. Rainfall and Runoff Erosivity (R Factor)

Erosion is caused by the driving force of rainfall [75]. The R factor represents the effect of precipitation impact on soil erosion [76]. The relationship between the R factor and precipitation have been determined in many regression analyses [31]. It requires precise and ongoing rainfall data [73]. The study area lacks climatic data records because it was a battlefield during the Iran–Iraq war; therefore, monthly TRMM (3B43–V7) data were used. The TRMM data are one of the main types of input data used to assess and compute the R factor [72,77]. For data validation, the TRMM data were compared with the observed metrological data and showed a good correlation with the observed metrological data (Figure 3). Moreover, TRMM data were proved by several researchers, such as [23,24,78,79,80], to be an applicable source of rainfall data for the Zagros region. However, the appropriateness of using TRMM data in the study area was evaluated by comparing 264 months of data with their corresponding data from the observed precipitation dataset from the Sulaymaniyah meteorological station (Figure 3A). This figure shows a strong direct relationship with a significant p-value < 0.05 and a coefficient of determination (R²) of 0.77. The comparison between monthly means of the 22 years (1998–2019) for the TRMM and Sulaymaniyah meteorological station shows a strong direct correlation with R² of 0.94 (Figure 3B). The average annual precipitation of the DLB ranges between 352.78 mm·yr⁻¹ (in the northeast) and 692.85 mm·yr⁻¹ (in the northwest).

The TRMM pixels were vectorized in a points format, which was interpolated with a spatial resolution of 30 m, using the Inverse Distance Weighting (IDW) method. The pixel size was resampled to obtain the exact pixel size of the DEM. The resulting maps were used to estimate the R factor in six different models, which are Equation (3) [12,13,25], Equation (4) [13,23], Equation (5) [34,35,46,48,81], Equation (6) [29,82], Equation (7) [50,83], and Equation (8) [47,84,85].

The relationship between the R factor and precipitation alters widely based on the different climatic zones [31]. Therefore, the six equations (

Table 3. Rainfall and runoff erosivity (R) factor models were tested in this study.

Method	The Article Used within Iran–Iraq–Turkey	Note	Equation
$\mathrm{R}=79+0.363 {\mathrm{P}}_{\mathrm{A}}$	[25]		(3)
$\mathrm{R}=81.5+0.38 {\mathrm{P}}_{\mathrm{A}}$	[23]	340 < PA < 3500 mm	(4)
$\mathrm{R}=\frac{95.77-6.081\mathrm{M}\mathrm{I}\mathrm{F}+0.447\mathrm{M}\mathrm{I}{\mathrm{F}}^2}{17.2}$	[34,35,46,48]	F > 55 mm	(5)
$\mathrm{R}=\left(4.17\mathrm{M}\mathrm{I}\mathrm{F}\right)-152$	[29]		(6)
$\mathrm{R}=0.264{\mathrm{MIF}}^{1.5}$	[50,83]		(7)
$\mathrm{R}=3.7628\mathrm{M}\mathrm{I}\mathrm{F}-3.532$	[47,84,85]		(8)

) were successfully tested in the Zagros countries (i.e., Iraq, Iran, and Turkey), which have similar climatic conditions to precisely select the R factor equation.

Where $\mathrm{R}$ is the runoff erosivity factor in MJ·mm·ha⁻¹·h⁻¹·y⁻¹, ${\mathrm{P}}_{\mathrm{A}}$ is the average annual precipitation in mm, and MIF is the Modified Fournier index (Equation (9)), which was defined by Arnoldus [82]. The mean MIF was 80.79 mm, which is >55 mm.

$$\mathrm{M}\mathrm{I}\mathrm{F}=\frac{{\sum }_{\mathrm{I} = 1}^{12}{\mathrm{P}}_{\mathrm{m}}{}^2}{\mathrm{P}\mathrm{A}}$$

(9)

where PA is the average annual precipitation (mm), Pm is the average monthly precipitation (mm), and R is the rainfall and runoff erosivity factor (MJ·mm·ha⁻¹·h⁻¹·y⁻¹).

3.2.2. Soil Erodibility (K Factor)

The K factor expresses the potential soil vulnerability to erosion by the R factor [86]. In addition to soil texture, coarse fragments, structure, permeability, and organic matter play an effective role in the K factor value, where the increase in organic matter will decrease its susceptibility to separation [87].

The K factor was estimated by using Equations (10) and (11) [15,73]. The widely applied RUSLE was used to extract the K factor, and the HWSD dataset [88,89] was used to obtain soil texture and soil organic carbon information. Soil organic matter could be estimated from soil organic carbon using a conversion factor [90]. Based on the more reliable hypothesis that carbon represents 58% of the soil OM, the conversion factor will be 1.724 [90].

$$\mathrm{K}=\left[\frac{2.1\times {10}^{-4} {\mathrm{M}}^{1.14 }\left(12-\mathrm{O}\mathrm{M}\right)+3.25\left(\mathrm{s}-2\right)+2.5\left(\mathrm{p}-3\right)}{100}\right]\ast 0.1317$$

(10)

$$\mathrm{M}=({\mathrm{m}}_{\mathrm{s}\mathrm{i}\mathrm{l}\mathrm{t}}+{\mathrm{m}}_{\mathrm{v}\mathrm{f}\mathrm{s}})\ast \left(100-{\mathrm{m}}_{\mathrm{c}}\right)$$

(11)

$$\mathrm{O}\mathrm{M}=1.724·\mathrm{O}\mathrm{C}$$

(12)

where mc is clay fraction content%; m_silt silt fraction content%; m_vfs very fine sand fraction content%; OM is the organic matter content%; OC is the soil organic carbon; and s is the soil structure class, and p is the permeability class, and they can be obtained in

Table 4. Soil structure classes derived from the European Soil Database.

Structure Class (s)	Value	Size	Soil Database
Very fine granular	1	1–2 mm	G (good)
Fine granular	2	2–5 mm	N (normal)
Medium or coarse granular	3	5–10 mm	P (poor)
Blocky, platy, or massive	4	N10 mm	H (peaty topsoil)

and

Table 5. Soil permeability classes are estimated from major soil textural classes.

Permeability Class	Value	Texture
Fast and very fast	1	Sand
Moderate fast	2	Loamy sand, sandy loam
Moderate	3	Loam, silty loam
Moderate low	4	Sandy clay loam, clay loam
Slow	5	Silty clay loam, sand clay
Very slow	6	Silty clay, clay

, respectively [15,73].

3.2.3. Slope Length (L Factor) and Slope Steepness (S Factor)

The slope-length factor reflects the impact of the terrain on the SL. Several equations have been suggested to estimate the L factor, such as [14,15,87]. Equations (13)–(15), suggested by Moore and Burch [91], were applied. Accordingly, the slope map (in percent) was extracted and classified into four groups, i.e., <1%, 1–3%, 3–5%, and ≥5%. Each group of the slope has its constant (m) (Equation (15)). Meanwhile, the S factor is calculated by using three slope (in percent) classes, which are <9%, 9–18%, and ≥18% (Equation (16); [92]).

$$\mathrm{L}={(\frac{\mathrm{\lambda }}{22.12848})}^{\mathrm{m}}$$

(13)

$$\mathrm{\lambda }=\mathrm{FA}\ast \mathrm{Ps}$$

(14)

$$\mathrm{m}=\{\begin{array}{c}0.2 \mathrm{\theta }<1\%\\ 0.3 1\%\le \mathrm{\theta }<3\%\\ 0.4 3\%\le \mathrm{\theta }<5\%\\ 0.5 \mathrm{\theta }\ge 5\%\end{array}$$

(15)

$$\mathrm{S}=\{\begin{array}{c}10.8·\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\theta }+0.03 \mathrm{\theta }<9\%\\ 16.8·\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\theta }-0.05 9\%\le \mathrm{\theta }<18\%\\ 21.9·\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\theta }-0.96 \mathrm{\theta }\ge 18\%\end{array}$$

(16)

where $\mathrm{L}$ is the slope length, $\mathrm{S}$ is the slope steepness factor, $\mathrm{\lambda }$ is the horizontal projection of slope length (m), $\mathrm{m}$ is a constant based on the value of slope gradient (Equation (15)), $\mathrm{FA}$ is stream flow accumulation, $\mathrm{Ps}$ is the pixel size of DEM, and $\mathrm{\theta }$ is the slope gradient in percent.

3.2.4. Cover and Management (C Factor)

The C factor is the fraction of SL from an area with specific vegetation to the corresponding SL from a continuous fallow area [15]. It is one of the changeable erosion factors affected by human action [86]. In this study, three models, which have been used in the Zagros countries, were used to estimate the C factor. Equation (17) is the most common model [93], which has been tested in most of the articles surrounding the DLB, such as Iran [48,49], Turkey [27,33,94], and Iraq [23,26]. The second model to estimate the C factor used the land-cover maps. The land-cover model is applied in Turkey [28,29,31,32], Iran [25,46], and Iraq [25]. Nine land-cover classes of MCD12Q1 MODIS data exist in the DLB, which have different C-factor values (

Table 6. C factors weights for different land-cover classes in the DLB (Model 2).

Name	C Factor	References
Open Shrublands	0.10	[102]
Savannas	0.05	[102]
Grasslands	0.01	[102]
Permanent Wetlands	0	[13]
Croplands	0.3	[12,13,102]
Urban and Built-Up Lands	0	[13,102]
Cropland/Natural Vegetation Mosaics	0.3	[12,13,102]
Barren	0	[13,102]
Water Bodies	0	[12,13]

). The MCD12Q1 MODIS was used to estimate the C factor [77]. The third model was applied by using Equation (18) [95]. This model has been tested in Zagros countries, as well [35,44]. Equations (17) and (18) depend on MODIS NDVI to estimate the C factor, which has been widely applied [96,97,98,99,100].

$$\mathrm{C}=\mathrm{E}\mathrm{x}\mathrm{p}\left[-\mathrm{\alpha }·\frac{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}{\left(\mathrm{\beta }-\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\right)}\right]$$

(17)

$$\mathrm{C}=\frac{-\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}+1}{2}$$

(18)

where $\mathrm{\alpha }$ and $\mathrm{\beta }$ are constants, which are 2 and 1, respectively [101].

3.2.5. Support Practice (P Factor)

The P factor represents the fraction of the SL by an area with specific conservation soil to the corresponding loss upslope and downslope [14]. In DLB, Equation (19) was used. This equation was modified by Othman et al. [23] from the Wener Equation [103,104] to estimate the P factor.

$$\mathrm{P}=0.2+0.03·\mathrm{\theta }$$

(19)

where $\mathrm{P}$ is the support practice factor, and $\mathrm{\theta }$ is the slope angle in percent.

3.3. Sediment Delivery Ratio (SDR)

The SDR is the amount of the SL that actually reaches the specific outline at a specific time [6,20,105]. The range of the SDR is between 0 and 1 [106]. The SDR has been considered for a long time to be a constant number [21]; however, several models have been suggested to estimate the SDR. Most famous models are based on the nonlinear regression between the SDR and the basin area (Equation (20)) [105,107]. The area of the basin is the most affecting factor in determining the SDR for these models. This study considered the suggested models that were only used for estimating the SDR in basins that have an area close to that of DLB.

$$\mathrm{S}\mathrm{D}\mathrm{R}=\mathrm{\alpha }·{\mathrm{A}}_{\mathrm{b}}{}^{-\mathrm{\beta }}$$

(20)

where the SDR is the sediment delivery ratio, the α and β are coefficients (Table 7), and the A_b is the basin area in km².

Borselli et al. [22] suggested a new approach to estimating the SDR (SDR1; Equation (21)), which depends on calculating the IC (Equation (22)). The IC depends on the topographical information and can be calculated in an ArcGIS environment [22]. Topographic surface roughness (RI) was used to estimate the average weighting factor $\left(\overline{\mathrm{W}}\right)$ (Equation (23)), which was suggested by [69]. The RI can be estimated by using Equation (24) [108].

Table 7. The α and β coefficients used to estimate the SDR in the DLB.

α	β	References	Unit of the Area	Model No.
0.4724	0.125	[32,94,105]	km2	SDR2
1.817	0.132	[23,107]	km2	SRD3
2.945	0.205	[107]	km2	SDR4
0.51	0.11	[77,109,110]	mi2	SDR5

Table 7. The α and β coefficients used to estimate the SDR in the DLB.

α	β	References	Unit of the Area	Model No.
0.4724	0.125	[32,94,105]	km2	SDR2
1.817	0.132	[23,107]	km2	SRD3
2.945	0.205	[107]	km2	SDR4
0.51	0.11	[77,109,110]	mi2	SDR5

$$\mathrm{S}\mathrm{D}\mathrm{R}=\frac{\mathrm{S}\mathrm{D}{\mathrm{R}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{1+{\mathrm{e}}^{\left(\frac{\mathrm{I}{\mathrm{C}}_0-\mathrm{I}{\mathrm{C}}_{\mathrm{i}}}{{\mathrm{k}}_{\mathrm{b}}}\right)}}$$

(21)

where SDR is the sediment delivery ratio; SDR_max is the maximum theoretical sediment delivery ratio, which is suggested to be 0.8 [6,21]; and IC₀ and K_b are calibration parameters and equal to 0.5 and 2, respectively [12]. Both parameters can be determined by the S-shape of the sigmoid function relationship between the sediment delivery ratio and the IC [6]. The IC_i is the index of connectivity for a specific cell.

$$\mathrm{I}\mathrm{C}=\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}\left(\frac{{\mathrm{D}}_{\mathrm{u}\mathrm{p}}}{{\mathrm{D}}_{\mathrm{d}\mathrm{n}}}\right)=\mathrm{l}\mathrm{o}{\mathrm{g}}_{10}\left(\frac{\overline{\mathrm{W}}\overline{\mathrm{S}}\sqrt{\mathrm{A}}}{{\sum }_{\mathrm{i}}\frac{{\mathrm{d}}_{\mathrm{i}}}{{\mathrm{W}}_{\mathrm{i}}{\mathrm{S}}_{\mathrm{i}}}}\right)$$

(22)

where IC is the Index of Connectivity; Dup is the upslope components; D_dn is the downslope components; $\overline{\mathrm{W}}$ is the average weighting factor of the contributing area; $\overline{\mathrm{S}}$ is the average slope of the upslope contributing area (m/m); A is the area of the contribution (m²); and d_i, W_i, and S_i, are the flow path to the main downstream channel, the weighting factor, and the slope gradient of the ith cell, respectively.

$${\mathrm{W}}_{\mathrm{i}}=1-(\frac{\mathrm{R}\mathrm{I}}{\mathrm{R}{\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}})$$

(23)

$$\mathrm{R}\mathrm{I}=\sqrt{\frac{{\sum }_{\mathrm{i} = 1}^{\mathrm{m}\mathrm{w}{\mathrm{s}}^2}{({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}})}^2}{\mathrm{m}\mathrm{w}{\mathrm{s}}^2}}$$

(24)

where W_i is the weighting factor, RI is the topographic surface roughness, mws is the moving window size (normally 5-cells; [108]), xi is the one specific cell value of the residual topography within the moving window, and $\overline{\mathrm{x}}$ is the mean of the cells values within the moving window.

3.4. Reservoir Sedimentation (RSed)

Equation (25) [15] was used to estimate the RSed reached to DL in t·ha⁻¹·yr⁻¹. The absolute RSed for the six stages of the DLB areas were calculated in tons (Figure 2), where the DLB was changed with time as a result of the dams constructed within the basin (Table 2). Moreover, the amount of the RSed in tons was estimated until 2 May 2008, to validate the result of each scenario (see Section 4.2)

$$\mathrm{R}\mathrm{S}\mathrm{e}\mathrm{d}=\mathrm{R}\mathrm{U}\mathrm{S}\mathrm{L}\mathrm{E}·\mathrm{S}\mathrm{D}\mathrm{R}$$

(25)

where RSed is reservoir sedimentation, RUSLE is the soil loss, and SDR is the sediment delivery ratio.

3.5. Validation

For any model, validation is the most significant procedure to check the accuracy of the results [111]. The information about the RSed volume within the DL (

Table A2. The amount of the reservoir sedimentation gained for the period between November 1961 and 2 May 2008, within the Darbandikhan [112].

Elevation	VOLUME DIFF (MCM)	Elevation	VOLUME DIFF (MCM)	Elevation	VOLUME DIFF (MCM)	Elevation	VOLUME DIFF (MCM)
434	186.5196177	449	368.2601237	464	454.2837	479	444.557
435	202.4116618	450	376.7507805	465	460.0191	480	454.4112
436	217.8700461	451	383.6461686	466	462.8374	481	446.2946
437	232.2080715	452	390.8162663	467	464.5166	482	434.6787
438	246.3789168	453	397.2368461	468	465.7219	483	435.0195
439	259.793523	454	402.6708382	469	467.6152	484	434.0626
440	273.5651565	455	408.1217914	470	468.8454	485	423.8966
441	285.4514375	456	412.3380251	471	469.2708	486	404.414
442	297.1203966	457	416.9450171	472	469.4618	487	368.7144
443	308.6056451	458	422.0896716	473	471.0426	488	323.7879
444	319.7122029	459	427.3771193	474	472.0941	489	267.4594
445	330.0812299	460	433.8615028	475	474.5954	490	209.1261
446	340.4472993	461	438.4430503	476	471.8698	491	145.1009
447	349.8499081	462	443.3719036	477	467.5566	492	59.9097
448	359.4951166	463	448.2591828	478	451.6446	493	0

in Appendix B) offered by the ELC, Electroconsult; MED, Ingegneria; and SGI, Studio Galli Ingegneria companies [112] was used to validate the results of the 90 estimated scenarios of the RSed. Reference [112] used the historical topographic maps [113] produced before building the dam and a recent sub-water topographic map created by the topo-bathymetric survey by echo-sounding bathymetry data for the DL to estimate the amount of the RSed for the period between November 1961 and 2 May 2008.

The survey for the DL was implemented by [112] in clear, calm weather and water conditions. It was carried out between 25 April to 2 May 2008, at the dam sites, utilizing two small boats; an Echo sounder (Single beam Sonar Bathy500 Dual Frequency-200 kHz (10 degrees) Transducer) rod, and GPS (Positioning Trimble R6) antenna were mounted to the vessel via a wooden board manufactured directly at the site [112].

To convert the weight of the 90 scenarios (in tons) to size (in km³), three samples from the siltation within the DL were collected to find the average density of the siltation. Since the lake’s water level dropped as a result of the reduction in rainfall, the sediments exposed in the lake became a typical location to sample the actual sedimentation of the lake. The scenarios estimated the RSed in tons, while the siltation measured by the [112] is estimated in a million cubic meters (MCM).

All the results of the RUSLE factors, RUSLE, SDR, and RSed were presented for the periods between November 1961 to 2 May 2008. Finally, the RSed within the DL will be predicted by the end of 2019 using the best scenario out of the ninety scenarios tested. Moreover, the error rates were calculated using Equation (26):

$$\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}\%=|\left(\frac{\mathrm{i}\mathrm{n} \mathrm{s}\mathrm{i}\mathrm{t}\mathrm{u}-\mathrm{R}\mathrm{S}\mathrm{e}\mathrm{d}}{\mathrm{i}\mathrm{n} \mathrm{s}\mathrm{i}\mathrm{t}\mathrm{u}}\right)·100|$$

(26)

where in situ means the results obtained by the surveys by echo-sounding bathymetry, and RSed is the outcome of the erosion model.

In addition, the only available survey for the DL, carried out between 25 April and 2 May 2008, is used indirectly to validate the estimation of the C and R factors, and the best scenario resulted from the best RUSLE and SDR combination for estimating the RSed.

The error density plots display the relationship between the position and dispersion of the estimations for the factors (R and C). For this target, 984,988 pixels (~5% of the total pixels) were randomly selected and utilized. This test was repeated more than one time. The uncertainty plots display two standard deviations (SDs) of these estimates in the y-axis against their mean value in the x-axis, using the chosen pixels. The uncertainty plots allow the interpreter to determine the locations of the dispersion values for the estimated scenarios [80].

4. Results

The ninety scenarios of the RSed models for the DLB were estimated until 2 May 2008 and compared with the in situ surveys by echo-sounding bathymetry. The best RSed scenario resulting from this comparison was Scenario #18 (Table A2 in Appendix B). This scenario used Equation (5) to estimate the R factor, Equation (18) to estimate the C factor, and the IC model to estimate the SDR. Therefore, Equations (5) and (18) were reported to exhibit the result of the R factor and C factor, respectively, instead of other equations. In this section, a brief of all scenarios that were tested is given, and the focus will be on the optimum scenario result. The following subsections show the results of calculating RSed, SDR, RUSLE, and its factors.

4.1. Estimation RUSLE and Its Factors

Table 8. R-factor parameters of the used models.

R Factor	Minimum	Maximum	Mean	SD
Equation (3)	215.80	332.54	290.14	25.44
Equation (4)	224.71	346.92	302.53	26.64
Equation (5)	83.69	335.47	210.40	64.18
Equation (6)	106.75	347.05	242.71	60.30
Equation (7)	129.04	345.63	245.28	54.81
Equation (8)	229.95	446.79	352.64	54.39

shows six models of the R factor. The R factor resulting from Equation (5) varies from 83.69 MJ·mm·ha⁻¹·h⁻¹·y⁻¹ to 335.47 MJ·mm·ha⁻¹·h⁻¹·y⁻¹, with an average value of 210.4 MJ·mm·ha⁻¹·h⁻¹·y⁻¹. The R-factor map was sliced into five classes, which decrease toward the east. The areas with very high class (>275 MJ·mm·ha⁻¹·h⁻¹·y⁻¹) are located in the northwest, which is almost all located within Iraqi areas. Meanwhile, the regions with very low class are located in the eastern part of the DLB, within Iranian areas (Figure 4A). Figure 5 shows the distributions of the ~5% random selected pixels from R-factor values for the six equations. The R-factor distribution in Equation (8) has a higher performance than the others, while the performance of the R factors for Equations (5), (7), and (6) looks the same (Figure 5). Equations (3) and (4) have intermediate behavior between the two groups (i.e., Equations (5), (7), and (6); and Equations (3) and (4)).

Similarly, four types of soil were exposed in the study area. The majority of the DLB is loam texture, followed by clay and clay loam (

Table 9. Soil types, textures, and the K-factor values.

Soil Type	Texture Class	Sand%	Silt%	Clay%	K Factor
Lithosols	Loam	43	34	23	0.048767
Chromic Vertisols	Clay	16	29	55	0.023007
Haplic Xerosols	Clay loam	23	33	44	0.056780
Calcic Xerosols	Clay loam	40	37	23	0.063365

). Figure 4B is the K-factor map, where the K factor value is low (0.023007 t·ha·h·ha⁻¹·MJ⁻¹·mm⁻¹) in the northwestern part of the DLB, and the K factor value is high (0.063365 t·ha·h·ha⁻¹·MJ⁻¹·mm⁻¹) in small patches in the north and northeastern parts of the DLB (Figure 4B).

The slope gradient in the DLB reaches 77.6°, and the average is 16.2°. The majority of the high-slope areas are located in the central part of the DLB. This variation between the slopes, coupled with the variation in the flow accumulation, is the cause of the high fluctuation of the LS factor values. Likewise, the LS factor was sliced into five classes. Approximately 41.2% of the study area has LS < 1; almost all of these areas are located in the western part. The highest LS values (>100) cover 2.35% of the study area, while the average of the LS factor is 6.84 t·ha⁻¹·y⁻¹ (Figure 6A).

Similar to the R factor, the three calculated models of the C factor are stated in

Table 10. C-factor parameters of the used models.

C Factor	Minimum	Maximum	Mean	SD
Equation (17)	0.029	1	0. 618	0. 13
Land cover	0	0.3	0.091	0.127
Equation (18)	0.213	0.579	0. 396	0. 034

. The C factor ranged between 0.21 and 0.58 for Equation (18), which is classified into five classes (Figure 6B). The very high values (0.45–0.58) are located in the landslide and highly eroded areas (Figure 7A), which are distributed in small patches within the DLB. In general, the eastern part of the DLB has a C factor more than the western part (Figure 6B), while the lower values (<0.3) were presented in the forest and the agricultural areas (Figure 7B). Figure 8 shows the distributions of the ~5% random selected pixels from C-factor values for the three models (i.e., two equations and land-cover based), which have significant differences. The C-factor distribution in Equation (18) is higher than and the land-cover based one, while Equation (17) is higher than both (i.e., Equation (18) and land-cover based; Figure 8).

Similarly, Figure 9A shows the P factor, which is subdivided into five classes and ranges from 0.2 to 0.86. The very high P-factor values are located in the rough topography, which is mostly in the central part of the DLB. At the same time, the northwestern part of the study area shows low P-factor values.

Figure 10 shows the distributions of one set of random selected pixels (~5% of the total data of the study area) from RUSLE model values for the eighteen scenarios, where more than one set was tested. More detailed information about the eighteen scenarios of RUSLE models of the R- and C-factor combinations can be found in Table A1 in Appendix A. It shows that there are two main groups (i.e., Group 1, RUSLE 1 to 12; and Group 2, RUSLE 13 to 18). The SL or RUSLE map was classified into five groups, ranging from 0 t·ha⁻¹·y⁻¹ to 83,628.8 t·ha⁻¹·y⁻¹, with an average value of 8.3 t·ha⁻¹·y⁻¹ (Figure 9B). The areas with very low and low SL (0–3 t·ha⁻¹·y⁻¹) represent >59% of the study area. They can be observed in the northwestern (near Sulaymaniyah city) parts of the DLB, in addition to the eastern parts (south of Muchesh and Sanandaj cities). Severe (12%) and very severe (15%) SL areas are located in the central part of the DLB. The extremely severe SL areas (>1000 t·ha⁻¹·y⁻¹) cover 0.15 of the DLB.

4.2. Sediment Delivery Ratio (DRr), Reservoir Sedimentation (RSed), and the Model Validation

The five SDR models (i.e., the IC model and Table 7) show fluctuation in their results.

Table 11. SDR ranges for the five models used in the DLB for the period between 1961 and 2008.

Model No.	Minimum	Maximum	Mean	SD
1	0.125	0.128	0.126	0.0014
2	0.509	0.519	0.511	0.0059
3	0.402	0.420	0.410	0.0074
4	0.172	0.176	0.174	0.0017
IC (Equation (22))	0.013	0. 147	0.0327	0.0076

shows the statistical parameters of the SDR models. For the best scenario (IC model), the SDR ranges from 0.013 to 0.147, with an average value of 0.0327. The IC model map was classified into five classes (Figure 11A). The very high and high classes areas are located within and near the main valleys. In contrast, the very low and low classes areas can be observed in the northwestern (near Sulaymaniyah city) parts of the DLB, in addition to the eastern parts (south of Muchesh and Sanandaj cities).

The average density of the three samples collected from the field trip (1.631 kg/L) was used to convert the RSed scenarios from tons to km³. Figure 12 shows the distribution of part of the scenarios tested in this study (more than 20 km³·yr⁻¹ and less than 50 km³·yr⁻¹). The amount of Scenario #18, the best scenario (22.294 km³), makes it the closest model to the siltation measured by [112] within the DLB for the period from November 1961 to 2 May 2008, which is 22.223 km³. The average siltation per year is 0.46458 km³·yr⁻¹. This scenario is followed by Scenarios #66, #61, and #62.

Almost all of the study areas are provided by the DLB with the sedimentation of <5 tons/year. Some of the areas, specifically in the central part of the DLB, provided the DLB with sedimentation between 5 and 50 tons/year. Very few areas (could be considered outliers) supplied the DLB with sedimentation >50 tons/year, which is normally located in the main streams (Figure 11B).

4.3. RUSLE, Its Factors, and Reservoir Sedimentation in the Present Day

The RUSLE for the DLB and the RSed within the DL were estimated by the end of 2019, using Scenario #18 (the best scenario out of the ninety scenarios tested). The RUSLE map was classified into five groups, ranging from 0 t·ha⁻¹·y⁻¹ to 82,725.2 t·ha⁻¹·y⁻¹, with an average value of 8.2 t·ha⁻¹·y⁻¹ (Figure 13A). The classes within the RUSLE maps show the same distribution for the old (2008) and new (2019) maps.

The average RSed within the DLB for the period from 2 May 2008 to 2 December 2019 is 0.3836 km³.yr⁻¹ (Figure 13B). The total amount of the RSed within the DL for the period from November 1961 to the end of 2019 is 27.05 km³. The prediction of the RSed to the end of 2022 will be 28.201 km³.

5. Discussion

5.1. RUSLE-SDR and Its Factors

Two main insights must be considered when comparing predictions with in situ data. First, the model is valid to be applied to a specific area, which has its own circumstances. Therefore, almost all of the available Zagros literature was reviewed in the Scopus dataset [23,25,26,27,28,29,31,32,33,34,35,44,46,47,48,49,50,83,84,85,94]. Second, the successful prediction scenario must have an acceptable correlation with the in situ measurement. In this study, the final output of the 90 RSed scenarios was verified by using the bathymetry survey of the DL, which mirrors the verification of RUSLE factors and the SDR. When the bathymetry survey is compared with obtained RSed, Scenarios #18 and #66 give approximately similar results (the best; error < 1%), followed by Scenarios #61 and #62 (better; error < 10%), and then Scenarios #13 and #14 (good; error < 20%).

Among the six equations of the R factor, the equations that used MIF [82] give better results than other models, particularly Equations (5), (7), and (6), best, better, and good, respectively. Scenarios #18 and #66 used Equation (5) to estimate the R factor, while Scenarios #61 and #62 used Equations (7) and (6), respectively. Meanwhile Equation (8), which also depends on MIF, did not show a promising result (

Table A3. The results of the scenarios examined in the Darbandikhan Lake for the period between November 1961 and 2 May 2022.

Scenarios	C Factor	R Factor	SDR	RSed (km3)	Error%	Scenarios	C Factor	R Factor	SDR	RSed (km3)	Error %
1	Equation (17)	Equation (7)	SDR1	36.96916	66.3554	46	Land cover	Equation (3)	SDR5	5.181852	76.68248
2	Equation (17)	Equation (6)	SDR1	37.03418	66.64798	47	Land cover	Equation (4)	SDR5	5.417935	75.62015
3	Equation (17)	Equation (8)	SDR1	52.05412	134.2353	48	Land cover	Equation (5)	SDR5	4.169816	81.23648
4	Equation (17)	Equation (3)	SDR1	42.38818	90.74013	49	Equation (18)	Equation (7)	SDR5	34.57858	55.59816
5	Equation (17)	Equation (4)	SDR1	47.80158	115.0996	50	Equation (18)	Equation (6)	SDR5	34.63592	55.85619
6	Equation (17)	Equation (5)	SDR1	32.41343	45.85533	51	Equation (18)	Equation (8)	SDR5	48.8311	119.7323
7	Land cover	Equation (7)	SDR1	3.045926	86.29381	52	Equation (18)	Equation (3)	SDR5	39.88271	79.46591
8	Land cover	Equation (6)	SDR1	3.062053	86.22124	53	Equation (18)	Equation (4)	SDR5	41.76512	87.93646
9	Land cover	Equation (8)	SDR1	4.221344	81.00462	54	Equation (18)	Equation (5)	SDR5	30.20864	35.93412
10	Land cover	Equation (3)	SDR1	3.415393	84.63127	55	Equation (17)	Equation (7)	SDR4	287.1974	1192.343
11	Land cover	Equation (4)	SDR1	3.575052	83.91283	56	Equation (17)	Equation (6)	SDR4	287.0661	1191.752
12	Land cover	Equation (5)	SDR1	2.718805	87.76581	57	Equation (17)	Equation (8)	SDR4	409.8523	1744.271
13	Equation (18)	Equation (7)	SDR1	26.49843	19.23876	58	Equation (17)	Equation (3)	SDR4	336.9121	1416.051
14	Equation (18)	Equation (6)	SDR1	26.56477	19.53728	59	Equation (17)	Equation (4)	SDR4	352.692	1487.058
15	Equation (18)	Equation (8)	SDR1	37.20219	67.404	60	Equation (17)	Equation (5)	SDR4	248.1158	1016.482
16	Equation (18)	Equation (3)	SDR1	30.24063	36.07807	61	Land cover	Equation (7)	SDR4	24.96615	12.34374
17	Equation (18)	Equation (4)	SDR1	31.67841	42.54786	62	Land cover	Equation (6)	SDR4	25.08741	12.88939
18	Equation (18)	Equation (5)	SDR1	22.29416	0.320209	63	Land cover	Equation (8)	SDR4	34.49158	55.20668
19	Equation (17)	Equation (7)	SDR3	358.5325	1513.34	64	Land cover	Equation (3)	SDR4	27.89489	25.52261
20	Equation (17)	Equation (6)	SDR3	358.3681	1512.6	65	Land cover	Equation (4)	SDR4	29.16574	31.24124
21	Equation (17)	Equation (8)	SDR3	511.6561	2202.372	66	Land cover	Equation (5)	SDR4	22.44697	1.00783
22	Equation (17)	Equation (3)	SDR3	420.5986	1792.627	67	Equation (18)	Equation (7)	SDR4	186.1391	737.5966
23	Equation (17)	Equation (4)	SDR3	440.2978	1881.271	68	Equation (18)	Equation (6)	SDR4	186.4481	738.9871
24	Equation (17)	Equation (5)	SDR3	309.742	1293.79	69	Equation (18)	Equation (8)	SDR4	262.8601	1082.829
25	Land cover	Equation (7)	SDR3	31.16664	40.24497	70	Equation (18)	Equation (3)	SDR4	214.6905	866.0734
26	Land cover	Equation (6)	SDR3	31.31803	40.9262	71	Equation (18)	Equation (4)	SDR4	224.8236	911.6708
27	Land cover	Equation (8)	SDR3	43.05787	93.75363	72	Equation (18)	Equation (5)	SDR4	162.6163	631.7477
28	Land cover	Equation (3)	SDR3	34.82278	56.69703	73	Equation (17)	Equation (7)	SDR2	99.73186	348.7777
29	Land cover	Equation (4)	SDR3	36.40928	63.83603	74	Equation (17)	Equation (6)	SDR2	99.6861	348.5717
30	Land cover	Equation (5)	SDR3	28.02177	26.09355	75	Equation (17)	Equation (8)	SDR2	142.3258	540.4437
31	Equation (18)	Equation (7)	SDR3	232.3718	945.6365	76	Equation (17)	Equation (3)	SDR2	116.9967	426.4667
32	Equation (18)	Equation (6)	SDR3	232.7572	947.3707	77	Equation (17)	Equation (4)	SDR2	122.4763	451.1241
33	Equation (18)	Equation (8)	SDR3	328.1499	1376.623	78	Equation (17)	Equation (5)	SDR2	86.1599	287.706
34	Equation (18)	Equation (3)	SDR3	268.0158	1106.029	79	Land cover	Equation (7)	SDR2	8.669506	60.98859
35	Equation (18)	Equation (4)	SDR3	280.6658	1162.952	80	Land cover	Equation (6)	SDR2	8.711617	60.7991
36	Equation (18)	Equation (5)	SDR3	203.0056	813.4932	81	Land cover	Equation (8)	SDR2	11.97725	46.10426
37	Equation (17)	Equation (7)	SDR5	53.35229	140.0769	82	Land cover	Equation (3)	SDR2	9.686521	56.41218
38	Equation (17)	Equation (6)	SDR5	53.3278	139.9667	83	Land cover	Equation (4)	SDR2	10.12783	54.42636
39	Equation (17)	Equation (8)	SDR5	76.13834	242.6105	84	Land cover	Equation (5)	SDR2	7.794709	64.92504
40	Equation (17)	Equation (3)	SDR5	62.58831	181.6375	85	Equation (18)	Equation (7)	SDR2	64.63809	190.8612
41	Equation (17)	Equation (4)	SDR5	65.5197	194.8283	86	Equation (18)	Equation (6)	SDR2	64.74529	191.3436
42	Equation (17)	Equation (5)	SDR5	46.09182	107.4059	87	Equation (18)	Equation (8)	SDR2	91.28039	310.7474
43	Land cover	Equation (7)	SDR5	4.637795	79.13065	88	Equation (18)	Equation (3)	SDR2	74.55308	235.4771
44	Land cover	Equation (6)	SDR5	4.660323	79.02928	89	Equation (18)	Equation (4)	SDR2	78.07189	251.3112
45	Land cover	Equation (8)	SDR5	6.40729	71.1682	90	Equation (18)	Equation (5)	SDR2	56.46937	154.1033

in Appendix C,

Table 12. The results of the best scenarios, which have error rates less than 20% in the Darbandikhan Lake for the period between November 1961 and May 02, 2022.

Scenarios	C Factor	R Factor	SDR	RSed (km3)	Error %
18	Equation (18)	Equation (5)	SDR1	22.29	0.32
66	Land cover	Equation (5)	SDR4	22.445	1.01
61	Land cover	Equation (7)	SDR4	24.97	12.34
62	Land cover	Equation (6)	SDR4	25.09	12.89
13	Equation (18)	Equation (7)	SDR1	26.50	19.24
14	Equation (18)	Equation (6)	SDR1	26.57	19.54

, and Figure 5). Moreover, the standard deviations of the R factor in the study area are lower than those of Azari et al. [47]. The results of the equations of the R factor are located within the range stated by Azari et al. [47]. Therefore, this study agrees with them to use the MIF models because of the impossibility of extracting the R factors in more suitable equations.

Scenario #18 used Equation (18), while Scenarios #66, #61, and #62 used the land-cover map to estimate the C factor. Scenario #18 used the IC model to estimate the SDR, while Scenarios #66, #61, and #62 used Equation (18), with α = 2.945 and β = 0.205 (Table A3 in Appendix C). This study shows that almost all (68 scenarios) of the prediction models are higher than the in situ measurement conducted by [112]. Eighteen scenarios are less than the in situ measurement.

This study recommends using SDR₁ (IC) because it is changeable and can be applied for different area sizes and depends on the influence of topography on sediment connectivity [69], in contrast to SDR₄, which depends on the area size and is affected by the change in landform. Moreover, the Zagros region is a mountainous area, varies in rock types and strengths, and has high tectonic activities [114], with medium to high erosional risk [115,116,117], confirming that this area suffers from severe erosion with low delivery for the soil. In addition, the results of RUSLE 3 correspond with the results of the research reported in the literature [25,31,32,45,48].

The scenarios of RSed, which have high accuracy, can be classified into two groups: scenarios with high SDR and low RUSLE and scenarios with low SDR and high RUSLE. The best SDRs are SDR₁ and SDR₄ because they comes three times for the best six scenarios (present of 50%), which has an accuracy >80 (Table 12). These two SDRs show variation in the delivery rate, as SDR₁ ranges between 1.3% and 14.7%, while the mean SDR₄ is 40.2%. Therefore, the SDR models can integrate these RUSLE models to obtain suitable results for the RSed, especially in inaccessible areas or areas with no verification. However, RUSLE 3 (includes C-factor Equation (18)) works better if it is coupled with SDR₁ (Scenario #18), while RUSLE 15 (includes land-cover-based model) works better with SDR₄ (Scenario #66; Table A3 in Appendix C and Table 12). In other words, the best combinations are represented by Scenario #18 (the best) and then Scenario #66 (better).

5.2. R Factor, C Factors, and RUSLE Uncertainties

Despite the variations and high uncertainty between the three models of the C factor, it appears that two models are valid in DLB: Equation (18) and land-cover models. Therefore, this study confirms that the selection of the C factor must consider the SDR model because of the high variation in the C-factor models (Figure 14B). This appears clearly in Figure 10, where the land-cover-based models (RUSLE 13-RUSLE18) have fewer distributions (low C-factor values) than the other RUSLE models. Moreover, increasing the number of NDVI scenes from one to three scenes per month will decrease the uncertainty of the C factor.

Table 12 emphasizes that the R factor does not play a big role in the discrepancy in the obtained results, especially for the six best scenarios of the RSed, where the same equation (i.e., Equation (5)) was used for Scenarios #18 and #66; Scenario #66 comes from land-cover-based models, and Scenario #18 comes from Equation (18). This is confirmed by Figure 14A, which displays that the variation between the six models is small, contrary to what appears in Figure 14B. However, the SDR models can be integrated with these RUSLE models to obtain suitable results for the RSed, especially in inaccessible areas or areas with no verification. Equation (18) works better if it is coupled with SDR₁ (Scenario #18), while the land-cover model works better with SDR₄ (Scenario #66; Table A3 in Appendix C and Table 12). In other words, the best combinations are represented by Scenario #18 (the best) and then Scenario #66 (better).

One of the variabilities of the results of the SRed is the uncertainty in the RUSLE models (Figure 15), which represents the distributions of the ~1% and ~5% randomly selected pixels from RUSLE model values for the eighteen scenarios. In both tests, i.e., Figure 15A,B, the pixels have low variation for the SL with high values, while the areas with low values have high variation. In other words, the uncertainty decreased with the SL values. The uncertainty of the RUSLE factors and SDR models is migrated to the RSed scenarios. The accuracy of the RUSLE factors and SDR models is significant; however, this requires field data for verification to evaluate these items [118], which requires carefully dealing with the local precipitation patterns and vegetation density [2], which is not an easy task [119].

5.3. Implications of This Study

Although 13 dams have been built within the DLB since 1978, the risk of RSed is present, as there was a decrease in the RSed amount from 0.46458 km³·yr⁻¹ to 0.3836 km³·yr⁻¹, i.e., 17.43% from the total amount. In comparison, the present catchment area for the DLB is decreasing to 36.2% of the original catchment (Table 1). It means that almost all of the areas (rugged Zagros mountains), which suffer from high and very high SL, are within the remaining catchments of the DL. Therefore, this study highly recommends keeping monitoring the RSed amounts in DLB and doing another bathymetry survey for this purpose. As a result, the new survey will emphasize our obtained findings, whose methodology can be applied and modified for other related and similar environments worldwide.

6. Conclusions

This study used the Revised Universal Soil Loss Equation (RUSLE)–Sediment Delivery Ratio (SDR) approach supported by multiple remote-sensing data sources to estimate the annual soil loss (SL) and reservoir sedimentation (RSed) in the Darbandikhan Lake Basin (DLB). For this purpose, eighteen RUSLE models coupled with five SDR models were used to estimate ninety scenarios of the RSed. The Digital Elevation Model (DEM), Tropical Rainfall Measuring Mission (TRMM), Harmonized World Soil Database (HWSD), and AQUA eMODIS NDVI V6 data were implemented for this end. The echo-sounding bathymetry survey was used to validate the results.

The research considered error factors, such as gravel and sand quarries, which may have a probable effect on decreasing the actual RSed. To check the scenario’s uncertainty, ~5% of the total pixels was used to investigate the relationship between the potion and the dispersion of the R and C factors estimation as a density plot.

The most proper model for the R and C factors and SDR was proposed in this research. Consequently, the most proper RUSLE-SDR model suitable for RSed modeling in the Zagros area was chosen based on the results from this research. Typical RUSLE (average of 8.3 t·ha⁻¹·y⁻¹) and RSed (average of 0.46458 km³·yr⁻¹) for the Zagros region were proved by the results of this research. In addition, this research recommended avoiding the use of the Index of Connectivity (IC) model for calculating SDR and land cover for calculating the C Factor.

This study recommended using the Modified Fournier index (MIF) equations for estimating the R factor. It can be concluded that the scenarios of RSed, which have high accuracy, can be classified into two groups: scenarios with high SDR and low RUSLE and scenarios with low SDR and high RUSLE. For future work, we recommended increasing the number of scenes (three scenes per month) when calculating the C factor to reduce the uncertainty.