Torrential Flood Water Management: Rainwater Harvesting through Relation Based Dam Suitability Analysis and Quantification of Erosion Potential

In this study, a relation-based dam suitability analysis (RDSA) technique is developed to identify the most suitable sites for dams. The methodology focused on a group of the most important parameters/indicators (stream order, terrain roughness index, slope, multiresolution valley bottom flatness index, closed depression, valley depth, and downslope gradient difference) and their relation to the dam wall and reservoir suitability. Quantitative assessment results in an elevation-area-capacity (EAC) curve substantiating the capacity determination of selected sites. The methodology also incorporates the estimation of soil erosion (SE) using the Revised Universal Soil Loss Equation (RUSLE) model and sediment yield at the selected dam sites. The RDSA technique identifies two suitable dam sites (A and B) with a maximum collective capacity of approximately 1202 million m3. The RDSA technique was validated with the existing dam, Gomal-Zam, in the north of Sanghar catchment, where RDSA classified the Gomal-Zam Dam in a very high suitability class. The SE estimates show an average of 75 t-ha−1y−1 of soil loss occurs in the study area. The result shows approximately 298,073 and 318,000 tons of annual average sediment yield (SY) will feed the dam A and B respectively. The SE-based sediment yield substantiates the approximate life of Dam-A and Dam-B to be 87 and 90 years, respectively. The approach is dynamic and can be applied for any other location globally for dam site selection and SE estimation.


Introduction
Flooding is a global issue and is becoming a prime concern of research comprehensively. Flooding enriches grounds with minerals and sediments transported by water. These sediments increase the fertility of the soil and replace long-standing with new soil [1]. Floods are a great natural way to recharge groundwater [1,2]. Despite such benefits, when they cause human and economic loss, they are viewed as "catastrophes". For better flood management researchers around the globe are dealing with procedures to utilize the flood water and control loss of assets. These studies presented structural and non-structural management of flood. The structural management considers physical measures including the construction of dams, leaves, floodwalls, and elevated buildings, cleaning of water bodies, and flood proving properties [3,4]. The second possible solution is non-structural, which includes planning for disaster, flood plain zoning, and early warning systems. These measures cannot completely utilize flood resources, this is the explanation integrated water resource management (IWRM) is standing out enough to be noticed [4].
Flash floods are a potential source of freshwater. However, flash floods are seldom considered as source of water due to the unavailability of management resources. Flash floods are different from traditional riverine floods in their extent and properties [5][6][7][8][9][10][11]. These floods last for a short time but are intense in the extent of destruction they cause. Resources from these floods cannot be as proficiently saddled as those from riverine floods.
techniques of GIS and RS [17]. The implementation of RUSLE for the estimation of SE based on the integration of GIS and RS was also used in different studies [33].
Torrential flooding and erratic behavior of rainfall in western catchments of Dera Ghazi Khan have created a myth for disaster risk reduction strategies. The unplanned urbanization and arid nature of the area further limited the availability of freshwater resources. Suitable site identification for dam construction is the utmost need for torrential flood management and to meet the water demands of the area. The study aims to provide effective management of hill torrents' flow, converting the disastrous energy into a useful source by identification of suitable sites for storage structures. To select a group of indicators that exceedingly affect the dam suitability, different literature was investigated. The selected indicators follow hydrologic and engineering investigations in the spatial domain. The methodology adopted is based on the relation of selected indicators with reservoir and dam wall suitability. Focusing on the main problem of the area, the study also utilizes the GIS and RS technique to quantify the sediment yield and calculating the life of proposed structures.

