New Insight on Soil Loss Estimation in the Northwestern Region of the Zagros Fold and Thrust Belt

: Soil loss is one of the most important causes of land degradation. It is an inevitable environmental and socio-economic problem that exists in many physiographic regions of the world, which, besides other impacts, has a direct bearing on agricultural productivity. A reliable estimate of soil loss is critical for designing and implementing any mitigation measures. We applied the widely used Revised Universal Soil Loss Equation (RUSLE) in the Khabur River Basin (KhRB) within the NW part of the Zagros Fold and Thrust Belt (ZFTB). The areas such as the NW Zagros range, characterized by rugged topography, steep slope, high rainfall, and sparse vegetation, are most susceptible to soil erosion. We used the Digital Elevation Model (DEM) of the Shuttle Radar Topography Mission (SRTM), Tropical Rainfall Measuring Mission (TRMM), Harmonized World Soil Database (HWSD), and Landsat imagery to estimate annual soil loss using the RUSLE model. In addition, we estimated sediment yield (SY) at sub-basin scale, in the KhRB where a number of dams are planned, and where basic studies on soil erosion are lacking. Estimation of SY will be useful in mitigation of excessive sedimentation affecting dam performance and watershed management in this region. We determined the average annual soil loss and the SY in the KhRB to be 11.16 t.ha − 1 .y − 1 and 57.79 t.ha − 1 .y − 1 , respectively. The rainfall and runoff erosivity (R factor), slope length (L factor), and slope steepness (S factor), are the three main factors controlling soil loss in the region. This is the ﬁrst study to determine soil loss at the sub-basin scale along with identifying suitable locations for check dams to trap the sediment before it enters downstream reservoirs. The study provides valuable input data for design of the dams to prevent excessive siltation. This study also aims at offering a new approach in relating potential soil erosion to the actual erosion and hypsometric integrals.


Introduction
Soil loss by water erosion is one of the most influential factors in land degradation, which cause significant environmental and socio-economic consequences that are critical to human welfare [1,2]. The evolution of the earth's landscape results from interplay among three main factors, which are lithology, climate, and tectonics [3,4]. These factors control, directly or indirectly, the erosional rate and soil losses [5]. Geomorphological processes, including erosion, continually modify the landscape, reduce area relief, and generate a vast amount of sediment that directly impact agricultural, ecological, and water resources.
Based on their lithological characteristics, different rock units, exposed in the Zagros region, offer varying levels of resistance to erosional forces, thereby exhibiting different erosional patterns [3]. In terms of climate, the northwestern part of the Zagros region is characterized by a heavy rainfall season preceded by a long dry season [6,7]. Collectively, the lithological and meteorological features sculpt the land, creating distinctive features, such as rugged topography with steep slopes [8], that occur in the northwestern Zagros [9]. Tectonics has also been a dominant factor in modifying the landscape of the Zagros Fold and Thrust Belt (ZFTB). For instance, high topography and steep slopes to the northeast of the Dezful and Kirkuk embayments result from intense thrusting caused by low strain at the two embayments [9,10]. This intense thrusting and the resulting steep slopes cause rivers to incise deep topography, yielding large volume of sediments, which can potentially affect the safety of the downstream dams and shorten reservoir life.
Hypsometry-the relation between forms of drainage basins (catchments) and their slopes [11]-is a function of river action, which, in turn, is controlled by lithology, climate, and tectonics. An individual basin can be represented by different intervals of elevation that enclose different areas. These intervals of elevations and areas form the hypsometric curve, which can be integrated into the Hypsometric Integral (HI) [12]. The amount of residual terrain above the lowest horizontal plane of a basin is the HI of a basin and can be used as a proxy for the erosional stage and landforms development [12,13]. HI has been used in several studies to express the relation between uplift and erosion and it gives reliable results on the degree of dissection to describe the erosional cycle through the youthful, mature, and old stages [14]. This allows the determination of active and inactive erosional regions.
Soil erosion removes the most fertile topsoil containing rich organic matter and nutrients, adversely impacting soil productivity [15,16]. An estimated 84% of land degradation worldwide results from erosion [17,18]. Due to the global nature of the problem and its direct impact on agricultural productivity, availability of a reliable and cost-effective method for estimating soil loss is of critical importance. Many models have been proposed to quantify soil loss that can be grouped into: (1) physical models [19,20], which are based on principles of water erosion; (2) conceptual models, that define significant human and natural factors, which affect the processes of erosion and sedimentation [21]; and (3) empirical models, which identify the different interacting erosion variables. Of these, the empirical models are easy to apply and require less computational time [22,23].
The popular Universal Soil Loss Equation (USLE) model was developed to estimate the amount of soil or sediment (mass of soil/sediment loss per unit area over time) reaching the outlet of a basin undergoing water erosion [24]. The USLE was purposed to foresee the mean rates of long-term erosion but not for prediction of soil loss events [22] where there is high fluctuation of soil losses from individual events [25] and such fluctuation cannot not be accurately included in the USLE model. This limitation has been removed in the Revised Universal Soil Loss Equation (RUSLE) by taking into consideration the rainfall energy [15]. Recent research, integrating RUSLE in the Geographic Information System (GIS) and remote sensing techniques, aided by the ever-increasing power and efficiency of computers [26], enables estimating soil loss for large areas quickly and at a reasonable cost [27][28][29][30][31][32][33][34][35][36][37][38] in different climatic zones [39].
The RUSLE model is widely used for soil loss estimation in many countries of the world due to its flexibility, direct application at field-scale, and value in establishing conservation measures and policies [40]. It is designed to estimate the average annual soil loss (A) based on five factors [41]; these factors are the rainfall erosivity factor (R), the slope length and steepness factor (LS), the conservation practice factor (P), the soil erodibility factor (K), the cover management factor (C). [42] stated that soil erosion can be classified into two types: potential erosion (PE) and actual erosion (AE). The R, K, and LS factors were examined as PE, whereas the C and P factors, in addition to R, K, LS factors, were considered as actual soil erosion [42]. We used the RUSLE model to estimate AE.
Soil erosion is the primary process for sediment delivery to rivers and streams. The amount of sediments exported at the basin outlet over a given period of time constitutes the sediment yield (SY) [43]. A barrier in the flow path, such as a dam, causes sediment deposition in the dam's reservoir in the form of bedload and suspended load [44]. Gradual accumulation of sediments leads to reduced reservoir storage capacity that lowers its useful life [45][46][47]. Mosul Dam, downstream of the present study area, is a good example of excessive reservoir siltation problem [7,48].
Several investigators have used the RUSLE model to estimate soil loss in neighboring countries (i.e., Iran [49][50][51][52][53][54][55] and Turkey [56,57]), and have applied RUSLE to determine the Sediment Delivery Ratio (SDR) [55,58]. However, RUSLE and SY models have not been used in the Kurdistan Region, north of Iraq, for soil loss estimation. Thus, our study serves to fill this critical knowledge gap and provides a reliable estimate of soil loss and SY in this part of the Zagros Range. In addition, our study also investigated the relationship between PE and HI at a basin scale.
This study also presents a new approach to relate PE-represented by the R, K, and LS factors-with rate of AE. It offers a useful tool to predict the effectiveness of soil loss control measures in a virgin area with no, or minimum human intervention. The resulting data should be valuable in implementing erosion control measures, such as ground cover management (C) and terracing and/or contour farming (P).

