Soil Erosion and Sediment Load Management Strategies for Sustainable Irrigation in Arid Regions

Abstract

Soil erosion is a serious environmental issue in the Gomal River catchment shared by Pakistan and Afghanistan. The river segment between the Gomal Zam dam and a diversion barrage (~40 km) brings a huge load of sediments that negatively affects the downstream irrigation system, but the sediment sources have not been explored in detail in this sub-catchment. The analysis of flow and sediment data shows that the significant sediment yield is still contributing to the diversion barrage despite the Gomal Zam dam construction. However, the sediment share at the diversion barrage from the sub-catchment is much larger than its relative size. A spatial assessment of erosion rates in the sub-catchment with the revised universal soil loss equation (RUSLE) shows that most of the sub-catchment falls into very severe and catastrophic erosion rate categories (>100 t h⁻¹y⁻¹). The sediment entry into the irrigation system can be managed both by limiting erosion in the catchment and trapping sediments into a hydraulic structure. The authors tested a scenario by improving the crop management factor in RUSLE as a catchment management option. The results show that improving the crop management factor makes little difference in reducing the erosion rates in the sub-catchment, suggesting other RUSLE factors, and perhaps slope is a more obvious reason for high erosion rates. This research also explores the efficiency of a proposed settling reservoir as a sediment load management option for the flows diverted from the barrage. The proposed settling reservoir is simulated using a computer-based sediment transport model. The modeling results suggest that a settling reservoir can reduce sediment entry into the irrigation network by trapping 95% and 25% for sand and silt particles, respectively. The findings of the study suggest that managing the sub-catchment characterizing an arid region and having steep slopes and barren mountains is a less compelling option to reduce sediment entry into the irrigation system compared to the settling reservoir at the diversion barrage. Managing the entire catchment (including upstream of Gomal Zam dam) can be a potential solution, but it would require cooperative planning due to the transboundary nature of the Gomal river catchment. The output of this research can aid policy and decision-makers to sustainably manage sedimentation issues in the irrigation network.

1. Introduction

Soil erosion in catchments occurs in various forms such as sheet, rill, gully, river bed and bank erosion, and landslides that contribute sediments to the water bodies. The rate of erosion is primarily determined by the erosive events (e.g., short duration and high-intensity rainfall events), soil type, and characteristics of the terrain [1]. The impacts of accelerated soil erosion processes can be severe, not only through land degradation and fertility loss but through a conspicuous number of off-site effects such as sedimentation, siltation, and eutrophication of waterways or enhanced flooding [2]. Soil erosion rates are exacerbated for the arid and semi-arid regions due to barren mountains with scattered vegetation that provide direct exposure to heavy rainfall. For example, 50% of the soil loss occurs in Pakistan during the monsoon season [3] causing huge sediment load diverted into the irrigation canals. Also, landslides and debris flow from the catchment increase the sediment load in the river flows. The eroded sediments finally deposit in reservoirs, stream channels, irrigation canals, and water conveyance structures and reduce the capacities of these water bodies to perform their prime functions and often requires costly treatments [1,4].

Soil erosion and the fate of eroded sediments are widely recognized as one of the major environmental concerns worldwide. The results of soil erosion rate for the year 2012 published in a global study [5] show that South America, Africa, and Asia had respectively 8.3%, 7.7%, and 7.6% of the continental area in high erosion rate classes (>10 t ha⁻¹yr⁻¹). In the global map of soil erosion rates in [5], the northern parts of Pakistan (where the study area of this research lies) fall in the range of vulnerable to high erosion rates (>50 t ha⁻¹yr⁻¹).

Sustainable irrigation requires appropriate soil erosion and sediment load management options and strategies both at the catchment level and the diversion structures. The catchment management options such as improved vegetation, slope stabilization, etc. are often considered preferred solutions to limit soil erosion rates as opposed to the structural measures for its control. However, in some cases, it becomes impractical to achieve anticipated benefits of the catchment management interventions amid factors such as; steep terrain, transboundary nature, urgency and severity of the challenges, etc. Another critical factor to consider is the feasibility of implementing such interventions.