Study Area
Approximately, 200 hill torrents are counted that originate from the Suleman mountain range affecting Dera Ghazi Khan (D G Khan, Pakistan), Dera Ismael Khan, and Rajanpur districts. Among these torrents, 13 are classified as belonging to a major category as shown in Figure 1. These torrents are usually classified as small, medium, and major. The classification is based on average annual peak flows (small ≤ 5000; 5000 < Medium < 15,000; Major > 15,000) cfs [34]. The major torrents affecting D. G. Khan include Kaura, Vehowa, Sanghar, Sori Lund, and Vidor. The characteristics of torrents including; variation in peak flows, and high sediment yield hardened the management of these torrents [3].  The steep slope mountainous catchment is primarily majorly composed of loam, clay loam, and sandy clay loam. Desert plain and barren land cover a major part of Sanghar catchment. However, low-density shrub covers also exist in some parts of the catchment.
The torrential stream debouching onto the "pachad/piedmont area" from the Sanghar catchment. The stream locally named "Rod Kohi" is categorized as non-perennial and remains active in monsoonal season. The total length of the torrential stream is approximately 158 km [3,4,34]. The discharge data (Table 1) collected from Punjab Irrigation Department (2010-2019) shows that event concentration is clustered in the monsoon season (June, July, and August). The area is arid and water scarcity has not only limited agriculture but also socioeconomic activities. During the monsoon season, flash flooding occurs and runoff generated by these torrential streams is far more than the local manageable capacity of farmers. The flows during monsoon season remain unutilized and frequently damage the low-lying urban settlements and patchy crops in piedmont plains [34]. The Sanghar is the most massive torrent in D G Khan District and the second-largest in the belt of D G Khan and Rajanpur districts with an approximate area of 4900 km 2 . The area is divided into three physiographic regions including (i) Mountain range (catchment area), (ii) Pachad/Piedmont plain (plain between foothill and canal command area), and (iii) Canal command area (CCA).
The steep slope mountainous catchment is primarily majorly composed of loam, clay loam, and sandy clay loam. Desert plain and barren land cover a major part of Sanghar catchment. However, low-density shrub covers also exist in some parts of the catchment.
The torrential stream debouching onto the "pachad/piedmont area" from the Sanghar catchment. The stream locally named "Rod Kohi" is categorized as non-perennial and remains active in monsoonal season. The total length of the torrential stream is approximately 158 km [3,4,34]. The discharge data (Table 1) collected from Punjab Irrigation Department (2010-2019) shows that event concentration is clustered in the monsoon season (June, July, and August).

Geospatial Dataset
The data used for site selection were collected from different online and as well as from government organizations ( Table 2). For dam site suitability analysis, digital elevation model (DEM) is the most important dataset. In this study, DEM product of 5 m spatial resolution was used to generate by-products indicators including multi-resolution valley bottom flatness (MRVBF), terrain roughness index (TRI), closed depression (CD), valley depth (VD), gradient to downslope (GD), and length slope factor (LS). The DEM product was collected from Punjab Irrigation Department.

Indicator's Relation and Suitability
The spatial suitability of storage structures is directly associated with morphometry and the terrain of the area. The indicators including stream ordering (SO), multiresolution valley bottom flatness index (MRVBF), terrain roughness/ruggedness index (TRI), slope, closed depression Index (CDI), and Valley depth (VD) thematic layers are used for reservoir suitability. However, TRI, slope, and downslope difference gradient (GD) were used for dam wall suitability analysis. The selection of indicators is solely based on the criterions defined by Munir et al. [4] and Stephens [31]. The criterions are a set of rules defined in the literature that should be followed for dam site suitability analysis. Each selected indicator represents certain criteria.
Stream ordering (1st criteria) was performed using the Strahler order scheme. It was of utmost importance to classify the streams in order to identify the length and position of the main channel in the study area. MRVBF index was developed by Gallant and Dowling [35]. The index classifies the flat bottom of valleys as low areas. Higher values of MRVBF represents the flat valleys along the stream. According to the defined relations, flat areas along the main stream were identified using MRVBF index (2nd criteria). TRI (3rd Criteria) quantifies heterogeneity-abruptness in the terrain over a specified length [36], is calculated as Slope (4th criteria) of the area also impacts the site suitability analysis. The selected site should have gentle/low main slope values for reservoir and steep cross slope are suitable for dam wall [4].
Generally, the topography of the area commands the choice for the type of dam. A perfect gorge (narrow valley with steep cross-sectional slope) favors a concrete dam. However, rolling plains with wider valleys favors an embankment dam. The VD (5th criteria) and CD (6th criteria) indicators are utilized for the reservoir suitability, to examine the reservoir's ability for water storage. CD indicator also represent the contour rule defined by [4,31]. However, the GD index quantifies downslope controls on local drainage [37]. The relationships of selected indicators with site suitability were defined based on a literature review [4,31,[38][39][40][41] engineering as well as hydrologic investigations, and field expert's counseling.
Buffer analysis is performed on the main ordered stream to develop a rectangular box unit. The box unit is 250 m in length along the mainstream and the width of the buffer is selected as 100 m; however, it may vary for different areas depending upon the topography of the mainstream. The average values of each indicator are extracted for n-box units along the mainstream and are classified into five classes based on their relationship as shown in Figure 2. The combination of yellow with 0% is the weakest relation with suitability and a red with a 100% effect gives a strong relation with site suitability. The classification of very low (VL) to very high (VH) is user-dependent and it can be replaced with numeric values ranging from 1 to 5. raphy of the mainstream. The average values of each indicator are extracted for n-box units along the mainstream and are classified into five classes based on their relationship as shown in Figure 2. The combination of yellow with 0% is the weakest relation with suitability and a red with a 100% effect gives a strong relation with site suitability. The classification of very low (VL) to very high (VH) is user-dependent and it can be replaced with numeric values ranging from 1 to 5. The average extracted values of each indicator further analyzed using the simple additive weighting (SAW) technique for final decision making of storage structures. The technique solves the weighted performance of each indicator for individual box unit. SAW requires a normalization process either performed as a decision matrix or on the actual values. In this study, SAW technique was modified and performed on actual values of each indicator. The weight factor in SAW technique is divided into Inter weights (weights assigned to different indicators) and intra-weights (weights assigned to value/box units of individual indicators) categories. The technique follows the relation as under: r ij = Normalized values of indicators X ij = Actual thematic layer values W inter = Inter weights between selected indicators W rr intra = Intra weights for each box unit The W rr intra are calculated for all box units of selected indicators individually using the reciprocal rank technique. In this technique, weights are derived from the normalized reciprocals of ranks. The technique involves an ordinal number to each record/value rank, starting with the highest rank of 1 [42].
where W rr intra is intra-weight using rank reciprocal, i = individual value of the indicator, j = ordinal number, and N = total number of box units. Furthermore, W inter weights are assigned on equal weight strategy because the selected indicators influence equally on the selection of suitable site.
The final SAW output (V) for each box unit is classified in ArcGIS 10.5 using standard deviation (SD) to minimize the manual dependency of users. The SD strategy makes the technique a widespread methodology where the dataset is free of specialists for ar-rangement and classification [43]. In the SD classification scheme, the number of classes is dependent upon the data spread and distribution. The symbology of classes is userdependent. In this study, the number of classes generated through SD strategy further symbolized as very low (VL) to very high (VH). The combination of reservoir and dam wall suitable sites results in final output for storage structures.