Study Area
The study area is located at the distal northwestern part of the ZFTB between 36 • 55 28 -37 • 47 01 N latitudes and 42 • 34 44 -43 • 29 5 E longitudes, along the border between Iraq and Turkey. It encompasses an area of~3500 km 2 and represents the Khabur River Basin (KhRB). The area is characterized by rugged topography and high relief to the north, due to the presence of the High Zagros Fault (HZF), and subdued topography and low relief in the southern part, near Dohuk city ( Figure 1). The area is characterized by its remoteness and difficulty of access, frequent landslides and lacks useful in situ hydrological and pedological information. The area encompasses the Khabur River that is the source of life-sustaining water to the region. The Zakho meteorological station data indicate significant seasonal fluctuation in the region, characterized by long, dry summers and short, rainy winter conditions. Most of the precipitation (586 mm) occurs between October and May, with a maximum monthly average of 134.3 mm occurring in January. The monthly average temperature and evaporation occurs in July, 33.96 • C, and 354.7 mm, respectively, and the minimum temperature of 4.17 • C was recorded in January ( Figure 2).
A number of dams have been planned in the KhRB and, like the existing Mosul dam with a high rate of reservoir siltation at 78,000 t.ha −1 .y −1 [59]. Excessive soil erosion in the region would create a major problem for the proposed dams, as pointed out by [48]. The rate of soil loss determined by this study will, therefore, enable the decision makers to incorporate suitable soil conservation measures to mitigate reservoir siltation problem, and also to optimize agricultural productivity.

Data and Software
We used four scenes of the mosaicked Digital Elevation Model (DEM), with 30 m spatial resolution (1 arc-second), gathered by Shuttle Radar Topography Mission (SRTM) [60]. DEM was used to extract the drainage network and calculate the slope and runoff flow accumulations. Due to the lack of in situ meteorological data, we employed Tropical Rainfall Measuring Mission (TRMM) data [61] to produce the rainfall map. The TRMM (3B43-V7) is a raster structure with 0.25° × 0.25° spatial resolution that collects monthly precipitation data [62]. The Harmonized World Soil Database (HWSD) [63] was used to determine soil erodibility, or the K factor. It has a cell size of 30 arc-second (~1 km). We