Significant sediment entry from river flows may reduce the discharge capacity of irrigation canals up to 40% due to the settlement of coarse particles [6]. Many remedial measures have been suggested by different researchers to reduce the sediment load in canals, e.g., 41% using silt excluder [7], up to 40% using submerged vanes [8], and 40–45% using the vortex tube sediment extractor [9]. These structural measures only remove the coarser sediment in the canals; therefore, settling basins are assumed more appropriate options of sediment control for both coarser and finer particles [10]. Dredging of deposited sediment from irrigation networks also requires considerable effort and cost that farmers often have to put in. Similarly, the control of sediment load at the diversion structure depends on the feasibility of the proposed solution [11]. Hence, before making investment decisions, it is more practical to identify the targeted areas (i.e., hotspot analysis) in the case of opting for catchment management options. In the case of remedial measures (e.g., a hydraulic structure proposed to control sediment entry) the simulation model can be useful to predict the impacts of the proposed intervention on sustainable irrigation supplies.

Pakistan and Afghanistan share three main rivers, namely the Kabul, Kurram, and Gomal rivers. The Gomal river is the smallest among these three in terms of average annual inflow, but it is an essential source of livelihood for the downstream users in the Khyber Pakhtunkhwa (KP) province of Pakistan. The combined average flow of Kabul, Kurram, and Gomal rivers is 23.58 Gm³ yr⁻¹, out of which the Gomal river contributes 4.1% (0.974 Gm³ yr⁻¹). The physical infrastructures in the Gomal river catchment on the Pakistan side comprise a storage dam, two powerhouses, a diversion barrage, and a canal irrigation system to irrigate more than 77,000 ha of agricultural land.

Afghanistan is the upstream riparian on all shared rivers with Pakistan, including the Gomal river. The eroded sediment from the catchment in Afghanistan finds its fate in the reservoir of the Gomal Zam (GZ) dam and periodically sluiced in the downstream river reach. In the absence of a joint river management agreement between Pakistan and Afghanistan, integrated catchment management intervention cannot be introduced effectively. There is also no data-sharing mechanism between the two catchment-sharing countries.

The Gomal river brings a considerable amount of sediments to the diversion barrage, as shown in Figure 1. This sediment is a combination of soil eroded by rainfall in the sub-catchment below the GZ dam and the sediment generated in the catchment upstream of the GZ dam. Hydro-meteorological monitoring is very limited in the catchment area. Furthermore, the sediment flushing operation at the GZ dam is ad hoc, and its concentration in the outflows is not monitored. All this makes the catchment of the Gomal river a data-sparse catchment. When sediment-laden flow enters into the irrigation canal network, the sediments settle down in the channels and impede the conveyance capacity of the canals. This deprives the downstream users of their due share and gives rise to social unrest and conflicts.

The decisions to manage erosion rates in the catchment and eroded sediments in the water conveyance system are informed by the erosion/sediment assessment studies. These assessments help identify catchment areas with high soil erosion susceptibility and a clear understanding of the system’s operation. The revised universal soil loss equation (RUSLE) is commonly used to assess soil erosion rates. Soil loss predictions are frequently unrealistic because the methods used to estimate the soil loss equation’s factors are empirically derived. Therefore, the methodology section emphasizes the selection of suitable methods to determine RUSLE’s factors. Few studies have identified the impact of RUSLE factors on soil loss. For example, both ref [12,13] concurred in their findings that the topographic factor is the most significant factor among others, affecting the soil erosion rates. On the other hand, ref [14] has found that improving the support practice factor by adopting terracing showed a reduction in sediment yields. Previous studies show that no one rule fits all when investing in catchment management interventions.

Many studies that employ the RUSLE equation limit their analysis to improve the estimation of potential soil erosion. The input factors in RULSE can be estimated by using values from the literature or adapted for empirical and statistical data in combination with geographic information system (GIS) software [15], and this remains the focus of most of the studies.