Quantitative Analysis
The capacity and storage characteristics are required for the effectiveness of selected structures in flood scenarios. Elevation-Area-Capacity (EAC) curves are generated using (i) elevation surface area method designed for reservoirs and (ii) Trapezoidal rule implication for capacity calculation at each contour area [4,44,45]. The mathematical expression of the trapezoidal rule is given in the following equation: where h 2 , h 1 , n = value of the largest, smallest and total number of contour.

Sediment Analysis and Dam Life
Sediment erosion (SE) analysis is performed using the RUSLE equation to calculate the life of proposed structures in the area. The technique is an empirical approach for predicting the long-term rate of erosion potential expressed as a function of five indicators including (i) rainfall erosivity (R), (ii) soil erodibility (K), (iii) slope length and steepness (L and S), (iv) support practices (P), and (v) cover management practices (C) inputs [46][47][48]. The Sanghar catchment area mostly consists of a barren mountain with little or no cropping system. Therefore, the average SE values are mainly or solely depending upon the LS, R, and K factors.
The multi-source gathered and created parameters of RUSLE equation with various spatial resolution goes through bilinear resampling measure. The cycle ascertains the estimations of every pixel by averaging the values of encompassing four pixels. The resampling cycle permits to downscale the precipitation raster to coordinate the spatial resolution with other datasets.

Rainfall Erosivity (R) Factor
The intensity and extent of every single rainfall event are expressed as R factor. Satellite rainfall data on daily time scales of GPM (global precipitation measurements) product is used for the calculation of R factor. The product is widely used due to its enhanced accuracy. The gauge data of nearby stations (D.G. Khan) is used to validate the satellite data. In this study, R-factor was calculated using the proposed scheme of [49].
where coefficient k and exponent a are model parameters, and P d is daily precipitation. The R is calculated with daily precipitation for all days P d > 0 in a month is summed. However, coefficient k varies both spatially and temporally [50]. Ref. [51] showed that exponent a is close to 2, while a number of empirical approaches substantiate the value of a between 1.5 and 2 [52].

Soil Erodibility (K) Factor
The soil susceptibility to erosion due to rainfall and generated runoff is justified by the K factor. The factor is calculated from different empirical equations [53,54] incorporating the soil properties like percentage organic matter, soil texture, and permeability.
where a, b, and c are percentage organic matter, soil structure code, and soil permeability code, respectively. However, the percentage content of organic matter (OM) is related to organic carbon (OC) content. OM content is calculated using the relation defined by [55,56]. The percentage of sand, silt, clay, and organic carbon is collected from soilgrid and Harmonized World Soil Database (HWSD). The factors b and c are calculated using soil texture classification [57]. The classification is based on the texture of sand, silt, and clay content as shown in Table 3. Factor m in the above equation represents the texture percentage of soil as:

Topographic LS Factor
The LS-factor explains how topography affects the SE process. The L and S factors describe the impact of slope length and slope steepness respectively. The required SE model is calculated through the combined L and S factors (LS-factor). The slopes greater than 9% are more prone to soil loss [58]. The methodology adopted the Renard et al., 1997 [59] algorithm for the estimation of the S-factor based on slope gradient.
where θ is a slope in degrees The L-factor is calculated using the model proposed by Desmet and Govers 1996 [60]. The model extended Foster and Wischmeier's 1974 [61] approach to a 2-dimensional terrain [62].
where A i,j−in the area is contributing to the inlet grid cell (i, j) measured in m 2 . G is cell size in meters. α i,j is the aspect of the grid cell (i, j). However, m is related to the ratio of the rill to interill erosion (β). and θ is a slope in degrees. The value of m ranges from 0 to 1 [62].

Cover (C) and Control Practice (P) Factor
The impact of landcover on soil erosion is estimated using C-factor. The satellite data of Landsat-8-Operational Land Imager and Thermal Infrared Sensor (OLI/TIRS) was used for pixel-based LULC (Land use land cover) classification. The available cloud free Landsat images were acquired for 3 March 2020. The Sanghar catchment was covered in two Landsat images of Paths 152, and 151 and Row 39. The images consist of nine spectral bands with a spatial resolution of 30 m. The band 2 to 7 (visible to short wave infrared) were stacked and mosaiced to generated a multiband composite satellite image of the study area. The training sample from mosaiced image were collected and clustered by analyzing the identical spectral behaviour of pixels. The sample pixels with identical spectral behavior were labelled to the particular LULC class. The multiband composite image with different clustered training samples was further processed in Erdas Imagine 2014 using ML (maximum likelihood) algorithm for LULC classification. The ML algorithm assumes that the statistics for every class in each band are normally distributed and calculates the probability that a given pixel belongs to a selected class. Each pixel is assigned to the LULC class that has the highest probability. The values of the C-factor against each LULC class were collected from [63,64] and were then joined with classified raster. The Table 4 shows the assigned C-factor values for each LULC class. Table 4. C-factor values.

Sr. No
Land Cover C-Factor The exercises/practices to minimize erosion are described by P-factor. The Sanghar catchment is a hilly catchment and has never been practices for erosion management; subsequently P-factor is taken as 1 for the whole catchment.

Sediment Yield (SY) and Sediment Delivery Ratio (SDR)
The SY is the gross SE delivered to a specific location. The SDR estimates the extent to which SE is stored in an area. SDR is the ratio of SY to SE of a basin. The SDR was calculated using empirical relation developed by Sharda and Ojasvi 2016 (Equation (16)) [64]. The relation was modified by Swaarnkar et al., 2018 (Equation (17)) [48]. and where A is the area of basin/watershed/catchment. The Sanghar catchment was further sub-divided into smaller 133 sub-catchments according to the drainage density of the area. The SY was calculated at the sub-catchment level. The division into smaller sub-catchment allows the SY estimates at a different location using SDR and SE of individual sub-catchment. The process was followed to calculate the SY at the proposed dam sites.
The area under study is an ungauged catchment with no sediment measurements. The methodology was verified with the adjacent catchment of Vidor. A project on Vidor torrent was initiated in 2016 for sediment measurements by the Punjab Irrigation Department, Pakistan. With the limited number of records, the SY results were compared with the observed SY records.

