**1. Introduction**

It is important to understand the sediment transport and river hydraulics in river systems for a variety of disciplines, such as hydrology, geomorphology, and risk management, including reservoir management. The sediment yield from a catchment is dependent on several parameters, including the topography, terrain slope, rainfall, temperature, and soil type of the catchment area [1]. On the other hand, the yield of sediment fluxes is a combination effect of weathering, land sliding, glacial, and fluvial erosions [2]. Sediment yield from these effects is quite complex [3] and the sediment transport in rivers varies seasonally. The hydrology of Nepal is primarily dominated by the monsoons, characterized by higher precipitation during the summer monsoon from June to September, contributing about 80% of the total annual precipitation [4]. Dahal and Hasegawa [5] reported that about 10% of the total precipitation occurs in a single day and 50% of the total annual precipitation occurs within 10 days of the monsoon period, responsible for triggering landslides and debris flows. The main natural agents for triggering landslides in the Himalayas are the monsoon climate, extremities in precipitation, seismic activities, excess developed internal stress, and undercutting of slopes by streams [6]. The sediments are transported by mountain streams in the form of a suspended load, as well as a bedload [7], depending on the intensity of the rainfall and number of landslide events that occurred within the catchment area [8]. Dams constructed to regulate flood magnitudes limits the downstream transportation of all suspended sediments [9]. However, the annual deposition of sediment in reservoirs decreases the capacity of reservoirs, which compromises the operability and sustainability of dams [10]. Basin morphology and lithological formation governs the amount of sediment crossing a stream station at a certain timepoint, which is generally acted upon by both active and passive forces [11].

Outbursts of glaciers and the failure of moraine dams trigger flash floods [6,12–14], which is one of the main causes of large boulder transportation in high gradient rivers in mountain regions. Different hydraulic parameters, such as shear stress, specific stream power, and flow velocity, can be combined in different ways to form sediment transport predictors [15,16]. Shear stress is a well-known hydraulic parameter that can easily determine the ability of rivers to transport coarse bedload material [17,18]. Similarly, flow competence assessments of floods related to the largest particle size transported are described by the mean flow stress, specific stream power, and mean velocity [19,20]. A number of studies have demonstrated the relationships of shear stress [20–24], specific stream power [20,23,24], and flow velocity [20,21,23–26] of rivers with the size of the boulder movement in the river. It is important to perform this study in Kali Gandaki River as this river originates from the Himalayas and there is limited research on sediment transport by this river, which is crucial in Nepal due to differences in the terrain within a short distance.

In this study, relationships between the fluvial discharge and hydraulic parameters, such as the shear stress, specific stream power, and flow velocity, were generated to derive a lowest boundary equation for the maximum size of particles transported by fluvial discharge in the Kali Gandaki River at a point 5 km upstream of the hydropower dam. The equation was used to calculate the maximum size of particles transported by fluvial discharge during 2006 to 2011. Additionally, it explored the nature of hysteresis loops, developed a hysteresis index, quantified the annual suspended sediment load (ASSL) transport, developed different suspended sediment transport models for Kali Gandaki River, and applied them to predict the suspended sediment rate as well as the average ASSL transport.

#### **2. Materials and Methods**

#### *2.1. Study Site Description*

The Kali Gandaki River is a glacier-fed river originating from the Himalaya region, Nepal [27]. The basin has a complex geomorphology and watershed topography with rapid changes in elevation, ranging from about 529 m MSL to 8143 m MSL. It flows from north to south in the higher Himalayan region before flowing eastward through the lower Himalayan region, entering the Terai plains of Nepal and connecting with Narayani River, which ultimately merges with the Ganges River in India. The snowfall area is separated, with elevation ranges less than 2000 m MSL having no snow cover, 2000 to 4700 m MSL having seasonal snow, 4700 to 5200 m MSL having complete snow cover except for 1 or 2 months, and elevations greater than 5200 m MSL having permanent snow [4]. The Kali Gandaki catchment basin covers a 7060 km<sup>2</sup> area, comprised of elevations of 529~2000 m MSL covering 1317 km<sup>2</sup> (19% coverage); 2000~4700 m MSL covering 3388 km<sup>2</sup> (48% coverage); 4700~5200 m MSL covering 731 km<sup>2</sup> (10% coverage); and elevations greater than 5200 m MSL covering 1624 km2 (23% coverage). Figure 1a shows the different altitude areas' coverage map showing river networks, with the locations of meteorological stations, created in ArcGIS 10.3.1 (ERSI Inc., Berkeley, CA, USA) software. The