This paper contributes to the science and application related to soil erosion as it is a critical problem worldwide, and especially in dryland areas. The paper first identifies the sources of sediment in a segment of the Gomal river (between GZ dam and the diversion barrage) in Pakistan. The incoming sediment limits supply in the irrigation system downstream and creates huge maintenance problems. A novel contribution of this study is its solution-oriented approach to the prevailing problem of sedimentation in the water conveyance system. There could be many possible solutions, but this paper has analyzed in detail a catchment management option (preventive) and a structural option (remedial). The analysis consists of erosion rates estimation and their severity level in the sub-catchment of the Gomal river using a computer-based empirical model. The study is focused on the sub-catchment because it is not transboundary and, therefore, management options can be introduced, monitored, and evaluated with relative ease. The analysis expands on predicting the impact of improved vegetation in the sub-catchment (i.e., crop management factor of RUSLE) as a catchment management option. Finally, a sediment transport model is used to determine the sediment trapping efficiency of a hypothetical settling reservoir at the diversion barrage.

2. Materials and Methods

2.1. The Study Area

Afghanistan and Pakistan share the catchment of the Gomal river. The Gomal river originates in the mountains of Ghazni province of Afghanistan and enters the South Waziristan district of KP, Pakistan. The total catchment area at Gomal Zam (GZ) dam is 33,950 km², about 70% of the area lies within Pakistan, and the remaining lies in Afghanistan.

The Government of Pakistan constructed the GZ dam (32°5′55″ N 69°52′53″ E) in South Waziristan district (formerly the South Waziristan Agency of the Federally Administered Tribal Area) with the financial assistance of the United States Agency for International Development (USAID). The construction of the dam started in August 2001 and was completed in April 2011. The dam impounds the Gomal River in a reservoir with 1.41 km³ storage capacity. The primary purpose of the dam is to supply water to irrigate 77,295 ha of land through a newly built irrigation system in the Tank and Dera Ismail Khan districts of Khyber Pakhtunkhwa province of Pakistan. The additional benefits of the dam are flood control and a modest amount of hydroelectric power generation (17.4 MW).

The Gomal river traverses around 40 km downstream of the dam site where a barrage is constructed near the town of Kot Murtaza, to divert the river flows into the irrigation system. The sub-catchment of this river segment (herein referred to as sub-catchment) is shown in Figure 2. The sub-catchment undergoes erosion and imparts sediments into the river stream, which find their way to the diversion barrage and enters into the canal irrigation system. This not only challenges the management of the irrigation system but also places significant demands on the operations and maintenance budget of the irrigation system.

Figure 2 shows the location and pictures of the GZ dam and diversion barrage. The diversion barrage is 140 m long and 7 m high. Three gates divert water into the main canal; in addition to that, four under sluice gates are also provided as shown in Figure 1a. The barrage is designed to pass a flood of 5200 m³s⁻¹ safely. The primary purpose of the GZ dam is to irrigate 77,295 ha of land in two districts of Khyber Pakhtunkhwa province i.e., D I Khan and Tank districts. For this purpose, an irrigation system consisting of the diversion barrage, 60.5 km long main canal, and 206 km of distribution canals have been constructed. The irrigation system is operational since November 2014, except for the portion fed by a branch canal named Waran Canal. A flow and sediment monitoring station is located approximately 2.5 km upstream of the diversion barrage, as shown in Figure 2.

At the diversion barrage, the river inflow is a combination of the flows released from the GZ dam and runoff from the sub-catchment. The area of this sub-catchment is ~450 km² constituting only 1.3% of the entire catchment of the Gomal river. This sub-catchment constitutes the study area for detailed analysis in this research.

The sub-catchment is mostly barren, with no essential utilities for human settlements. The community living in the command area of the irrigation system depends entirely on the water from the Gomal river not only for irrigated agriculture but in some parts also for domestic uses.

The existing dam operation is such that water is released from the tailrace of two hydroelectric power units at the GZ dam into the Gomal River. The main spillway has never been operated because water levels in the reservoir have not yet reached the spillway crest level. A low-level outlet is provided in the main dam to flush the sediments, which is operated on an ad hoc basis. This bottom outlet at an elevation of 680 m above mean sea level (dead storage level is 711 m) was designed for flushing out deposited sediments in the reservoir during monsoon, but one of the hydroelectric power units has never worked and, therefore, part of the water demand (8–10 m³s⁻¹) is supplied through this bottom outlet.

