Topographic Factor (LS)

The topographic factor *LS* reflects the effect of slope length and steepness on soil erosion from the catchment. It considers sheet, rill, and inter-rill erosion by water, and is a ratio of expected soil loss from a field slope relative to the original USLE unit plot [32]. In RUSLE, this relationship was extended to the more complex hill slopes for better estimation of the topographic effect by modifying the steepness factor that is more sensitive than the length factor [1]. Some studies extend the *LS* factor to topographically complex terrains using a method that incorporates contributing area and flow accumulation [35]. The flow accumulation tool built into a GIS allows the estimation of the upslope contributing area using a digital elevation model (DEM) that, in turn, is used to calculate *LS* factor to account for more topographically complex landscapes [35,36]. The flow accumulation method to calculate slope length and steepness explicitly accounts for convergence and divergence of flow, which is important while considering soil erosion over a complex landscape [37]. However, a high-resolution DEM is needed for an accurate representation of the topography to calculate *LS* factor. Benavidez et al. [25] reported that the method of flow accumulation performed significantly better at the sub-watershed or field scale. We have used Equation (5) in this research to determine the *LS* factor, developed by [36] and widely used globally as well as for regional studies.

$$L\mathcal{S} = (m+1)(\frac{\mathcal{U}}{L\_0})^m(\frac{\operatorname{Sim}\mathcal{\beta}}{\mathcal{S}\_0})^n\tag{5}$$

where *U* is upslope contributing area per unit width (m2/m), *L*<sup>0</sup> is the length of the unit plot (m), *S*<sup>0</sup> is the slope of the unit plot (%), *β* is the slope (%), *m* is a factor for sheet erosion and *n* is a factor for rill erosion. The value of *m* ranges from 0.4 to 0.6 whereas that of *n* ranges from 1.0 to 1.3, depending on the prevailing erosion types. For this study, the value of *m* and *n* were used as 0.4 and 1.3, respectively as suggested by [36]. *L*<sup>0</sup> and *S*<sup>0</sup> are constants, and their values are 22.1 m and 9% in Equation (5).