elevations of Kali Gandaki River decrease from 5039 m MSL in the higher Himalayas to 529 m MSL at Setibeni, 5 km upstream of the hydropower dam (Figure 1b). This encompasses a wide variation in mean rainfall, ranging from less than 500 mm year−<sup>1</sup> in the Tibetan Plateau to about 2000 mm year−<sup>1</sup> in the monsoon-dominated Himalayas [8].

**Figure 1.** (**a**) Map of Kali Gandaki River catchment area; (**b**) longitudinal profile of Kali Gandaki River.

The main physiographic characteristics of the Kali Gandaki River basin at the hydropower station are shown in Table 1.

The discharge of this river varies seasonally and is dependent on the rainfall received by its tributaries' catchments in addition to the amount of snow melting from the Himalayas. A dam (27◦58 44.88 N, 83◦34 49.68 E) was constructed in 2002 for a hydropower project with a 144 MW power generation capacity, at Mirmi, Syangja.


**Table 1.** Main characteristics of the river basin.

#### *2.2. Data Collection and Acquisition*

The department of Hydrology and Meteorology (DHM), Nepal established a gauge station (28◦00 30 N, 83◦36 10 E) in 1964 (www.dhm.gov.np) and it operated until 1995. The gauge station was not operated during the hydropower dam construction period (1997–2002). The bed level of the dam increased yearly due to the trapping of bedload as well as a suspended sediment load by the dam, which reduced the sediment load downstream. The cross-sectional areas of different years were calculated from area–discharge regression equations obtained from historical discharge rating DHM data (1964–1995). Sedimentation lowers the reservoir capacity of the dam annually.

#### *2.3. Analysis of Shear Stress, Specific Power, and Flow Velocity*

Historical discharge and cross profile elevations data sourced from Nepal Electricity Authority (NEA), Nepal were used to calculate the bed shear stress, specific power developed, and flow velocity by using the following common equations [28–30]:

$$\text{Bed shear stress}, \ \tau\_b = \rho \times \underline{g} \times \mathbb{R} \times \text{i.} \tag{1}$$

The mean available power supply over a unit of bed area is calculated by:

$$
\omega = \frac{\Omega}{w\_t} = \frac{\rho \times g \times Q \times i}{w\_t},
\tag{2}
$$

where *wt* represents the width of the flow, and Ω is the available stream power supply or the time rate of the energy supply to the unit length of the stream in w.m−<sup>1</sup> and is given by:

$$
\Omega = \rho \times \mathbb{g} \times Q \times i. \tag{3}
$$

The flow velocity is calculated by Manning's formula:

$$v = \frac{1}{n} \mathcal{R}^{2/3} \dot{t}^{1/2},\tag{4}$$

where <sup>τ</sup>*<sup>b</sup>* is the bed shear stress (N·m<sup>−</sup>2), <sup>ρ</sup> is the density of water (1000 kg·m<sup>−</sup>3), *g* is the acceleration due to gravity (9.81 m·s<sup>−</sup>2), *R* is the hydraulic radius (m), *i* is the slope of the river bed (m·m<sup>−</sup>1), <sup>ω</sup> is the mean available specific stream power per unit area (w·m<sup>−</sup>2), *<sup>Q</sup>* is the observed discharge (m3·s<sup>−</sup>1), *v* is the flow velocity (m·s<sup>−</sup>1), and *n* is Manning's constant. Manning's constant, *n*, in a steep natural channel is calculated by the equation proposed by Jarrett [31]:

$$n = 0.39 S^{0.38} (3.28R)^{-0.16}.\tag{5}$$

#### *2.4. Development of Di*ff*erent Models for Suspended Sediment Predictions*

The daily suspended sediment load transported by the river in the catchment area is a key indicator to visualize the sediment losses from the higher Himalayas and to assess the reservoir management in hydropower projects. Different researchers have developed multiple linear regression (MLR) and nonlinear multiple regression (NMLR), sediment rating curve (SRC), and artificial neural networks (ANNs) models for the prediction of the daily suspended sediment load [32,33].