The mean annual rainfall in the sub-catchment is 276 mm yr⁻¹. The rainfall events are typically of high intensity, generate significant runoff, and bring with them eroded sediments. Erosion in this sub-catchment is often considered as a leading source of sediment [16] and, therefore, efforts for watershed management and conservation practices are emphasized as a potential solution to this issue of soil erosion [17]. The sediments from the catchment upstream of the GZ dam are trapped in the reservoir behind the dam and flushed through the low-level outlet. This sediment contribution is not quantified and monitored systematically; therefore, the only data available for analysis are at the monitoring station near the diversion barrage.

2.2. Data and Methods

2.2.1. Annual Sediment Yield in Gomal River

The long-term annual sediment yield data (1960–2015) was analyzed to find the cumulative trend at the monitoring station before and after the construction of the GZ dam. The monitoring station is located 2.5 km upstream of the diversion barrage site, and hence the sediment yield at this location includes the contribution from the entire catchment of the Gomal river, including the sub-catchment. The monitoring station is maintained by the Water and Power Development Authority (WAPDA) of Pakistan, which is responsible for the collection and analysis of sediment data. Moreover, it is noteworthy that only suspended sediment samples are collected by WAPDA using depth-integrated samplers (DH-48).

2.2.2. Assessment of Erosion Rates

The universal soil loss equation (USLE) was developed in 1965 [18] to estimate the rate of soil loss. This empirical equation was improved and revised in 1995, since then it has been known as the RUSLE [1] and has been applied in numerous catchments worldwide. RUSLE is not intended to estimate the sediment yield (as the amount of sediment reaching or passing a point of interest in a given period) instead; its application is limited to calculate the annual soil loss rate. However, the sediment delivery ratio (SDR) method can be used to link soil erosion with sediment delivery to the stream. With the RUSLE model, the average annual rate of soil loss for a catchment of interest can be predicted for any number of scenarios in association with cropping systems, management techniques, and erosion control practices [19]. Therefore, many studies have applied RUSLE to estimate soil loss rates and identification of high-risk areas, e.g., [20,21,22]. It has become the most widely used approach during an 80-year history of erosion modeling applied in 109 countries [23]. Alewell et al. [23] point out the limitations of USLE-type modeling that it does not address larger rills or gully erosion but is restricted to sheet/interrill and small rill erosion only. Due to these limitations, some researchers, e.g., reference [24] have criticized RUSLE for using it in predicting sediment delivery ratios.

In this research, RUSLE is applied to the sub-catchment of the Gomal river to assess erosion rate and test an improved vegetation scenario.

$$A=RKLSCP$$

(1)

where A denotes average annual soil loss (t ha⁻¹ y⁻¹), R is rainfall-runoff erosivity factor (MJ mm ha⁻¹ h⁻¹ yr⁻¹), K is soil erodibility factor (Mg ha h ha⁻¹ MJ⁻¹ mm⁻¹), LS is topographic factor, C is crop management factor (ranges from zero to 1), and P is conservation support practice factor (ranges from zero to 1).

Soil loss predictions using RUSLE yield varying results because the methods used to estimate the factors in Equation (1) are empirically derived. The regional applicability of the RUSLE requires the sub-factors to be adjusted and modified based on the specific characteristics of the study site [25]. Some studies have attempted to evaluate the impact of different RUSLE factors on soil loss. For example, reference [12] found the land slope (LS) as the factor most significant factor affecting soil erosion. The authors in [12] conclude that severe erosion occurs along the drainage channels due to steep bank slopes. Similarly, ref [13] found that the land slope (LS) factor is that with most influence on soil erosion, followed by P, K, C, and R. The support practice P may also be an important factor when the management of soil erosion from the catchment is considered [14].