Dam Suitability
The selected indicators show the skewed distribution for 269 box units along the mainstream. The positively skewed nature of indicators restricts the recurrence of high values in the catchment. The mean and standard deviation of all indicators substantiates high frequency of low to medium range of values. The skewed nature of indicators was observed due to the high ridges and heterogeneity in the area. The statistics of all indicators are shown in the Table 5. Based on the defined relations of indicators followed by inter and intra weights incorporated in the SAW technique, two sites were identified as most suitable for flash flood water management. The sites were selected from five suitable combinations of reservoir and dam wall as shown in Figure 3.
As per the Strahler classification scheme, the Sanghar hill torrent is composed of five stream orders with 5th as the main channel. Dam-A is located upstream of Dam-B on the fifth ordered stream (first criteria). The two suggested sites are approximately 22 km apart. The proposed site for Dam-B is considered as the main dam with a much higher capacity as compare to Dam-A. The reservoir areas of selected sites contain high to very high values of the MRVBF index (second criteria). However, TRI results substantiate low values in the reservoir area and high values for the dam wall (third criteria). The behavior of TRI shows that reservoir comprises of homogeneous topography as shown in Figure 4.
The reservoir areas of both sites show very low value (<5%) along the reservoir and high to very high values (100% < slope < 120%) for the dam wall respectively (fourth criteria). The sites were selected with large catchment areas to accumulate maximum flash flood water. Among the 269 box units on main order stream 25, and 60 area classified in very high, and high classes respectively for the reservoir suitability. However, 18, and 60 box units are classified in very high, and high class respectively for dam wall suitability. The topographical summary of the proposed sites is shown in Table 6.
CD, VD, and contour along the main channel were analyzed to identify the best gorge in the area (fifth, and sixth criteria). A gorge is classified as a narrow valley with steep cross slopes located between hills or mountains. The selected sites represent a narrow opening with a large reservoir area having the capability to accumulate maximum water with minimum dam wall construction. This will allow the selected sites to accumulate maximum water with narrow opening. The mountainous topography of the area substantiates natural existence of perfect gorges along the main ordered stream. The cross-sectional profile of the dam's wall is shown in Figure 5.  As per the Strahler classification scheme, the Sanghar hill torrent is composed of five stream orders with 5th as the main channel. Dam-A is located upstream of Dam-B on the fifth ordered stream (first criteria). The two suggested sites are approximately 22 km apart. The proposed site for Dam-B is considered as the main dam with a much higher capacity as compare to Dam-A. The reservoir areas of selected sites contain high to very high values of the MRVBF index (second criteria). However, TRI results substantiate low values in the reservoir area and high values for the dam wall (third criteria). The behavior of TRI shows that reservoir comprises of homogeneous topography as shown in Figure 4.  The reservoir areas of both sites show very low value (<5%) along the reservoir and high to very high values (100% < slope < 120%) for the dam wall respectively (fourth criteria). The sites were selected with large catchment areas to accumulate maximum flash flood water. Among the 269 box units on main order stream 25, and 60 area classified in very high, and high classes respectively for the reservoir suitability. However, 18, and 60 box units are classified in very high, and high class respectively for dam wall suitability. The topographical summary of the proposed sites is shown in Table 6. CD, VD, and contour along the main channel were analyzed to identify the best gorge in the area (fifth, and sixth criteria). A gorge is classified as a narrow valley with steep cross slopes located between hills or mountains. The selected sites represent a narrow opening with a large reservoir area having the capability to accumulate maximum water with minimum dam wall construction. This will allow the selected sites to accumulate  maximum water with narrow opening. The mountainous topography of the area substantiates natural existence of perfect gorges along the main ordered stream. The cross-sectional profile of the dam's wall is shown in Figure 5.

EAC Curve
The contour analysis of suitable sites was performed for EAC curve formation. The surface area and capacity at individual contour elevation was calculated. The selected site for dam-A consists of a maximum of 39 contours with a 3 m interval (39 × 3 = 117 m) having a base contour elevation of 513 m.