#### 2.4.1. Multiple Linear Regression

Multiple linear regression assumes that the sediment load transported by a river is in the linear form. A dependent variable, the suspended sediment load, *Qst* , depends on two independent variables, the daily average discharge of a river (*Qwt*) and the average rainfall (*Rt*) of a catchment area, and the model is expressed in the form of a regression equation [32,33]:

$$Q\_{\mathbb{S}\_l} = \beta\_0 + \beta\_1 Q\_{w\_l} + \beta\_2 R\_l. \tag{6}$$

The different linear models were created by considering *Qwt* , a day lag discharge; *Qwt*−<sup>1</sup> and *Rt*, and a day lag rainfall, *Rt*−1, were the input variables and the performance of different models was also evaluated.

#### 2.4.2. Nonlinear Multiple Regression

The suspended sediment transported by the river shows a dynamic state in a nonlinear form so that it is expressed in the form of a polynomial equation [32,33]:

$$Q\_{s\_l} = \beta\_0 + \beta\_1 Q\_{w\_l} + \beta\_2 R\_l + \beta\_{11} Q\_{w\_l}{}^2 + \beta\_{22} R\_l{}^2 \tag{7}$$

Different nonlinear models were also created and their performance was evaluated separately.

#### 2.4.3. Sediment Rating Curve

SRC is expressed [34] in the form of:

$$Q\_{\mathbb{S}^t} = aQ\_{\mathbb{W}^t}^{\
b} \tag{8}$$

where *Qst* is the suspended sediment load (kg·s<sup>−</sup>1), *Qwt* is the daily average discharge of river, and a and b are coefficients that depend on the characteristics of a river.

#### 2.4.4. Artificial Neural Networks

An artificial neural network is capable of solving complex nonlinear relationships between input and output parameters, which consists of three different layers known as an input, hidden, and output layer, respectively [33]. MATLAB (R2016a) software was used to develop different artificial neural networks, where the input consisted of the average daily river discharge (*Qwt* ), a day lag discharge (*Qwt*−<sup>1</sup> ), and average daily rainfall (*Rt*), a day lag rainfall (*Rt*−1), where the output consisted of the average daily suspended sediment load (*Qst*). Out of 2191 data sets, 70% of the data was used for training, 15% for validation, and 15% for testing in the ANNs.

#### *2.5. Model Performance*

The performance of different models was assessed in terms of the root mean square (RMSE), percent BIAS (PBIAS), RMSE-observations standard deviation ratio (RSR), coefficient of determination (R2), and Nash–Sutcliffe efficiency (NSE) [32,35,36]:

$$RMSE = \sqrt{\frac{\sum\_{i=1}^{n} \left(Q\_{s\_{\bar{v},i}} - Q\_{s\_{\bar{p},i}}\right)^2}{N}} \,. \tag{9}$$

The lower the *RMSE* value, the better the model's performance:

$$PBIAS = \left\{ \frac{\sum\_{i=1}^{n} \left\{ Q\_{s\_{0i}} - Q\_{s\_{pj}} \right\}}{\sum\_{i=1}^{n} Q\_{s\_{0i}}} \right\} \times 100,\tag{10}$$

where the optimal PBIAS value is 0.0, and positive values indicate a model underestimation bias and negative values indicate a model overestimation bias:

$$RSR = \frac{RMSE}{STD EV\_o} = \frac{\left\{\sqrt{\Sigma\_{i=1}^n \left\{Q\_{s\_{0,i}} - Q\_{s\_{p,i}}\right\}^2}\right\}}{\left\{\sqrt{\Sigma\_{i=1}^n \left\{Q\_{s\_{0,i}} - \overline{Q}\_{s\_{p,i}}\right\}^2}\right\}}.\tag{11}$$

The optimal value for RSR is 0.0; the lower the RSR, the lower the RMSE and the better the model's performance.

Coefficient of determination:

$$R^2 = \left\{ \frac{\sum\_{i=1}^n \left\{ Q\_{s\_{o,i}} - \overline{Q}\_{s\_{o,i}} \right\} \left\| Q\_{s\_{p,i}} - \overline{Q}\_{s\_{p,i}} \right\|}{\sqrt{\sum\_{i=1}^n \left\{ Q\_{s\_{o,i}} - \overline{Q}\_{s\_{o,i}} \right\}^2 \sum\_{i=1}^n \left\| Q\_{s\_{p,i}} - \overline{Q}\_{s\_{p,i}} \right\|^2}} \right\}^2. \tag{12}$$

The optimal value for R<sup>2</sup> is 1.0; the higher the value of R2, the better the model's performance:

$$\text{NSE} = \left\{ 1 - \frac{\sum\_{i=1}^{n} \left\{ Q\_{s\_{\upsilon}j} - Q\_{s\_{\upsilon}i} \right\}^2}{\sum\_{i=1}^{n} \left\{ Q\_{s\_{\upsilon}j} - \overline{Q}\_{s\_{\upsilon}i} \right\}^2} \right\}. \tag{13}$$

The optimal value for NSE is 1.0 and values range from −∞ to 1. Values between 0.0 and 1.0 are taken as acceptable levels of performance whereas negative values indicate that the mean observed value is a better predictor than the predicted value, which indicates an unacceptable performance. Here, *Qso*,*<sup>i</sup>* and *Qsp*,*<sup>i</sup>* are the observed and predicted suspended sediment and *Qso*,*<sup>i</sup>* and *Qsp*,*<sup>i</sup>* are the average observed and average predicted suspended sediment, respectively.

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

The average discharge of Kali Gandaki River from 2003 to 2012 was 306 m3·s−<sup>1</sup> with a minimum discharge of 40.73 m3·s−<sup>1</sup> during winter in 2009 and a maximum discharge 2824.50 of m3·s−<sup>1</sup> during summer in 2009. The maximum discharge showed a decreasing trend from 2003 to 2006 whereas an increasing trend from 2007 onwards was observed, as shown in Figure 2a. The yearly transported sediment load increased nearby upstream river bed level elevations of the reservoir (Figure 2b) and sediment deposited into the reservoir decreased the reservoir's capacity. The effects of climate change in the higher Himalayas appeared in the form of uneven patterns of increasing rainfall, glacial rate erosion, and permafrost degradation, resulting in an increase in landslides and debris flows [2], which also reflects the temporal and spatial variation of the water balance components in the Kali Gandaki basin [37]. The amount and intensity of rainfall around its catchment affected the discharge rating curve [27].

**Figure 2.** Yearly (**a**) discharge (NEA 200–2012) (**b**) cross profiles (2002–2011) of Kali Gandaki River.

#### *3.1. Relationship of Shear Stress, Specific Stream Power, and Flow Velocity with Discharge*

The calculated shear stress, specific stream power, and flow velocity of Kali Gandaki River at the discharge gauge station, which was about 5 km upstream from the dam within limited data from 2003 to 2011, was related as:

$$
\pi\_b = 3.143 \times Q^{0.621} \left( R^2 = 0.72 \right) \tag{14}
$$

$$
\omega = 27.40 \times Q^{0.612} \left( R^2 = 0.65 \right) \tag{15}
$$

$$v = 0.108 \times Q^{0.519} \left( R^2 = 0.95 \right). \tag{16}$$

The highest shear stress, specific stream power, and flow velocity were observed during 2008 whereas the lowest were observed during 2007. These parameters were directly related with the hydraulic radius in the case of shear stress and flow velocity, whereas the fluvial discharge in the case of specific power (Equations (1), (2), and (4)). The sedimentation process increased the bed level elevation, changing the cross geomorphology of the bed (Figure 2b). These parameters followed nearly the same trends during the remaining years. The shear stress, specific power, and flow velocity of the river increased the function of the fluvial discharge, as shown in Figure 3a,b and Figure 4a, respectively.

**Figure 3.** Relationship of fluvial discharge and (**a**) shear stress and (**b**) specific power.

#### *3.2. Relationship of Particle Sizes and Fluvial Discharge*

The hydraulic parameters were the shear stress, specific stream power, and flow velocity depict transportation of different particle sizes. When subjected to the same fluvial discharge, the specific power showed an increase of 327 mm to 2062 mm particle size whereas the flow velocity depicted an increase of 37 mm to 1794 mm. The shear stress exhibited an increase of 147 mm to 1492 mm particles, which covered the lowest maximum particle sizes compared to the specific power and flow velocity (Figure 4b). These three parameters were derived from the fluvial discharge, as summarized in the lowest boundary equation form of the fluvial discharge as shown in Figure 4b:

$$d\_{mm} = 0.4 \times Q^{1.093} \text{ (25 mm} \le d \le 840 \text{ mm)}.\tag{17}$$

Equation (17) predicted that from the 2003 to 2011, the discharge during monsoons was capable of transporting an 840 mm particle size. Hydraulic parameters, such as the bed shear stress, specific stream power, and flow velocity, have gained wider acceptability among different researchers [20–26] regarding their useful contribution to the derivation of the relationship between particle sizes and hydraulic parameters. The shear stress and particle size relationship of this study was compared with Costa's [20] average of <sup>τ</sup>*<sup>b</sup>* = 0.163*d*1.213, lower boundary of <sup>τ</sup>*<sup>b</sup>* = 0.056*d*1.213 for 50 mm <sup>≤</sup> *<sup>d</sup>* <sup>≤</sup> 3290 mm; Komar's [21] <sup>τ</sup>*<sup>b</sup>* = 0.164*d*1.21 for 50 mm <sup>≤</sup> *d* <sup>≤</sup> 5000 mm; Lenzi's [22] <sup>τ</sup>*<sup>b</sup>* = 86.629*d*0.25 for 20 mm <sup>≤</sup> *<sup>d</sup>* <sup>≤</sup> 1000 mm; O'Connor's [23] average of <sup>τ</sup>*<sup>b</sup>* = 0.0249 *<sup>d</sup>*1.12 for 270 mm <sup>≤</sup> *d* <sup>≤</sup> 6240 mm; and Williams [24] lower boundary of <sup>τ</sup>*<sup>b</sup>* = 0.17*d*1.0 for 10 mm <sup>≤</sup> *<sup>d</sup>* <sup>≤</sup> 3300 mm (Figure 5).

**Figure 4.** Relationship of the fluvial discharge, (**a**) flow velocity, and (**b**) particle size (boulder).

**Figure 5.** Relationship of shear stress and particle size (boulder) and a comparison with different researchers.

For a comparative study of the specific stream power, the particle size relationship of this river was compared with Costa's [20] average of ω = 0.030*d*1.686 and lower boundary of ω = 0.009*d*1.686 for 50 mm <sup>≤</sup> *<sup>d</sup>* <sup>≤</sup> 3290 mm; O'Connor's [23] average of <sup>ω</sup> = 0.002*d*1.71 and lower boundary <sup>ω</sup> = <sup>30</sup> <sup>×</sup> 1.00865*d*0.1*<sup>d</sup>* for a particle size of 270 mm <sup>≤</sup> *<sup>d</sup>* <sup>≤</sup> 6240 mm; and Williams' [24] lower boundary of <sup>ω</sup> = 0.079*d*1.3 for a particle size of 10 mm ≤ *d* ≤ 1500 mm (Figure 6).

**Figure 6.** Relationship of the specific power and particle size (boulder) and a comparison with different researchers.

The flow velocity and particle size relationship of this study was compared with Costa's [20] average of *<sup>v</sup>* = 0.20*d*0.455 and lower boundary of *<sup>v</sup>* = 0.14*d*0.455 for 50 mm <sup>≤</sup> *d* <sup>≤</sup> 3290; U.S.B.R.'s [20] *<sup>v</sup>* = 0.187*d*0.50 for 1 mm <sup>≤</sup> d <sup>≤</sup> 600 mm; Komar's [21] *<sup>v</sup>* = 0.197*d*0.46 for 8 mm <sup>≤</sup> d <sup>≤</sup> 5000 mm; O'Connor's [23] average of *<sup>v</sup>* = 0.074*d*0.60 for 270 mm <sup>≤</sup> d <sup>≤</sup> 6240 mm; Williams' [24] lower boundary of *<sup>v</sup>* = 0.065*d*0.50 for 10 mm <sup>≤</sup> <sup>d</sup> <sup>≤</sup> 1500 mm; Bradely and Mears' [25] *<sup>v</sup>* = 0.163*d*0.5 for 50 mm <sup>≤</sup> *<sup>d</sup>* <sup>≤</sup> 3290 mm; and Helley's [26] *<sup>v</sup>* = 0.1545*d*0.499 for 1 mm <sup>≤</sup> *<sup>d</sup>* <sup>≤</sup> 600 mm (Figure 7).