Assessment of soil loss using RUSLE approach requires a careful selection of methods determining its factors. A recent study by [25] reviews the different factors of USLE and RUSLE, and analyses how various studies around the world have adapted the equations to local conditions.

Rainfall-Runoff Erosivity Factor R

The factor R considers the erosion due to rainfall and runoff, which is the function of long-term rainfall kinetic energy and maximum 30-min intensity [1,26]. For this study, precipitation data were collected from Tank meteorological station (32°12′23″ N, 70°23′30″ E) maintained by the Water and Power Development Authority (WAPDA). The meteorological station does not have an automatic recording rain-gauge, and hence high-intensity rainfall data of considerable length was not available. The rainfall data for the period of 2014–2016 was collected to estimate R using an appropriate empirical equation for the study area. Many regional studies have suggested empirical equations (see in [25]) that can be adopted to estimate R using monthly rainfall data. The advantage of considering these empirical equations is that they do not necessitate long-term and high-resolution rainfall data and are suitable for use in data-sparse catchments like the Gomal river catchment.

Table 1. Different empirical relationships for calculation of rainfall-runoff erosivity factor R.

Serial Number	Equation	Geographic Location	Author(s)
1	$\mathrm{R}=4.79\mathrm{MFI}-142$	Morocco and other locations West Africa	Arnoldus (1980) [28]
2	$\mathrm{R}=4.79\mathrm{MFI}-143$
3	$\mathrm{R}=0.66\mathrm{MFI}-3$
4	$\mathrm{R}=0.5\times \mathrm{p}$	Africa	Roose (1975) [27]
5	$\mathrm{R}=-823.8+5.213\mathrm{p}$	USA	Fernandez et al. (2003) [29] originally developed by the USDA-ARS
6	$\mathrm{R}=0.1215{ \times \mathrm{MFI}}^{2.2421}$	Turkey	İrvem et al. (2007) [30]
7	$\mathrm{R}=839.1{ \times \mathrm{e}}^{0.0008\mathrm{p}}$	India	Nakil (2014) [31]

lists down a few of these empirical equations, which we considered to calculate R but finally adopted the empirical equation by [27] in this research. A comparison of R calculated with the empirical equations (as shown in Table 1) with the R from the global rainfall erosivity dataset is presented in the Results and Discussion section. The global rainfall erosivity dataset is collected from the Joint Research Centre, European Soil Data Centre (ESDAC) (https://esdac.jrc.ec.europa.eu/, accessed on 30 June 2019).

In Table 1, MFI is the modified Fournier index calculated as in Equation (2).

$$\mathrm{MFI}={{\displaystyle \sum }}_{i=1}^{12}\frac{{\mathrm{p}}_{\mathrm{i}}^2}{\mathrm{p}}$$

(2)

where p_i is the monthly rainfall in mm, and p is the annual rainfall in mm.

Soil Erodibility Factor K

The soil erodibility factor K represents the susceptibility of soil to erosion under standard plot conditions. Soil erodibility is scaled from 0–1, depending on the soil texture, e.g., high values are used for silt to fine sand and low values for coarse sand due to its resistance to erosion. To calculate soil erodibility factor K, ref [32] has proposed a relationship shown in Equation (3) that considers the soil textural information, organic matter, information about the soil structure, and permeability. Although Equation (3) is based on the soil data of regional conditions in the United States, in many studies outside the United States researchers have used Equation (3) for soil erodibility factor calculations [25].

$$K=\frac{0.1317\left[2.1\times {10}^{-4}{M}^{1.14}\left(12-a\right)+3.25\left(b-2\right)+2.5\left(c-3\right)\right]}{100}$$

(3)

where M is a parameter based upon soil texture and estimated using Equation (4) as adopted from [33], a is the organic matter content (%), b is a code related to soil structure (1 for very fine granular, 2 for fine granular, 3 for coarse granular, 4 for massive), and c is a code related to soil permeability (1 to 6, fast to very slow draining characteristics).

$$M=\left({\mathrm{m}}_{\mathrm{silt}}+{\mathrm{m}}_{\mathrm{vfs}}\right)\times \left(100-{\mathrm{m}}_{\mathrm{clay}}\right)$$

(4)

where m_clay is a fraction of clay contents with particle size <0.002 mm, m_silt is a fraction of silt contents with particle size ranging from 0.002 to 0.05 mm, and m_vfs is a fraction of very fine sand contents with particle size ranging from 0.05 to 0.1 mm. All fractions are expressed as a percentage.

Soil maps help to define the codes related to soil structure and soil permeability. Soil erodibility can also be identified based on the regional studies, e.g., reference [34] have reported the erodibility factor K values for the Potohar region, Pakistan. In Pakistan, the Soil Survey of Pakistan and the Geological Survey of Pakistan prepare maps that can be used to define erodibility classes. Moreover, some satellite products are also available at a coarser scale for harmonized soil type information. Satellite data of different soil fractions are also available on 250 m grid resolution. However, the erodibility factor for the current study was calculated using the soil texture information extracted from the Food and Agriculture Organization (FAO) soil datasets (www.fao.org, accessed on 23 May 2019). Soil erodibility values were assigned based on the texture and organic matter of the soil. Three soil layers were extracted from the FAO soil database. Soil classes were identified from clay loam to loam with organic matter of more than 2% because the soils in the catchment are very rich in organic carbon. The values assigned for clay loam and loam were 0.26 and 0.28, respectively.

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.

$$LS=\left(m+1\right){(\frac{U}{L_0})}^m(\frac{Sin\beta }{S_0}{)}^n$$

(5)

where U is upslope contributing area per unit width (m²/m), L₀ is the length of the unit plot (m), S₀ 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₀ and S₀ 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,45,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].