EAC Curve
The contour analysis of suitable sites was performed for EAC curve formation. The surface area and capacity at individual contour elevation was calculated. The selected site for dam-A consists of a maximum of 39 contours with a 3 m interval (

EAC Curve
The contour analysis of suitable sites was performed for EAC curve formation. The surface area and capacity at individual contour elevation was calculated. The selected site for dam-A consists of a maximum of 39 contours with a 3 m interval (

R-Factor
The effect of the R-factor for the estimation of SE is averaged out over the selected period. The areas with a higher altitude, particularly mountainous catchments, receive high rainfall. The study area is a part of the mountainous range with varied topography from upstream to downstream. The R-factor shows that the rainfall rate is higher around the ridge in the center of the catchment with a height > 2000 m. The value of the R-factor ranges from 1073.75 to 2418.49 from low to high, respectively. Low rainfall values show

R-Factor
The effect of the R-factor for the estimation of SE is averaged out over the selected period. The areas with a higher altitude, particularly mountainous catchments, receive high rainfall. The study area is a part of the mountainous range with varied topography from upstream to downstream. The R-factor shows that the rainfall rate is higher around the ridge in the center of the catchment with a height > 2000 m. The value of the R-factor ranges from 1073.75 to 2418.49 from low to high, respectively. Low rainfall values show low erosivity in the eastern part of the catchment. The R-factor gradually reduces in the downstream relatively low elevated areas. The R-factor varies considerably for different storm intensities.

LS-Factor
Due to the mountainous nature of the catchment, LS-factor shows a high influence on SE results. The results show higher values of cells on the ridge and adjacent areas (13 ≤ LS ≤ 20.5). The higher values for ridge areas can be due to the steep slope (S-factor). However, relatively higher values of LS-factor were also observed for the cells close to the main channel (fifth order). This is mainly due to the large values of stream accumulation/contributing area (L-factor). Low values of LS-factor (<2) in the patchy pattern were observed for the rest of the catchment. The variation for the other parts of catchments was observed between 2 and 12.

K Factor
The soil map shows the presence of clay loam in the eastern part and patches of sandy clay loam in the catchment; the rest of the catchment is mainly covered with loam. Typically, the loamy soil is considered more susceptible to erosion as compared to the other classes. K-factor calculation also involves the percentage of OM content in the soil. OM increases soil fertility and stability. The OM results show a patch of higher values (2 ≤ OM% ≤ 3.57) in the north of the catchment, however, small values of OM (<1) exist in patches; the rest of the catchment contains a spatially variable range of OM between 1 and 2. The percentage silt, sand, clay and OM content results in K-factor ranges from 0.016 to 0.034. The area is mostly of loam soil, spatial variation of high silt content along the loam soil reflects the high values of K-factor.

C Factor
The LULC suggested C-factor values mentioned in Table 4 show that 1 for bare-land and 0.97 for desert plain covers 1826.7 km 2 and 1209.2 km 2 areas, respectively. Higher values of the C-factor increase the vulnerability to soil erosion. However, the rate of soil erosion decreases with vegetation cover. The vegetation cover improves the physical and chemical properties with the increase in the percentage of OM and provides the barrier to eroded soil. Herbaceous cover (C-factor = 0.7) and shrubland (C-factor = 0.6) covers 673.4 km 2 and 985.8 km 2 area of the Sanghar catchment, respectively.

Soil Erosion (SE)
The results of different factors are combined using the RUSLE equation resulting in an estimation of SE in tons per hectare per year (t-ha −1 y −1 ). The skewed nature of input datasets including LS, OM, and slope tends the SE output in a positively skewed distribution. The generated raster has a mean of 75 t-ha −1 y −1 with a high standard deviation of 108.27. The minimum and maximum values of SE are 0 and 799 t-ha −1 y −1 , respectively.
The SE has spatial variability in low values (SE < 50) t-ha −1 y −1 due to the variability of terrain. The low range covers the maximum of the Sanghar catchment. Higher values (SE > 373) t-ha −1 y −1 of SE are observed in the central part of the catchment. The clustered high values of SE in the center of the Sanghar catchment is due to the concentration of high values of R, LS, and K-factors and associated heterogeneity in elevations. SE values (211 < SE < 373) t-ha −1 y −1 are trending in the adjacent cells of higher SE values. However, the value range of (50 < SE < 211) t-ha −1 y −1 covers the rest of the catchment as shown in Figure 8. The generated SE substantiates among the 8 defined ranges, the majority of the area lies in the range (SE < 50) t-ha −1 y −1 as shown in Figure 8. The average value of SE at individual sub-catchment shows strong positive relation with LS-factor. SE values also show a positive relationship with R-factor; however, the relationship is weak as compared to LS-factor. The soil composition also relates to the SE estimates. The average values of SE increase with the increase in percentage sand concentration and decrease with an increase in percentage silt concentration as shown in Figure  9. These two compositions of soil play a significant role in the calculation of the K-factor. The average value of SE at individual sub-catchment shows strong positive relation with LS-factor. SE values also show a positive relationship with R-factor; however, the relationship is weak as compared to LS-factor. The soil composition also relates to the SE estimates. The average values of SE increase with the increase in percentage sand concentration and decrease with an increase in percentage silt concentration as shown in Figure 9. These two compositions of soil play a significant role in the calculation of the K-factor.

Sediment Yield and Dam Life
The soil eroded from a catchment due to rainfall is termed as SY. The SY estimates at dam locations solely depend upon the empirical relation of SDR. For both dam locations, the two SDR equations show uncertainties. The SDR relation developed by Swarnkar et al. [64] gives better estimates of SY. The SDR overestimates the SY results for a larger catchment in the area. However, it was observed that the relation behaves well at the sub-catchment level. Based on records of rainfall, the annual average sediment of 298,073 tons and 318,000 tons will feed Dam-A and Dam-B, respectively.
According to the ICOLD [65] classification, a small dam has a height (2.5 m < H < 15 m) and volume of water stored in a reservoir (H 2 × √ V < 200 Mm 3 ). However, height restriction is not observed for dam classification in Pakistan, and dams of height (H > 40 m) have been constructed as small dams [66]. Based on usual practices in Pakistan lives of selected dams were calculated at 63 m and 48 m height for Dam-A and Dam-B, respectively. The result substantiates approximate lives of 87 and 90 years for Dam-A and Dam-B, respectively, as shown in Figure 10.

Discussion
The mountain range of Suleiman with widespread from northeast to southeast of Pakistan is prone to flash floods. The frequent occurrence of low to medium intensity flash floods in the monsoon season makes these catchments more vulnerable as compared to the other parts of Pakistan. The seasonal erratic patterns of rainfall further restrict the usual cultivation of different crops due to the water-scarce nature of the area. The historical cropping level of the hill torrents also indicates that only a small percentage of flood flows during high floods are utilized for agriculture and the remaining flood flows damage canals and canal command area.
Implementation of IMPs (integrated management practices) and BMPs (best management practices) using modern techniques of hydrologic investigations have proven to be a viable method for the management of these areas through dam site suitability and soil loss estimation.
The present study adopts a new methodology for dam site selection using different morphometric datasets including SO, MrVBF, TRI, Slope, CD, VD, and GD. The methodology proves a strong dimension in the dam site selection. The selected parameters/indicators reflect the maximum possible set of rules defined in the literature [3,4,18,67,68] for dam site suitability analysis. The index/indicator as a qualitative approach makes the methodology unique and reduce the user dependency for individual thematic layer generation. Validation of methodology was conducted with an existing dam (Gomal-Zem Dam) in the north of the Suleman mountain range. The adopted methodology identified the Gomal-Zem Dam in a very high suitable class. The designed methodology and selected sites are presented to Secretary and Additional Secretary PID (Punjab Irrigation Department), Pakistan, and found satisfactory. According to [62], both dams are classified in the small dam category. Several studies [3,4,67,68] identified suitable locations for the small dam, check dam, and percolation tank for water harvesting in Pakistan and India.

Discussion
The mountain range of Suleiman with widespread from northeast to southeast of Pakistan is prone to flash floods. The frequent occurrence of low to medium intensity flash floods in the monsoon season makes these catchments more vulnerable as compared to the other parts of Pakistan. The seasonal erratic patterns of rainfall further restrict the usual cultivation of different crops due to the water-scarce nature of the area. The historical cropping level of the hill torrents also indicates that only a small percentage of flood flows during high floods are utilized for agriculture and the remaining flood flows damage canals and canal command area.
Implementation of IMPs (integrated management practices) and BMPs (best management practices) using modern techniques of hydrologic investigations have proven to be a viable method for the management of these areas through dam site suitability and soil loss estimation.
The present study adopts a new methodology for dam site selection using different morphometric datasets including SO, MrVBF, TRI, Slope, CD, VD, and GD. The methodology proves a strong dimension in the dam site selection. The selected parameters/indicators reflect the maximum possible set of rules defined in the literature [3,4,18,67,68] for dam site suitability analysis. The index/indicator as a qualitative approach makes the methodology unique and reduce the user dependency for individual thematic layer generation. Validation of methodology was conducted with an existing dam (Gomal-Zem Dam) in the north of the Suleman mountain range. The adopted methodology identified the Gomal-Zem Dam in a very high suitable class. The designed methodology and selected sites are presented to Secretary and Additional Secretary PID (Punjab Irrigation Department), Pakistan, and found satisfactory. According to [62], both dams are classified in the small dam category. Several studies [3,4,67,68] identified suitable locations for the small dam, check dam, and percolation tank for water harvesting in Pakistan and India. For example, Jamali et al. 2014 [68] present the solution for site suitability for DAM in northern Pakistan using modern techniques of GIS and RS with spatial multi-criterion analysis (SMCA). Singh et al. 2008 [67] identify fourteen suitable sites for check dams using different thematic layers including; LULC, hydrologic soil group, slope, and DEM. The identification of suitable sites followed the guidelines provided by the Integrated Mission for Sustainable Development (IMSD). A latest research has been conducted in Iraq. This research utilized all possible variables including geological formation, soil type, fault line, tectonic line, altitude, slope, rainfall data, water discharge, land use/cover, road network, and material used for dam site selection [69]. Contrarily, very few attempts had been made which utilizes specific indices for dam site suitability analysis. However, the designed methodology grouped certain indices which represent multiple criteria for dam site suitability analysis. The set of indices set well on the varied topography of the areas. The selected indicators cover the maximum set of rules in the literature cited above. The relation between different indicators with the reservoir and dam wall suitability makes the methodology more robust. The selected indicators are generalized and hence could be implemented in any other area.
The efficacy of the designed structure is analyzed by annual average SE estimates. The estimation of SE was necessary because the mountain range of Koh e Suleman is prone to soil erosion due to its varied topography and barren nature of the land cover. The RUSLE equation involves R-rainfall, K-soil erodibility, LS-length slope, C-landcover management, and P-support practice factors to estimate the annual average soil loss of Sanghar torrent. The erratic behavior of rainfall mostly adopting sudden surge from cloud bursts and thunderstorm categorize the episodic nature of streams which remains active in monsoon season (July, August, and September) as shown in Table 2. The high variation in topography and sandy nature increase the soil loss probability in the study area. The SE results using the RUSLE equation shows an average soil loss of 75 t-ha −1 y −1 . The SY estimates using SDR and SE results show approximately 298,073 tons and 318,000 tons of annual average sediment yield will feed the dam-A and dam-B respectively. The SY results substantiate 87 y and 90 y lives of dam-A and dam-B at a height of 150 ft and 200 ft, respectively. The validation of methodology with limited records of Vidor shows sufficient accuracy of results. The average QOBSERVE and QRUSLE sediment record of Vidor torrent for the monsoon season of 2016 were 7925 t-day −1 and 7455 t-day −1 , respectively.
Focusing on the global climate change scenarios, if the frequency of events per year increases dramatically, it will increase the annual soil loss of the Sanghar catchment. The increased SE would decrease the associated dam life. The selected sites have the potential of upraising the dam height up to 385 ft and 325 ft for Dam-A and Dam-B, respectively. The Dam-B with a much higher capacity of 839 million m 3 is considered as the main dam for controlling the flash flood peaks. However, Dam-A could be beneficial in series to enhance the life of Dam-B. The strategy is flexible and could be implemented for climate change adaptation schemes. The enriched literature proves the applicability and efficacy of the RUSLE equation. The latest research was conducted by [23] for soil loss estimation using the RUSLE model on the Chitral river catchment in the North of Pakistan. Another study on Rawal Dam catchment in the North of Punjab province, Pakistan, was conducted by [70] and found GIS and RS techniques suitable for SE estimates (R2 = 0.76). Kusimi et al. 2015 [71] studied the effect of uncontrolled land use activities on SE in the Pra River, Ghana using the RUSLE equation. Djoukbala et al. 2018 [72] also evaluate the reliability and effectiveness of the RUSLE equation with the integration of GIS in Wadi El-Ham Algeria and found satisfactory results.
Pakistan is facing a shortage of water for several years, primarily because of the increase in population and mismanagement of accessible water resources. The water war scenario with neighboring countries accentuate the existing drastic nature of water shortage. There is average annual potential of 23 billion m 3 in D G Khan's hill torrent water resources of the country that has not yet been used to its productive potential [73]. To overcome the water shortage of the country climate change adaptation strategies are by far the most suitable practices that has been implemented around the globe. The structural management of flash flood by dam construction also fulfils the crop water requirement (CWR) in the piedmont plains of the area.

Conclusions
The study highlights the existing potential of hill torrents in Pakistan. The results promote the use of spate irrigation to improve the cropping intensity of the area. According to Ahmad et al., 2016 [73], the benefit cost ratio of spate irrigation for piedmont plains of D G Khan is higher than canal and ground water. Therefore, the study concludes the efficient utilization of hill torrents for spate irrigation through the construction of small dam that would manage the flash flood water in environmentally friendly way. Further studies should incorporate the CWR calculation which would help in the modelling, demarcation and management of cultural command areas in the downstream of selected dam sites. The adopted methodology confirms the applicability of IMPs in flash flood and erosion prone areas for the better management of natural hazards. The generated results substantiate effectiveness and reliability of RDSA technique for dam site suitability analysis and RUSLE for annual SE estimates. The study concludes the following testimonials: 1.
RDSA technique classified total of 269 box-units in very high to very low classes for reservoir and dam wall suitability. The changing behavior of climate over few decades causes flash flood scenarios more frequent in Pakistan. The proposed sites have the potential to enhance the storage capacity of dams to overcome the long-term forecasted climate change scenarios. The dam wall height of proposed sites can be increased to maximum of 384 ft and 325 ft for Dam-A and Dam-B, respectively.