The DEM of the Shuttle Radar Topography Mission (SRTM) with a 30 m resolution was used to determine catchment and sub-catchment boundaries and calculate slope length and slope steepness which are used in the RUSLE topographic factor *LS* (https: //earthexplorer.usgs.gov/, accessed on 6 May 2019). ArcGIS tool was used for spatial analysis for *LS* factor and catchment delineation using appropriate tools.

Crop Management Factor C and Conservation Support Practice Factor P

The crop management factor *C* reflects the effects of cropping and management practices on erosion. In RUSLE, its value is calculated as a product of five subfactors: prior land use, canopy cover, surface cover, surface roughness, and soil moisture. Factor *C* is expressed as a ratio of soil loss of a parcel with certain land use and a fallow condition [32,38]. Generally, *C* ranges between 1 and 0; higher values (≤1) are used for less vegetation, thus considering the land as barren whereas lower values (near to zero) indicate very strong cover and well-protected soil. Lanorte et al. [39] have reported that many authors have adopted simplified approaches to estimate *C,* e.g., by using land cover maps and assigning a *C* value to each class [40] or by applying remote-sensing techniques such as image classification [41,42] and vegetation indices [43]. Various studies [44–46] report mathematical functions to calculate *C* using the normalized difference vegetation index (NDVI), which is positively correlated with the amount of green biomass and indicates differences in green vegetation coverage [45]. In this study, the mathematical function developed by [45] has been used to estimate the crop management factor as shown in Equation (6). The NDVI calculations were made from Band 5 (NIR) and 4 (Red) of Landsat 8 image (https://earthexplorer.usgs.gov/, accessed on 29 May 2019) for July 2018, as shown in Figure 3. Most of the precipitation in the study area is received during the monsoon season (July to September) and the scattered vegetation dependent on the monsoon rainfall. Therefore, single imagery was used to calculate NDVI; otherwise, average values of different periods may be used where different cropping patterns prevail, as suggested by [13].

$$\mathcal{C} = -\exp[I(\text{NDVI}/r - \text{NDVI})] \tag{6}$$

where *C* is crop management factor, NDVI is normalized difference vegetation index, *l* and *r* are constants, and their values are taken as 2 and 1, respectively.

Conservation support practice factor *P* represents the ratio of soil loss for a unit area with a given erosion control method compared to the soil under the same conditions without any erosion control measures [25,26]. The value of the *P* factor ranges from 0 to 1 (0 indicates good conservation practice and 1 represents poor conservation practice). In this research, *P* was set equal to one across the sub-catchment because there are no erosion-control works in the sub-catchment to prevent soil erosion.

Finally, the annual soil erosion rate was calculated by applying RUSLE. Morgan [26] developed a coding system for soil erosion appraisal in the field where classes have been defined based on erosion rate. We use the coding system by [26] to explain erosion rates calculated in the sub-catchment.

geometry.

Multiblock option (SSIIM) 3D model [47] is used for sediment transport modeling to compute sediment trapped in the hypothetical settling reservoir. The simulation informs how efficient this structure could be in reducing sediment entry into the canal irrigation system. The SSIIM program solves the Navier Stokes equations over a structured mesh, using the k-epsilon model to represent the turbulent viscosities. The main strengths of SSIIM are the capability of modeling sediment transport with a movable bed in complex

Once the flow field is computed, the convection-diffusion 3D sediment equations

**Figure 3.** Normalized difference vegetation index (NDVI) for the sub-catchment.

### **Figure 3.** Normalized difference vegetation index (NDVI) for the sub-catchment. *2.3. Sediment Transport Modeling*

For the hypothetical settling reservoir, a grid of 22 × 6 × 6 (792 cells) was used. Inflow into the basin was taken as 24 m3 s−1; sediment discharge was considered as 10, 156, and 151 kg s−1 for sand silt and clay particle sizes, respectively. The water level of 2 m from the bed of the settling reservoir was considered as an outflow boundary condition. The flows and sediment concentration data of 2016 was used to calculate the amount of sediment deposition in the proposed settling reservoir. For this purpose, a sediment rating curve was developed to calculate the sediment concentration during missing days (Figure 4). The amount of bedload in river flows during rainfall events assumed as 10% of suspended sediment load as suggested by [48]. For improved estimates of the sediment load, About 70% of annual sediment load is diverted into the canals only during the *monsoon* season [9]. In flashy streams, significant proportions of sediment move as bedload, which is challenging to account for with conventional sampling methods for suspended sediments. The bed load fraction creates more problems when it enters irrigation canals due to the limited sediment transport capacity of the irrigation canals [10]. In our case, sediment enters into the irrigation system through a diversion barrage build across the Gomal river. The description of the diversion barrage is provided in the Study Area section. In this research, we simulated flows in a settling reservoir proposed immediately downstream of the diversion barrage and the main canal feeding from this settling reservoir rather than from a head regulator at the diversion barrage. The size of this hypothetical settling reservoir is taken as 100 <sup>×</sup> <sup>80</sup> <sup>×</sup> 2 m<sup>3</sup> . Sediment Simulation In Intakes with Multiblock option (SSIIM) 3D model [47] is used for sediment transport modeling to compute sediment trapped in the hypothetical settling reservoir. The simulation informs how efficient this structure could be in reducing sediment entry into the canal irrigation system. The SSIIM program solves the Navier Stokes equations over a structured mesh, using the k-epsilon model to represent the turbulent viscosities. The main strengths of SSIIM are the capability of modeling sediment transport with a movable bed in complex geometry.

> Once the flow field is computed, the convection-diffusion 3D sediment equations are solved for each of the sediment fractions considered. From these, the trapping efficiencies, as well as the depositional sediment pattern, are obtained.

> For the hypothetical settling reservoir, a grid of 22 × 6 × 6 (792 cells) was used. Inflow into the basin was taken as 24 m<sup>3</sup> s −1 ; sediment discharge was considered as 10, 156, and

151 kg s−<sup>1</sup> for sand silt and clay particle sizes, respectively. The water level of 2 m from the bed of the settling reservoir was considered as an outflow boundary condition. The flows and sediment concentration data of 2016 was used to calculate the amount of sediment deposition in the proposed settling reservoir. For this purpose, a sediment rating curve was developed to calculate the sediment concentration during missing days (Figure 4). The amount of bedload in river flows during rainfall events assumed as 10% of suspended sediment load as suggested by [48]. For improved estimates of the sediment load, the continuous records approach applied by [49] can be used, but it requires the hourly or daily sediment concentration data. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 12 of 23 the continuous records approach applied by [49] can be used, but it requires the hourly or daily sediment concentration data.

**Figure 4.** Sediment rating curve of the Gomal river at monitoring station during the year 2016. **Figure 4.** Sediment rating curve of the Gomal river at monitoring station during the year 2016.

### **3. Results and Discussion**

### **3. Results and Discussion**  *3.1. Annual Sediment Yield*

*3.1. Annual Sediment Yield*  Figure 5 shows the cumulative sediment yield at the monitoring station near the diversion barrage (see location in Figure 2). At this location, the sediment yield includes sediment contribution from the entire GZ catchment (including the sub-catchment). Before the GZ dam construction in 2011, the sediment load from the entire catchment was contributing to the monitoring station. There was no canal irrigation system, and the community had historic rights to use the river water through spate irrigation. Flow and sediment monitoring remained discontinued from 1996 to 2002. The trend lines for the periods of 1960–2010 and 2003–2010 are almost parallel, showing no substantial difference in the pattern of sediment load in the Gomal river. The trend line for the period of 2011–2015 shows a significant change in the sediment load. This is intuitive because the GZ reservoir has been filling during this period, and the sediment from the upstream catchment used to be trapped in the reservoir. The trend line pre and post-dam constructions show an 85% decrease in the annual sediment yield, suggesting that significant sediment yield is still contributing to the diversion barrage although the sub-catchment area comprised 1.3% of the total catchment area. This high sediment contribution can be attributed to the sediment releases from the dam during the flushing operation, and the Figure 5 shows the cumulative sediment yield at the monitoring station near the diversion barrage (see location in Figure 2). At this location, the sediment yield includes sediment contribution from the entire GZ catchment (including the sub-catchment). Before the GZ dam construction in 2011, the sediment load from the entire catchment was contributing to the monitoring station. There was no canal irrigation system, and the community had historic rights to use the river water through spate irrigation. Flow and sediment monitoring remained discontinued from 1996 to 2002. The trend lines for the periods of 1960–2010 and 2003–2010 are almost parallel, showing no substantial difference in the pattern of sediment load in the Gomal river. The trend line for the period of 2011–2015 shows a significant change in the sediment load. This is intuitive because the GZ reservoir has been filling during this period, and the sediment from the upstream catchment used to be trapped in the reservoir. The trend line pre and post-dam constructions show an 85% decrease in the annual sediment yield, suggesting that significant sediment yield is still contributing to the diversion barrage although the sub-catchment area comprised 1.3% of the total catchment area. This high sediment contribution can be attributed to the sediment releases from the dam during the flushing operation, and the sub-catchment as well during rainfall events.

Figure 6 presents a sediment rating curve for the selected years from 1994 to 2016 due to the availability of limited data. Suspended sediment concentration observed at the monitoring station ranges from 11 to 82,500 parts per million (ppm). The data after dam

toring station. Figure 7 presents the histogram of discharge and sediment concentration during 2015–2016. In Figure 7a, it is evident that the discharge remained 82% of the time around 24 m3s−1, which is about one-quarter of the irrigation demand (indent) of 100 m3s−1. On the other hand, Figure 7b shows that 50% of the time, sediment concentration

sub-catchment as well during rainfall events.

remained in the zero to 250 ppm range.

**Figure 5.** Cumulative sediment yield from Gomal river catchment during 1960–2015. **Figure 5.** Cumulative sediment yield from Gomal river catchment during 1960–2015.

A few high sediment concentrations were also observed at high discharges, i.e., >100 m3s−1, as shown in Figure 6. The figure shows high sediment concentration in 2016 as large red circles which correspond to the days when rainfall was observed at the Tank meteorological station. This confirms that the sub-catchment contributes high sediments in response to rainfall events as the sediment concentration was observed more than 4000 ppm during all rainfall events. However, sediment concentrations ranging from 61 to 6540 ppm in the same year 2016 (shown in green color) are clustered corresponding to 24 m3s−1. This sediment concentration is likely to be coming from water releases from the GZ dam, which is controlled to 24 m3s−1. It suggests that the dam also contributes to high Figure 6 presents a sediment rating curve for the selected years from 1994 to 2016 due to the availability of limited data. Suspended sediment concentration observed at the monitoring station ranges from 11 to 82,500 parts per million (ppm). The data after dam operation (2015–2016) gives more insight into the sediment and discharge at the monitoring station. Figure 7 presents the histogram of discharge and sediment concentration during 2015–2016. In Figure 7a, it is evident that the discharge remained 82% of the time around 24 m<sup>3</sup> s −1 , which is about one-quarter of the irrigation demand (indent) of 100 m<sup>3</sup> s −1 . On the other hand, Figure 7b shows that 50% of the time, sediment concentration remained in the zero to 250 ppm range. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 14 of 23

**Figure 6.** Sediment rating curve of the Gomal river at the monitoring station. **Figure 6.** Sediment rating curve of the Gomal river at the monitoring station.

Frequency (%)

 (**a**) (**b**) **Figure 7.** Histogram for data during 2015–2016; (**a**) discharge; (**b**) sediment concentration.

**Table 2.** Rainfall-runoff erosivity factor *R* using different approaches.

The European Soil Data Centre (ESDAC), Global dataset 538–1300 919

Rainfall-runoff erosivity factor *R (*MJ mm ha −1 h −1 yr −1)

0.0

Frequency (%)

81.6

2.1

0.0

7.1

5.0

0 25 50 100 200 400 600 800 1000 More

Discharge (m3s-1)

3.5

0.7

0.0

0.0

Rainfall at Tank station (mm) 226 212 134 191

Roose (1975) [27], Africa 113 106 67 95 Fernandez (2003) [29], USA 356 281 −128 170 İrvem (2007) [30], Turkey 412 531 272 405 Nakil (2014) [31], India 1006 994 934 978

Arnoldus (1980) [28], West Africa 38 59 8 35 Arnoldus (1980) [28], Western USA 37 58 7 34 Arnoldus (1980) [28], Northwest USA 22 25 18 22

 **Annual Data Average 2014 2015 2016** 

0.0

48.8

22.1

12.8

10.5

0 250 500 1,000 5,000 10,000 20,000 More

Sediment Concentration (ppm)

1.2

1.2

3.5

10

100

1,000

Sediment Concentration (ppm)

10,000

100,000

**Figure 6.** Sediment rating curve of the Gomal river at the monitoring station.

1 10 100 1,000 10,000

2015 2016 2016\_rainfall\_Tank Power (2010) Power (2007) Power (1994)

Discharge (m3s-1)

2010 2007 1994

**Figure 7.** Histogram for data during 2015–2016; (**a**) discharge; (**b**) sediment concentration. **Figure 7.** Histogram for data during 2015–2016; (**a**) discharge; (**b**) sediment concentration.

**Table 2.** Rainfall-runoff erosivity factor *R* using different approaches.  **Annual Data Average 2014 2015 2016**  Rainfall at Tank station (mm) 226 212 134 191 Rainfall-runoff erosivity factor *R (*MJ mm ha −1 h −1 yr −1) Arnoldus (1980) [28], West Africa 38 59 8 35 Arnoldus (1980) [28], Western USA 37 58 7 34 Arnoldus (1980) [28], Northwest USA 22 25 18 22 Roose (1975) [27], Africa 113 106 67 95 Fernandez (2003) [29], USA 356 281 −128 170 İrvem (2007) [30], Turkey 412 531 272 405 Nakil (2014) [31], India 1006 994 934 978 The European Soil Data Centre (ESDAC), Global dataset 538–1300 919 A few high sediment concentrations were also observed at high discharges, i.e., >100 m<sup>3</sup> s −1 , as shown in Figure 6. The figure shows high sediment concentration in 2016 as large red circles which correspond to the days when rainfall was observed at the Tank meteorological station. This confirms that the sub-catchment contributes high sediments in response to rainfall events as the sediment concentration was observed more than 4000 ppm during all rainfall events. However, sediment concentrations ranging from 61 to 6540 ppm in the same year 2016 (shown in green color) are clustered corresponding to 24 m<sup>3</sup> s −1 . This sediment concentration is likely to be coming from water releases from the GZ dam, which is controlled to 24 m<sup>3</sup> s −1 . It suggests that the dam also contributes to high sediment load during the flushing operation. Similar sediment concentration scatter was observed in the Upper Chenab Canal (UCC) and Marala Link Canal (MRLC) in Pakistan by [49], where average sediment concentrations were 450 and 500 ppm, respectively. However, the sediment concentration observed at the Gomal river monitoring station is much higher than UCC and MRLC. The above analysis is based on the suspended sediment concentration data only because the bed load at the monitoring station is not observed.

### *3.2. Erosion Rate Estimation*

Table 2 shows calculated values of rainfall-runoff erosivity factor *R* by using various empirical equations and the rainfall data on which calculations are based. The table also shows the range of *R* in the sub-catchment from the ESDAC global data set. The variation in calculated values of *R* is too wide to adopt a particular empirical equation. We have preferred the Roose equation [27] because researchers have used this equation to calculate *R* in recent studies to estimate soil erosion rate in Pakistan's Potohar catchments [34] and India's Pravara catchment [50]. The adopted value of *R* in this research is, therefore, 95 MJ mm ha−<sup>1</sup> h <sup>−</sup><sup>1</sup> yr−<sup>1</sup> (Table 2) which is based on rainfall data from the Tank meteorological station from 2014 to 2016.

Ideally, *R* is calculated as a sum of individual storm erosivity values (the rainfall intensity–kinetic energy relationship), as suggested by [1], whereby data of rainfall events should be of considerable length and high frequency. Such high temporal resolution rainfall (pluviograph) data needed for accurate computation of *R* is available only in a few regions of the world [51]. The ESDAC global dataset is based on the rainfall intensity–kinetic energy relationship and uses high-resolution rainfall data collected from 3540 stations across the globe [52]. The extensive list of stations, however, does not include many countries, and Pakistan is not an exception to this. In the case of the Gomal river sub-catchment, the limitation of low temporal resolution of rainfall data steered the analysis to consider empirical equations from regional studies. Table 2 shows significant variation amongst the computed *R* values using empirical equations and the global data set. The comparison in Table 2 is presented for the benefit of future studies in data-sparse catchments, where an empirical equation may yield a more realistic *R* than a value from the global data set. It may be particularly useful to consider empirical equations in those regions where high-resolution rainfall data was not available for global rainfall erosivity assessment by [52]. sity–kinetic energy relationship and uses high-resolution rainfall data collected from 3540 stations across the globe [52]. The extensive list of stations, however, does not include many countries, and Pakistan is not an exception to this. In the case of the Gomal river

sub-catchment, the limitation of low temporal resolution of rainfall data steered the

Ideally, *R* is calculated as a sum of individual storm erosivity values (the rainfall intensity–kinetic energy relationship), as suggested by [1], whereby data of rainfall events should be of considerable length and high frequency. Such high temporal resolution rainfall (pluviograph) data needed for accurate computation of *R* is available only in a few regions of the world [51]. The ESDAC global dataset is based on the rainfall inten-


**Table 2.** Rainfall-runoff erosivity factor *R* using different approaches. analysis to consider empirical equations from regional studies. Table 2 shows significant

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 15 of 23

The sub-catchment of the Gomal river is mostly hilly; thus, slopes vary from gentle to very steep. The average slope of the catchment is 33%; therefore, terrain may be considered steep. However, a significant portion of the sub-catchment exhibits less than a 40% slope. The topographic factor *LS* varies from 0 to more than 20, with a substantial proportion of the area (about 80%) having *LS* values less than 20 (Figure 8a). fore, it was used to calculate the topographic factor for this study. The computed crop management factor *C* (as shown in Figure 8b) is close to the published values reported by [34] for the regional catchment areas. The factor *C* > 0.83 indicates minimal vegetation cover, predominating in the sub-catchment and hence the soil is very vulnerable to erosion (Figure 8b).

**Figure 8.** Revised universal soil loss equation (RUSLE) factors estimation in the sub-catchment; (**a**) topographic factor *LS*; (**b**) crop management factor *C.*

Higher values >100 indicate more risk of soil erosion [53]. Moreover, the *LS* factor has been found to affect soil erosion rates more than any other RUSLE factor in the Appalachian hills of the United States [12]. In the Dobrov river catchment, Romania, researchers [54] have used different *LS* factors calculation methods and reported that the equations by [36] produced comparable erosion rates with the measured values; therefore, it was used to calculate the topographic factor for this study. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 16 of 23

The computed crop management factor *C* (as shown in Figure 8b) is close to the published values reported by [34] for the regional catchment areas. The factor *C* > 0.83 indicates minimal vegetation cover, predominating in the sub-catchment and hence the soil is very vulnerable to erosion (Figure 8b). **Figure 8.** Revised universal soil loss equation (RUSLE) factors estimation in the sub-catchment; (**a**) topographic factor *LS*; (**b**) crop management factor *C.*

> The results (Figure 9a) show that 46% of the sub-catchment area shows a very severe annual erosion rate (100–500 t h−<sup>1</sup> ) while 16% shows a catastrophic annual erosion rate (>500 t ha−<sup>1</sup> ). Very severe and catastrophic erosion rates have also been reported by different researchers for Pakistan, India, Ethiopia, and Brazil [34,55–57]. The main reasons for such high erosion rates are steep topography, sparse vegetative cover, and highly erodible soil type. The sub-catchment can be characterized as highly prone to erosion as compared with other catchments in Pakistan, e.g., the Simly dam catchment [58] or the Potohar Plateau [34,59]. Similarly, the soil erosion rates in the sub-catchment are much higher than reported catchments in the other region of the world, for example, in Nepal [60] and Ethiopia [61]. The results (Figure 9a) show that 46% of the sub-catchment area shows a very severe annual erosion rate (100–500 t h−1) while 16% shows a catastrophic annual erosion rate (>500 t ha−1). Very severe and catastrophic erosion rates have also been reported by different researchers for Pakistan, India, Ethiopia, and Brazil [34,55–57]. The main reasons for such high erosion rates are steep topography, sparse vegetative cover, and highly erodible soil type. The sub-catchment can be characterized as highly prone to erosion as compared with other catchments in Pakistan, e.g., the Simly dam catchment [58] or the Potohar Plateau [34,59]. Similarly, the soil erosion rates in the sub-catchment are much higher than reported catchments in the other region of the world, for example, in Nepal [60] and Ethiopia [61].

### *3.3. Scenario Testing in Revised Universal Soil Loss Equation (RUSLE) 3.3. Scenario Testing in Revised Universal Soil Loss Equation (RUSLE)*

Figure 8b shows computed *C* values in the Gomal river sub-catchment. The factor *C* can range from zero (for a non-erodible surface) to 1 for a bare plot (no vegetation). The sub-catchment represents mostly barren land with values >0.8. The minimum value of *C* in the sub-catchment is 0.3 that prevails sparsely in small pockets. We have tested a scenario to understand the impact of improved *C* on erosion rates. For this scenario, the value of *C* was set to 0.18 for the entire sub-catchment rather than computing *C* values with Equation (6) as shown in Figure 8a, and keeping all other factors of the RUSLE unaltered. The *C* used in the scenario is quite low, representing very well-covered and Figure 8b shows computed *C* values in the Gomal river sub-catchment. The factor *C* can range from zero (for a non-erodible surface) to 1 for a bare plot (no vegetation). The sub-catchment represents mostly barren land with values >0.8. The minimum value of *C* in the sub-catchment is 0.3 that prevails sparsely in small pockets. We have tested a scenario to understand the impact of improved *C* on erosion rates. For this scenario, the value of *C* was set to 0.18 for the entire sub-catchment rather than computing *C* values with Equation (6) as shown in Figure 8a, and keeping all other factors of the RUSLE unaltered. The *C* used in the scenario is quite low, representing very well-covered and managed land.

the European Union level and found the lowest mean values in Denmark and Hungary

In practical terms, achieving the *C* value of 0.18 in the sub-catchment by improving vegetation would be extremely difficult amid steep topography and less precipitation. However, we assume a meager value of *C* to see how well this potential catchment

management solution can work in reducing erosion rates.

(0.178, and 0.188, respectively).

To draw a parallel with the real cases, researchers [62] have mapped *C* at the European Union level and found the lowest mean values in Denmark and Hungary (0.178, and 0.188, respectively).

In practical terms, achieving the *C* value of 0.18 in the sub-catchment by improving vegetation would be extremely difficult amid steep topography and less precipitation. However, we assume a meager value of *C* to see how well this potential catchment management solution can work in reducing erosion rates.

The result in Figure 9b showed a decrease in the proportion of sub-catchment under catastrophic erosion rate, but the majority of the sub-catchment remained in the categories of high and severe erosion rates. A comparison of the erosion rates in various classes is shown in Table 3, which confirms that the crop management factor alone does not achieve a substantial reduction in soil erosion rates in the sub-catchment. Overall, the study area exhibits very high erosion rates. Global studies [5] show that management and the related land-use changes affect the spatial patterns and magnitude of accelerated soil erosion. The estimates of soil loss reduction derived from the implementation of conservation agriculture are encouraging in 40 countries, as reported by [5]. In contrast, the study area in this research has not shown a substantial reduction in soil erosion in response to a vegetation management scenario. These results also mirror the findings of [12], who reported that the Hocking river basin of the United States is covered with deciduous forests (about 50%) but is vulnerable to soil erosion due to steep slopes. Similarly, in a study on the Mangla River Basin [63], the authors reported that heavy rainfall increased the sediment load during the year 1992 despite increased vegetation cover compared to the previous years. Hence, vegetation alone in the study area would not help to significantly reduce erosion.


**Table 3.** Area of Gomal river sub-catchment in various erosion rate classes.

<sup>1</sup> Erosion rates categorization is used after Morgan [26].

Limiting management intervention to the sub-catchment can result in a partial reduction in sediments at the most, and even this may not turn out to be a viable option due to arid climatic conditions in this sub-catchment. However, the plantation campaign should be encouraged in the catchment upstream of the Gomal Zam dam without singling out any portion of the catchment that could be beneficial for environmental reasons. The geographic nature of the catchment is such that 30% of the catchment area is in Afghanistan, and two provinces share the part in Pakistan. Therefore, the catchment management option would be more effective if the two catchment-sharing countries jointly prepare a cooperative catchment management plan.

Moreover, check-dams can reduce the sediment load in the river basins with high erosion rates, as suggested by [64]. Check-dams in gullies could be one of the more effective ways to conserve soil and water, as observed in the Loess Plateau of China, where afforestation methods were not successful due to the arid climate and barren soils [65].

However, a few torrential storms could fill up the check dams due to high sediment load, as concluded by [66]. Also, the terracing at the catchment level could reduce the sediment yield by 4 to 5 times [14]. However, the use of appropriate management methods requires that whether the soil and climate conditions are feasible for the area or not. It could be inferred from the above discussion that any management option at the catchment scale would be challenging and have short-term impacts to control the sediment load in the river.

### *3.4. Sediment Transport in a Settling Reservoir*

A hypothetical settling reservoir provided immediately downstream of the existing diversion barrage is modeled. On the ground, a considerable area (256,691 m<sup>2</sup> ) is available downstream of the diversion barrage where a settling reservoir can be constructed. This area is currently not designated for productive uses, e.g., settlements, agriculture, etc. This makes the proposed settling reservoir a feasible option in the short term.

The sediment-laden flow would first enter into the settling reservoir where the sediments would settle down, and relatively clear water would outflow into the canal. The SSIIM model was set up using the hypothetical geometry and boundary conditions, as defined in the Data and Methods section. The results of the SSIIM model verify that the settling reservoir would be able to trap 16% of the incoming sediments. More precisely, it would be able to trap sand particles (coarse particles) almost entirely while the silt and clay would be trapped up to 25% and 0.5%, respectively. The overall trap efficiency is not as apparent as that of sand particles. However, the hopeful fact is that the settling reservoir is very efficient in preventing the entry of coarse particles, which are the most problematic particles among others. The model results reveal that 60,265 m<sup>3</sup> sediment will deposit in the settling reservoir every year that would otherwise have been entering the irrigation system.

The settling reservoir can provide an expeditious solution compared to catchment management options, e.g., improving vegetation, slope stabilization, etc. It would also minimize the dredging requirement from the pond area upstream of the existing diversion barrage. The modeling results also suggest that the settling reservoir should be operated such that a flow depth of 1.5 m should be maintained in the settling reservoir, which in turn would require regular dredging of the deposited sediments inside it.

The initial results of the model are encouraging and provide the basis for a quick solution. This modeling exercise was limited to only one geometry of the settling reservoir and a commonly prevailing flow condition. Prospective studies are needed to optimize settling reservoir size and use the transient sediment transport modeling approach.

### **4. Conclusions**

The irrigation system fed by the Gomal river is challenged by the massive amounts of sediment flowing into it. Limited sediment monitoring arrangement and hydroclimatic data make it a data-sparse catchment. Therefore, the sources of sediment are not clearly known. The soil erosion estimations and management option focuses on a sub-catchment of 450 km<sup>2</sup> (~1.2% of total catchment area) between the Gomal Zam dam and a diversion barrage ~40 km downstream of the dam.

The trends of cumulative sediment load pre and post-dam construction suggest an 85% decrease in the annual load. It can be inferred that the sediment share from the sub-catchment is ~15% that is larger compared to its relative size (1.3%) to the entire catchment. It means the sub-catchment could be considered as a hotspot within the Gomal river catchment. Moreover, results of the suspended sediment scatter also indicate that outflows during the flushing operation of the dam substantially contribute to the sediment load of the Gomal River that calls for a monitoring mechanism to compute the proportion of this contribution.

Soil erosion rates estimated using the revised universal soil loss equation (RUSLE) reveals that the different approaches and available data sets significantly impact the erosion rates. In data-sparse catchments like the Gomal river, the options become limited

to select a factor estimation method. In the case of rainfall-runoff erosivity factor *R,* a comparison of empirical equations and global data sets revealed colossal variation. The European Soil Data Centre (ESDAC) data set resulted in much higher values as compared with Roose's [27] relationship which is adopted for this study. In our case, the values of *R* adopted from a global dataset yield unrealistic erosion rates of catastrophic severity. Similarly, estimation of the cropping management factor with Knijff et al. [45] approach resulted in better values than the Durigon et al. (2014) [44] estimations. Therefore, the parameter estimation method for RUSLE empirical approach should be selected through careful consideration of regional literature rather than relying solely on global data sets, which in turn may give misleading results.

The results from RUSLE showed that most of the Gomal river sub-catchment area falls in very severe and catastrophic erosion rate categories (>100 t h−1y −1 ), and hence there are no apparent hotspots within the sub-catchment. High erosion rates can be attributed to steep topography, low vegetation, and highly erodible soil, as explained by RUSLE factors estimation.

A scenario was tested to understand the effectiveness of improved vegetation in the sub-catchment as an option to limit soil erosion. The results showed a small reduction in erosion rates. Half of the sub-catchment area still showed a high and severe erosion rate (10–100 t h−1y −1 ), and 17% of the area showed a very severe erosion rate (100–500 t ha−1y −1 ). This implies that improving vegetation in the sub-catchment can be beneficial for environmental reasons, but it may not be a single solution to prevent soil erosion from the catchment and consequently reducing sediment load in the river and that finds its way into the downstream irrigation system. For any catchment management option, e.g., improving vegetation, the entire catchment should be considered in the planning.

This study reveals that for catchments with such high erosion rates, complexities, and data limitations, improving one factor of RUSLE might not play a definitive role in managing land degradation. Sensitivity analysis of computed soil loss to changes in multiple factors of RUSLE and conditions could give better insight and poses an area of future work.

Reducing immediate negative impacts of high erosion rates remains a priority for catchment managers. In this context, the study finds a relatively quick solution by building a settling reservoir immediately downstream of the diversion barrage. The results from a 3D sediment transport model (SSIIM) suggest that 95% of sand particles and 25% silt particles can be trapped if the required water depth is maintained in the settling reservoir. The settling reservoir option, if adopted, would substantially reduce the sediment entry in the irrigation network. However, dredging would be required on an annual basis for the efficient operation of the settling reservoir.

Future research work should, therefore, seek to optimize the geometry of the settling reservoir and compare the results of various sediment transport models. It is also recommended to improve the sediment monitoring network at the Gomal river and the canal network for making informed management decisions. To better understand the sediment budget, future studies may use sediment yield assessment models for a better understanding of the sediment load contribution from the catchment area. The findings of this study are helpful for practitioners and to inform future investments to resolve the sedimentation problem. This research also provides the comparison of different methods to estimate the RUSLE factors applied in other parts of the world. The discussion around the selection of the methods can be helpful for the global audience to choose the appropriate methods of parameter calculation for catchments with similar characteristics.

**Author Contributions:** Conceptualization, M.T.B. and M.A.; Data curation, M.T.B.; Formal analysis, M.T.B. and M.A.; Funding acquisition, A.A.A.; Investigation, M.T.B. and M.A.; Methodology, M.T.B. and M.A.; Project administration, A.A.A.; Resources, A.A.A.; Software, M.A.; Supervision, A.A.A.; Validation, M.T.B. and A.A.A.; Writing—original draft, M.T.B.; Writing—review and editing, M.T.B., M.A. and A.A.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** U.S. Agency for International Development's Mission to Pakistan (USAID/Pakistan) under the terms of Award No. 72039118 IO 00003.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used for the study e.g., digital elevation model (DEM), satellite images, global rainfall-runoff erosivity factor, and soil type classification are available in the public domain and the web links are mentioned in the manuscript.

**Acknowledgments:** The International Water Management Institute (IWMI) is in receipt of financial support from the United States Agency for International Development (USAID) through Cooperative Agreement #72039118 IO 00003 which was used in part to support this study. The study design, data collection, analysis, and interpretation of the results are exclusively those of the authors and do not reflect the views or opinions of USAID, IWMI, or the CGIAR Research Program on Water, Land, and Ecosystems. The authors wish to acknowledge the contribution of Muhammad Humza in preparing the GIS maps used in this manuscript. The authors acknowledge NASA and USGS for providing Landsat images and SRTM DEM data free of cost. The authors also acknowledge the European Soil Data Centre (ESDAC) and FAO for providing free Global R and soil data, respectively.

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