$$C=-exp\left[l\left(\mathrm{NDVI}/r-\mathrm{NDVI}\right)\right]$$

(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.

2.3. Sediment Transport Modeling

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 × 80 × 2 m³. 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³ s⁻¹; sediment discharge was considered as 10, 156, and 151 kg s⁻¹ 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.

3. Results and Discussion

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 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 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³s⁻¹, which is about one-quarter of the irrigation demand (indent) of 100 m³s⁻¹. On the other hand, Figure 7b shows that 50% of the time, sediment concentration remained in the zero to 250 ppm range.

A few high sediment concentrations were also observed at high discharges, i.e., >100 m³s⁻¹, 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³s⁻¹. This sediment concentration is likely to be coming from water releases from the GZ dam, which is controlled to 24 m³s⁻¹. 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. 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

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⁻¹ h⁻¹ yr⁻¹ (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].

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).

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.

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).

The results (Figure 9a) show that 46% of the sub-catchment area shows a very severe annual erosion rate (100–500 t h⁻¹) while 16% shows a catastrophic annual erosion rate (>500 t ha⁻¹). Very severe and catastrophic erosion rates have also been reported by different researchers for Pakistan, India, Ethiopia, and Brazil [34,55,56,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)

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. 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. Area of Gomal river sub-catchment in various erosion rate classes.

Erosion Rate (t/ha)	Category 1	RUSLE		RUSLE with Improved C Factor		Change in the Area with Improved C Factor
		Area (km2)	Proportion (%)	Area (km2)	Area (%)	km2	Percent
<2	Very slight	133.3	29.3	135.1	29.6	1.8	1.3
2–5	Slight	0.4	0.1	5.3	1.2	4.9	1247.8
5–10	Moderate	1.0	0.2	9.7	2.1	8.7	848.2
10–50	High	15.4	3.4	111.9	24.6	96.5	625.6
50–100	Severe	21.8	4.8	114.6	25.1	92.8	425.6
100–500	Very severe	209.6	46.1	76.0	16.7	−133.6	−63.7
>500	Catastrophic	73.4	16.1	3.1	0.7	−70.3	−95.8
Total		454.9	100.0	455.6	100.0

¹ Erosion rates categorization is used after Morgan [26].

, 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.

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²) 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³ 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² (~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⁻¹y⁻¹), 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⁻¹y⁻¹), and 17% of the area showed a very severe erosion rate (100–500 t ha⁻¹y⁻¹). 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.