**Figure 7.** Relationship of the flow velocity and particle size (boulder) and comparison with different researchers.

The calculated values of the shear stress, specific stream power, and flow velocity were less than the observed values by Fort [6], who reconstructed the 1998 landslide dam located about 76 km upstream of the existing hydropower dam of Kali Gandaki River and estimated the hydraulic parameters with an exceptional dam breach discharge of 10,035 m3 s−1. This high discharge was responsible for the movement of a maximum boulder size of 4300 mm [6]. The higher shear stress, specific stream power, and flow velocity observed due to a higher fluvial discharge after the breaching of landslide dam were responsible for the transportation of larger sized boulders (Figures 5–7).

#### *3.3. Estimation of the Return Period by Gumbel's Distribution*

The flood return period from the historical data of DHM, Nepal can be forecasted by the Gumbel method [38] as *QT* = *Q* + *k*σ*n*, where *Q* is the mean discharge, *k* is the frequency factor, and σ*<sup>n</sup>* is the standard deviation of the maximum instantaneous flow, respectively. The frequency factor is given by *k* = *yt*−*yn sn* , where *yn* is the mean and *sn* is the standard deviation of Gumbel's reduced variate; *yt* is given by *yt* <sup>=</sup> <sup>−</sup>*ln*" *ln*" *<sup>T</sup> T*−1 ##. The observed highest flood in 1975 was 3280 m3·s−1. According to the Gumbel frequency of flood distribution, the highest flood will occur after a 40 year return period, as shown in Figure 8a, and the observed extreme discharge, as shown in Figure 8b.

**Figure 8.** (**a**) Gumbel flood return period, (**b**) extreme fluvial discharge.

#### *3.4. Boulder Movement Mechanisms in the Himalayas*

High gradient river hydraulics are strongly influenced by large boulders, with the diameters on the same scale as the channel depth or even the width [39]. Williams [24] mentioned that five possible mechanisms of boulder transport by high gradient river are by ice, mudflow, water stepwise creep by periodic erosion, undermining of stream banks, and avalanches. The bed forming material remains immobile during typical flows, and larger bed forming particles in steep gradient channels typically become mobile only every 50 to 100 years during a hydrologic event [40]. After that, the gravel stocked in low energy sites during lower floods is mobilized and travels as the bedload [40].

The failure of the mountain slope of Kali Gandaki catchment in 1988, 1989, and 1998 was due to an evolved rock avalanche and caused the damming of the Kali Gandaki River [2]. The shockwaves after the massive 7.8 Mw Gorkha earthquake, Nepal on 25 April 2015 and its aftershocks on 23 May 2015 created cracks in the weathered rocks and weakened the mountain slopes of this catchment, which brought rocks, debris, and mud down into the river [41,42]. The river was blocked about 56 km upstream from the hydropower dam by a landslide on 24 May 2015 for 15 h [41] (Figure 9a,b). The downstream fluvial discharge after the blockage was almost zero and a flash flood occurred after an outburst of the natural landslide dam (Figure 9c,d). Extreme flooding during the monsoon period due to high rainfall and a flash flood (Figure 9b), generated by the overtopping of landslide dams [42], was responsible for the noticeable transport of large boulders in the river bed of Kali Gandaki River.

**Figure 9.** (**a**) Natural landslide dam formation on 24 May 2015 (~56 km upstream of dam), (**b**) lake formation after the blockage of the river, (**c**) downstream fluvial discharge after the blockage of the river, and (**d**) extreme fluvial discharge after breaching of the landslide dam on 25 May 2015 (Source: http://kathmandupost.ekantipur.com/news/2015-05-24/blocked-kali-gandaki-river-flowsagain-with-photos.html).

The combination of fluid stress, localized scouring, and undermining of the stream banks may cause small near vertical displacements of large boulders [43]. Catastrophic events, such as natural dam breaks and debris flows, are responsible for larger translations of boulders in rivers [40,43].

#### *3.5. Quantification and Prediction of the Suspended Sediment*