Data and Software
We used four scenes of the mosaicked Digital Elevation Model (DEM), with 30 m spatial resolution (1 arc-second), gathered by Shuttle Radar Topography Mission (SRTM) [60]. DEM was used to extract the drainage network and calculate the slope and runoff flow accumulations. Due to the lack of in situ meteorological data, we employed Tropical Rainfall Measuring Mission (TRMM) data [61] to produce the rainfall map. The TRMM (3B43-V7) is a raster structure with 0.25° × 0.25° spatial resolution that collects monthly precipitation data [62]. The Harmonized World Soil Database (HWSD) [63] was used to determine soil erodibility, or the K factor. It has a cell size of 30 arc-second (~1 km). We

Data and Software
We used four scenes of the mosaicked Digital Elevation Model (DEM), with 30 m spatial resolution (1 arc-second), gathered by Shuttle Radar Topography Mission (SRTM) [60]. DEM was used to extract the drainage network and calculate the slope and runoff flow accumulations. Due to the lack of in situ meteorological data, we employed Tropical Rainfall Measuring Mission (TRMM) data [61] to produce the rainfall map. The TRMM (3B43-V7) is a raster structure with 0.25 • × 0.25 • spatial resolution that collects monthly precipitation data [62]. The Harmonized World Soil Database (HWSD) [63] was used to determine soil erodibility, or the K factor. It has a cell size of 30 arc-second (~1 km). We combined seasonal impacts on vegetation using four Landsat Operational Land Imager (OLI) scenes. These scenes in Level 1T format, acquired in 2019, were used in this study. Details are listed in Table 1. The spectral images are of good quality, geo-referenced, and orthorectified with a 30 m spatial resolution, and are available without any charge, from the Earth Resources Observation and Science (EROS) center in Sioux Falls, South Dakota, U.S.A. Only scenes with excellent quality and less than 12.5 % cloud cover were chosen for this study. The region of interest is positioned in zone 38 north and the path and row of the image are 170 and 34 respectively. All input data were reprojected to the Universal Transverse Mercator (UTM) projection/WGS 1984 datum, within zone 38N using the nearest neighbor resampling method to fit with the OLI scenes. Combining the powerful capabilities of remote sensing and GIS, and using multiple datasets from different sources [64], we applied the RUSLE model to estimate soil loss in the study area. We used the MATLAB based TecDEM 2.2 [65,66] software to extract the drainage network. Specific functions within the ArcGIS 10.3 [67], such as slope, runoff flow accumulation, interpolation, raster conversion, and raster calculator, were utilized for base map preparation. In addition, we used Environment for Visualizing Images (ENVI) software for Landsat imageries preparation and processing that included implementing Boolean operation, surface reflectance conversion, and Normalized Differential Vegetation Index (NDVI) determination. R-studio script was used for statistical analysis.

Hypsometric Integral (HI)
The HI is a useful metric to identify different stages of landscape development in its geomorphological erosion cycle [12,68,69]. We used HI to determine the stage of landscape development (erosional stages) to quantify the relationship between HI and amount of soil loss. The mature (highly eroded) areas have low HI values, while the youthful (slightly eroded) areas show high HI values [12,14,[69][70][71][72]. According to [73], the HI of a given area can be estimated using Equation (1): where H max , H min , and H mean are the maximum, minimum, and mean elevations, respectively. An HI value close to 1 means that the terrain removal is in the lowest stage (uplift > erosion) and the land surface is in its youthful stage, while low HI values (close to 0) refers to high amount of terrain removal (erosion > uplift) and the land surface is in its mature stage. This dimensionless parameter is very useful for comparing the erosional stage of different basins regardless of their size (areas).

RUSLE Model Description
The RUSLE model, proposed by [41], can been used to estimate annual soil loss as shown in Equation (2).
The following sections explain the method of calculating RUSLE factors.

Rainfall and Runoff Erosivity (R factor)
The R factor expresses the effect of precipitation on soil erosion [74,75], as soil particle detachment is directly related to precipitation [76]. Estimation of R factor requires the availability of accurate and continuous precipitation data [22,31].
Since the study area lacks reliable meteorological data, we used the monthly TRMM (3B43-V7) to calculate the yearly average precipitation (mm), due to its suitable temporal and spatial resolution [30], and because of its extensive application for assessing the R factor, as compared to other satellite data. The TRMM data covers 15 years between September 2002 and August 2017. The TRMM pixels have been changed in vectors (points) format, and these points were interpolated employing an Inverse Distance Weighting (IDW) method, to produce a raster format with a spatial resolution of 30 m.
Numerous approaches have been used to calculate the R factor by employing longterm precipitation data. [77] found that Equation (3) can be used to calculate the R factor, which is directly related to the average annual rainfall.
where R is the rainfall erosivity factor in MJ.mm.ha −1 .h −1 .y −1 and P A is the average annual precipitation in mm. Precipitation in the KhRB amounts to 562.85 mm·yr −1 in the south, increasing to 785 mm·yr −1 to the north. The suitability of applying TRMM data in the study area was evaluated by comparing it with the observed precipitation dataset of Zakho meteorological station. Sixty-two monthly precipitation measurements, covering the period from September 2002 to December 2007 were compared with the TRMM data ( Figure 3). In this figure, a strong direct relationship was found, the coefficient of determination (R 2 ) being 0.853, at a significant p-value of <0.05. The slope and intercept were 0.9124 and 9.5485, respectively.

Soil Erodibility (K) Factor
The K-factor representing soil vulnerability to erosion, is affected by rainfall and runoff, [30,78]. Soil susceptibility to erosion is affected by its physical and chemical characteristics, such as its texture, mineralogy, and structure. Five soil parameters (i.e., texture, organic matter content (OM), coarse fragments, structure, and permeability) are the main properties that influence the K-factor. With the soil OM being inversely proportional to soil erodibility, such that any increase in soil OM decreases its vulnerability to detachment [36]. For this study, dataset of the HWSD was used to obtain information about soil texture and soil organic carbon (OC). We used Tables 2 and 3 to estimate the soil structure class and the soil permeability class, respectively [22,79]. We estimated the K-factor of the soil types in the study area using Equations (4) and (5) [22,79]. s is the soil structure class (s = 1: very fine granular, s = 2: fine granular, s = 3, medium or coarse granular, s = 4: blocky, platy or massive; Table 2); p is the permeability class (p = 1: very rapid, . . . , p = 6: very slow; Table 3).

Soil Erodibility (K) Factor
The K-factor representing soil vulnerability to erosion, is affected by rainfall and runoff, [30,78]. Soil susceptibility to erosion is affected by its physical and chemical characteristics, such as its texture, mineralogy, and structure. Five soil parameters (i.e., texture, organic matter content (OM), coarse fragments, structure, and permeability) are the main properties that influence the K-factor. With the soil OM being inversely proportional to soil erodibility, such that any increase in soil OM decreases its vulnerability to detachment [36]. For this study, dataset of the HWSD was used to obtain information about soil texture and soil organic carbon (OC). We used Table 2 and 3 to estimate the soil structure class and the soil permeability class, respectively [22,79]. We estimated the K-factor of the soil types in the study area using Equation (4)  s is the soil structure class (s= 1: very fine granular, s = 2: fine granular, s = 3, medium or coarse granular, s = 4: blocky, platy or massive; Table 2); p is the permeability class (p = 1: very rapid, …, p = 6: very slow; Table 3).    Estimating soil OM from the organic carbon (OC) is a simple procedure but requires the use of a conversion factor [80]. However, there is no agreement on the conversion factor [81] to use for converting soil OC to OM [82]. According to [80,82], a more reliable value of OM can be obtained by using the conversion factor of 1.724, which is based on the assumption that soil OM contains 58% carbon, as shown in Equation (6).
The L factor can be computed by using a slope (in percent) and runoff flow accumulation. We divided the slope map into four classes, which are <1%, 1-3%, 3-5%, and ≥5%.
where, λ is the horizontal projection of slope length (m), m is a constant, based on the value of slope; FA is runoff flow accumulation, Ps is the cell size of DEM, and θ is the slope angle in percent.
where S is the slope steepness factor, θ is the slope angle in percent, and LS represents the slope length and slope steepness factors [87].

Cover and Management Factor (C factor)
As a dimensionless and influential factor used in the RUSLE model, the C factor refers to the amount of soil erosion related to the degree of vegetation/ground cover. Evaluation of vegetation cover is considered the most significant aspect of estimating the C factor [88].
We first performed the radiometric corrections to compute top of atmosphere reflectance (ρ) by using (Equation (12)).
where ρ is reflectance (at sensor), L λ represents at-sensor radiance, D is sun-earth distance in astronomical units, θ is the solar zenith angle, and E sun,λ is the solar spectral irradiance. Next step involved atmospheric corrections that were performed using the Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) model. The FLAASH [89] model is based on a MODTRAN4 model to decrease the molecular and particulate scattering, absorption, and adjacency effects to retrieve at-surface reflectance values (ρ) [90][91][92]. The surface reflectance (ρ) of the OLI imageries is used to compute NDVI following Equation (13) [93], which, in turn, was used to compute the C factor using Equation (14) [ [94][95][96].
where α, β are constants, which were deduced from [94,95] to be 2 and 1, respectively. These constants determine the shape of the NDVI-C curve. We omitted the bare rock in the soil loss assessment, affirmed by relatively low values of NDVI [32].

Support Practice Factor Related to Slope Direction (P)
The P factor describes soil loss ratio by a specific soil conservation or support practice to the corresponding loss up-and down-slope [23,97]. Lack of in situ data for conservation practices within the Zagros region prompted us to apply Equation (15) for estimating the P factor. We modified Equation (15) as suggested by Wener [33,37,98] and used the modified equation (Equation (16)) to calculate P.

Sediment Yield (SY)
We estimated the SY in the KhRB following [99]. First, we measured the areas of the sub-basins within the KhRB. Then, we calculated the SDR (Equation (17)). Finally, we estimated the SY amount for each sub-basin following Equation (18) [99]. Results of SY were used to estimate the amount of sediment transported by streams and rivers to the downstream area to determine the check dam locations.
SDR is sediment delivery ratio; A b is area of a basin (km 2 ); AE is the actual erosion (yearly rate of soil loss (t.ha −1 .y −1 )); SY is the sediment yield (t.ha −1 .y −1 ).

Work Procedure
We applied the RUSLE model to estimate soil loss for the KhRB. The KhRB was subdivided into a number of sub-basins ( Figure 4A). We estimated the mean of the AE, the PE, and the mean of the HI, for each sub-basin ( Figure 4B). Afterwards, we tested the relationship between AE and PE for them.
We classified various parameters in the RUSLE model into five classes, using the natural breaks technique, except the K factor that was classified into three classes based on the three soil types present in the KhRB. These classes are: very high, high, moderate, low, and very low. The natural break classification is best applied to familiar data ranges. It reduces the variance within classes and maximizes the variance between classes [72,100]. This characteristic makes the natural breaks useful when RUSLE and its factors histogram show distinct breaks. We classified various parameters in the RUSLE model into five classes, using the natural breaks technique, except the K factor that was classified into three classes based on the three soil types present in the KhRB. These classes are: very high, high, moderate, low, and very low. The natural break classification is best applied to familiar data ranges. It reduces the variance within classes and maximizes the variance between classes [72,100]. This characteristic makes the natural breaks useful when RUSLE and its factors histogram show distinct breaks.   Figure 4B displays the HI for the study area. The HI values range from 0.36 to 0.05; low values (<0.30) were obtained for the central and lower parts of the KhRB, indicating plain terrain with occasional slope. The HI values above 0.30 fall in the part of basin located inside Turkey. The topography of these areas possesses high relief, which might be controlled by the active thrust fault, HZF (Figure 1).

The R Factor
The value of R factor in the study area varied from~296 to~393 MJ.mm.ha −1 .h −1 .y −1 , with an average of~354 MJ.mm.ha −1 .h −1 .y −1 . The R factor map of KhRB, presented in ( Figure 5A), has been sliced into five groups, to show relatively lower R factor values in the south of the KhRB compared to the mountainous terrain in the north.

The LS Factor
The highest slope angle in the KhRB is 77°, the average being 17.5°. Nearly all high slope areas are located in the central and northern parts of the KhRB. Due to the difference

K Factor
The K factor data presented in Figure 5B show that the minimum and maximum soil erodibility is~0.023 t·ha·h·ha −1 ·MJ −1 .mm −1 in the southwestern part of the KhRB and 0.063 t·ha·h·ha −1 ·MJ −1 .mm −1 in the southeastern parts of the KhRB respectively. Major part of the basin has an intermediate to high K value (~0.0498) occurring in the central and northern part of the basin. Three types of soil occur in the KhRB: lithosols, calcic xerosols, and chromic vertisols. The textural description and the K-factor for each soil type is shown in Table 4 and Figure 5B. Variations in the K-factor values are possibly due to variations in the soil types' physical characteristics.

The LS Factor
The highest slope angle in the KhRB is 77 • , the average being 17.5 • . Nearly all high slope areas are located in the central and northern parts of the KhRB. Due to the difference in runoff flow accumulation, and slope distribution, the LS factor values range from 0 to 1619.6 with the average and standard deviation 12.05 and 12.74, respectively ( Figure 6A). The KhRB, generally, is characterized by moderate LS factor values: It is <10 in the lower part of the KhRB, but the middle and upper parts of the basin show relatively higher LS factor values, which indicate the possibility of greater soil erosion and consequently accelerated land degradation.

The C Factor
The C factor ranges from 0 to 0.5999 ( Figure 6B). Less than 16% of the KhRB is covered with vegetation, and most of the barren land occurs in the upstream and downstream parts of the basin ( Figure 7A) where C factor values are low (dark green and green colors). Presence of vegetation increases soil particles' adhesion and binding, making it less prone to erosion. The moderate, high, and very high (>0.2) values of the C factor (green, orange, and red colors) are mainly restricted to the midstream part of the KhRB. Four types of vegetative areas exist in the KhRB ( Figure 7B). Open shrublands (dominated by woody perennials (1-2 m height), 10-60% cover), savannas (tree cover 10-30% (canopy > 2 m)), grasslands (dominated by herbaceous annuals (<2 m)), and croplands (at least 60% of area is cultivated cropland) [101].

The P Factor
Values of the P factor within the KhRB range between 0 and 1 ( Figure 8A). Low values correspond to pixels with good conservation practice, while the high values denote

The P Factor
Values of the P factor within the KhRB range between 0 and 1 ( Figure 8A). Low values correspond to pixels with good conservation practice, while the high values denote areas of poor conservation practice [31]. P factor value depends on slope length and steepness, and farming method [102]. Therefore, an increase in the value of P factor is consistent with the rise in the LS factor value (Figures 6A and 8A). High and very high (0.35-1) values of P factor are mainly distributed in the middle part of the KhRB, while the low and very low values (0-0.22) are mainly distributed downstream of the KhRB ( Figure 8A).

Estimation of Soil Loss
We integrated the RUSLE modelʹs parameters within the GIS environment to evaluate the average annual soil loss in the KhRB and estimate it to be 11.16 t.ha −1 .y −1 (Figure

Estimation of Soil Loss
We integrated the RUSLE model's parameters within the GIS environment to evaluate the average annual soil loss in the KhRB and estimate it to be 11.16 t.ha −1 .y −1 ( Figure 8B). The minimum and maximum yearly soil losses are 0 t.ha −1 .y −1 and 3617.6 t.ha −1 .y −1 , respectively. Total annual soil loss for the whole study area is estimated at 31,237,283.358 t.y −1 . The low and very low range (<10 t.ha −1 .y −1 ) of the yearly soil loss values mainly occur in the northeastern and downstream reaches of the KhRB, while the moderate, high and very high range (>10 t.ha −1 .y −1 ) occurs in the midstream and upstream parts of the KhRB ( Figure 8B). Very high soil losses (>50 t.ha −1 .y −1 ) occur only in 6.19% of the KhRB.

Relationship of PE to RUSLE and the HI
Our study revealed a strong direct correlation between: (a) PE and soil loss, and (b) PE and HI. It is well recognized that heavy rain, and steeper slopes contribute to a large volume of runoff with high erosive energy and greater amount of soil loss. Sixtyone sub-basins included in the KhRB, comprise areas ranging between 3.5 and 974.6 km 2 (average: 57.5 km 2 ). Analyses of annual soil loss within the KhRB at the sub-basin scale reveal a moderate linear relationship between the mean of the PE and mean of soil loss. This direct relationship is clearly displayed in Figure 9, where the statistical correlation, the p-value, the t-test, the slope, and the intercept of the data were >0.511, 9.572 × 10 −11 , 7.8565, 0.02979, and 3.56121, respectively. Using the same procedure, we found a moderate linear relationship between the mean of HI and mean of soil loss (Figure 10), where the statistical correlation, the p-value, the t-test, the slope, and the intercept of the data were 0.522, 3.315 × 10 −10 , 7.6743, 0.0003911, and 0.0921817, respectively.

Relationship of PE to RUSLE and the HI
Our study revealed a strong direct correlation between: (a) PE and soil loss, and (b) PE and HI. It is well recognized that heavy rain, and steeper slopes contribute to a large volume of runoff with high erosive energy and greater amount of soil loss. Sixty-one subbasins included in the KhRB, comprise areas ranging between 3.5 and 974.6 km 2 (average: 57.5 km 2 ). Analyses of annual soil loss within the KhRB at the sub-basin scale reveal a moderate linear relationship between the mean of the PE and mean of soil loss. This direct relationship is clearly displayed in Figure 9, where the statistical correlation, the p-value, the t-test, the slope, and the intercept of the data were >0.511, 9.572 × 10 −11 , 7.8565, 0.02979, and 3.56121, respectively. Using the same procedure, we found a moderate linear relationship between the mean of HI and mean of soil loss (Figure 10), where the statistical correlation, the p-value, the t-test, the slope, and the intercept of the data were 0.522, 3.315 × 10 −10 , 7.6743, 0.0003911, and 0.0921817, respectively.

Estimation the Sediment Yield
In addition to SDR ( Figure 11A), we also calculated the average annual quantity of SY for the 61 sub-basins included in the KhRB. The SDR values range between 0.619 and 1.373, respectively. Estimate of the yearly average of SY in the KhRB is 57.79 t.y −1 . The minimum and maximum yearly SY are 6.9 t.y −1 and 165.12 t.y −1 , respectively. Total annual SY for the whole study area is estimated at 21,817.68 t.y −1 . The low and very low range (<10 t.ha −1 .y −1 ) of annual SY values mainly occur in the northeastern and downstream reaches of the KhRB, while the moderate, high and very high range (>10 t.ha −1 .y −1 ) occur in the midstream part of the KhRB ( Figure 11B).

Estimation the Sediment Yield
In addition to SDR ( Figure 11A), we also calculated the average annual quantity of SY for the 61 sub-basins included in the KhRB. The SDR values range between 0.619 and 1.373, respectively. Estimate of the yearly average of SY in the KhRB is 57.79 t.y −1 . The minimum and maximum yearly SY are 6.9 t.y −1 and 165.12 t.y −1 , respectively. Total annual SY for the whole study area is estimated at 21,817.68 t.y −1 . The low and very low range (<10 t.ha −1 .y −1 ) of annual SY values mainly occur in the northeastern and downstream reaches of the KhRB, while the moderate, high and very high range (>10 t.ha −1 .y −1 ) occur in the midstream part of the KhRB ( Figure 11B).

Discussion
It is well established that soil loss is controlled by the five factors of the RUSLE model (i.e., R, K, LS, C, and P). Change in one of these factors leads to a change in the quantity of soil loss. As can be seen in Equation (3), the R factor is directly related to the amount of precipitation [77]. Moreover, in the KhRB, precipitation is affected by elevation (Figure 1). The soil loss is a function of K factor, rainfall, and runoff [30,78]. The higher the K factor, the greater is the amount of soil loss and vice versa [103]. The same is true for the LS factor, where the amount of soil loss increases with increasing LS factor and vice versa [85,87]. Our results show moderate to high soil loss in areas with high R factor values ( Figure 5A), which is consistent with relatively steeper slope zones ( Figures 6A and 8B). Soil loss is a function of the products of all five factors in RUSLE, of which R, K, LS are natural factors

Discussion
It is well established that soil loss is controlled by the five factors of the RUSLE model (i.e., R, K, LS, C, and P). Change in one of these factors leads to a change in the quantity of soil loss. As can be seen in Equation (3), the R factor is directly related to the amount of precipitation [77]. Moreover, in the KhRB, precipitation is affected by elevation (Figure 1). The soil loss is a function of K factor, rainfall, and runoff [30,78]. The higher the K factor, the greater is the amount of soil loss and vice versa [103]. The same is true for the LS factor, where the amount of soil loss increases with increasing LS factor and vice versa [85,87]. Our results show moderate to high soil loss in areas with high R factor values ( Figure 5A), which is consistent with relatively steeper slope zones ( Figures 6A and 8B). Soil loss is a function of the products of all five factors in RUSLE, of which R, K, LS are natural factors and cannot be changed, P and C are factors that can be altered by human intervention. This is precisely what is done to minimize soil loss in agricultural land [104,105]. Therefore, by taking measures to reduce values of P and C factors by contour ploughing, no-till seeding, and terracing, etc. the overall soil loss can be reduced.
The C factor is relatively lower in the southern part of the KhRB ( Figure 6B) because of sparse vegetation, somewhat gentler slopes, and lower rainfall amounts-all of which contribute to lower rate of soil loss ( Figure 8B). On the other hand, the central and the northwestern parts of the KhRB, with lower values of the C factor, due to presence of vegetation, exhibit a relatively higher amount of soil loss due to relatively steeper slope and higher amounts of rainfall. Dynamics of the erosional process are greatly affected by the LS factor that controls runoff flow. Our results agree with the findings reported by [106] in which the altitude variance and topographic relief control the LS factor ( Figure 6A); and higher value of P factor ( Figure 8A) directly correlates with a higher rate of soil loss ( Figure 8B). Interestingly, both of which are related to the steeper slopes ( Figure 6A).
This study also concurs with [8], which stated that the main factor controlling soil loss "in the NW part of the ZFTB" is the LS factor ( Figure 6A). The next important factor that affects soil loss is the R factor ( Figure 5A). There is a narrow range in the values of the K factor (~0.023-0.063 t·ha·h·ha −1 ·MJ −1 .mm −1 ), C factor (0.2-0.87), P factor (0-1), and R factor in the sub-basins. The LS factor for~99% of the pixels in the sub-basins gave values of <166, except some that show an extremely high value of about 1500.
From the above discussion of the various RUSLE parameters, it can be stated that slope is the main proxy that controls soil loss in the region, where precipitation is strongly affected by the slope. Regional tectonism, evidenced by the presence of HZF in the study area (Figure 1), is believed to have contributed to development of steeper slopes.
Although the RUSLE model is widely used to quantify soil loss and may be more accurate than other methods [107], it has not been applied at pixel scale to estimate PE in the study area. Our study revealed direct relationship between potential soil loss and actual soil loss at the sub-basin scale, but a with moderate-strong value of 0.51 for the coefficient of determination, R 2 ( Figure 9). This may be due to the difference in precipitation, soil types, and vegetation density in the various sub-basins.
Overall, this study shows that soil loss in the study area increases with an increase in slope angle, rainfall, and lack of vegetation, and is more pronounced in regions of steeper slopes. This study agrees with [108], which reported that different slopes initiate different rates of soil loss, and progressively higher soil loss occurs on steeper slopes. Our study indicates that the direct relationship between PE and HI can be formed in quantitative relationships ( Figure 10). This relation confirms the finding of [14,71] that the sub-basins in youthful stage are more active than mature stage. This allows the determination of active and inactive erosional regions, which will be helpful to select the best sites for check dams.
Further studies should be carried out to ensure safety of existing dams in Iraq. Regarding potential siltation problems in the proposed dams in Iraq, research should be undertaken to determine annual volume of soil that will enter the reservoirs to ensure safety and longevity of the planned structures.

Estimation of the Dam Siltation in the KhRB
Estimation of soil loss at sub-basin scale provided reliable data to find suitable sites for locating main dams and the check dams. We applied RUSLE with GIS and RS at the subbasin scale for the first time in the KhRB to estimate soil loss. Since the planned dams are located in the sub-basins ( Figure 11B, Table 5), they would experience excessive reservoir siltation due to high rate of soil loss. To minimize this problem, we have recommended check dams to trap the sediments upstream of the main dams. Table 5. Features of the proposed dam (modified after [48]). The dam sites (4 and 9) are excluded because they are outside of the study area. One of the common and simple ways to mitigate excessive siltation is construction of small check dams. [48] had recommended suitable locations for dams in the KhRB ( Figure 11B and Table 5) based on lithology, tectonic zones, distance to active fault, distance to the lineaments, soil type, land cover, elevation, slope gradient, precipitation, stream width, curve number, distance to road, distance to towns and cities, and distance to villages [109][110][111]. These dam sites might be influenced by soil loss from the upstream areas leading to excessive reservoir siltation. Drawing upon our present research, we propose construction of five check dams to trap sediments transported from the upstream area of the KhRB (Figure 11). These dams will prevent~1757 t.y −1 of sediment, eroded from the upper reaches of the KhRB, from entering the main dams. The suggested check dams will prolong the lives of proposed dams ( Figure 11B), and also increase the useful life of the Mosul dam to the southwest of the KhRB.

Conclusions
The classical RUSLE model has been integrated with remote sensing and GIS to estimate soil loss and sediment yield in the KhRB. Results show the spatial distribution and rate of average annual soil loss and sediment yield within selected basins that ranges from 0 t.ha −1 .y −1 to 3617.6 t.ha −1 .y −1 , and from 6.9 t.y −1 to 165.12 t.y −1 , respectively. The northwestern part of KhRB experiences low rate of soil loss (A) and SY. These rates are mainly influenced by topography, characterized by steep slopes, and high rainfall. Our study also demonstrates that LS and R are the most dominant factors controlling soil loss in the selected area. Analysis of the potential and actual soil erosion for KhRB shows a direct relationship. We conclude that RUSLE can serve as an effective model for estimating sediment yield, and for mitigation its adverse effects on multi-purpose river valley development projects for water supply, agriculture, etc. We recommend locating check dams in the upstream reaches of KhRB to minimize reservoir siltation and ensure safety of the main dams. Our study showed a correlation between potential soil erosion and HI. We found that HI values at the sub-basin scale can be used to identify active and inactive erosional regions. Thus, HI could serve as a useful geomorphic index for locating check dams in basins with high rate of soil loss. Further studies should be performed to find out if the HI can be used for estimating actual soil loss at specific locations.
Author Contributions: Arsalan Ahmed Othman outlined the research, prepared and processed the data, accomplished the study, and wrote the manuscript. Ahmed K. Obaid processed the data and supported the writing, analysis, and discussion. Diary Ali Mohammed Amin Al-Manmi, Ahmed F. Al-Maamar, Ahmed T. Shihab, and Younus I. Al-Saady supported the analysis and RUSLE model. Veraldo Liesenberg and Syed E. Hasan supported the writing and discussion. All authors have checked and revised the manuscript. All authors have read and agreed to the published version of the manuscript.