Hydraulic Parameters for Sediment Transport and Prediction of Suspended Sediment for Kali Gandaki River Basin, Himalaya, Nepal

Sediment yield is a complex phenomenon of weathering, land sliding, and glacial and fluvial erosion. It is highly dependent on the catchment area, topography, slope of the catchment terrain, rainfall, temperature, and soil characteristics. This study was designed to evaluate the key hydraulic parameters of sediment transport for Kali Gandaki River at Setibeni, Syangja, located about 5 km upstream from a hydropower dam. Key parameters, including the bed shear stress (τb), specific stream power (ω), and flow velocity (v) associated with the maximum boulder size transport, were determined throughout the years, 2003 to 2011, by using a derived lower boundary equation. Clockwise hysteresis loops of the average hysteresis index of +1.59 were developed and an average of 40.904 ± 12.453 Megatons (Mt) suspended sediment have been transported annually from the higher Himalayas to the hydropower reservoir. Artificial neural networks (ANNs) were used to predict the daily suspended sediment rate and annual sediment load as 35.190 ± 7.018 Mt, which was satisfactory compared to the multiple linear regression, nonlinear multiple regression, general power model, and log transform models, including the sediment rating curve. Performance indicators were used to compare these models and satisfactory fittings were observed in ANNs. The root mean square error (RMSE) of 1982 kg s−1, percent bias (PBIAS) of +14.26, RMSE-observations standard deviation ratio (RSR) of 0.55, coefficient of determination (R2) of 0.71, and Nash–Sutcliffe efficiency (NSE) of +0.70 revealed that the ANNs’ model performed satisfactorily among all the proposed models.


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][13][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][21][22][23][24], specific stream power [20,23,24], and flow velocity [20,21,[23][24][25][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.

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 2 area, comprised of elevations of 529~2000 m MSL covering 1317 km 2 (19% coverage); 2000~4700 m MSL covering 3388 km 2 (48% coverage); 4700~5200 m MSL covering 731 km 2 (10% coverage); and elevations greater than 5200 m MSL covering 1624 km 2 (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 −1 in the Tibetan Plateau to about 2000 mm year −1 in the monsoon-dominated Himalayas [8]. 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 •

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)(1998)(1999)(2000)(2001)(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 . Sedimentation lowers the reservoir capacity of the dam annually.

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][29][30]: The mean available power supply over a unit of bed area is calculated by: where w t 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 −1 and is given by: The flow velocity is calculated by Manning's formula: where τ b is the bed shear stress (N·m −2 ), ρ is the density of water (1000 kg·m −3 ), g is the acceleration due to gravity (9.81 m·s −2 ), R is the hydraulic radius (m), i is the slope of the river bed (m·m −1 ), ω is the mean available specific stream power per unit area (w·m −2 ), Q is the observed discharge (m 3 ·s −1 ), v is the flow velocity (m·s −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]:

Development of Different 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].

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, Q s t , depends on two independent variables, the daily average discharge of a river (Q w t ) and the average rainfall (R t ) of a catchment area, and the model is expressed in the form of a regression equation [32,33]: The different linear models were created by considering Q w t , a day lag discharge; Q w t−1 and R t , and a day lag rainfall, R t−1 , were the input variables and the performance of different models was also evaluated.

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]: Different nonlinear models were also created and their performance was evaluated separately.

Sediment Rating Curve
SRC is expressed [34] in the form of: where Q s t is the suspended sediment load (kg·s −1 ), Q w t is the daily average discharge of river, and a and b are coefficients that depend on the characteristics of a river.

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 (Q w t ), a day lag discharge (Q w t−1 ), and average daily rainfall (R t ), a day lag rainfall (R t−1 ), where the output consisted of the average daily suspended sediment load (Q s t ). Out of 2191 data sets, 70% of the data was used for training, 15% for validation, and 15% for testing in the ANNs.

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 (R 2 ), and Nash-Sutcliffe efficiency (NSE) [32,35,36]: The lower the RMSE value, the better the model's performance: where the optimal PBIAS value is 0.0, and positive values indicate a model underestimation bias and negative values indicate a model overestimation bias: 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: The optimal value for R 2 is 1.0; the higher the value of R 2 , the better the model's performance: 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, Q s o,i and Q s p,i are the observed and predicted suspended sediment and Q s o,i and Q s p,i are the average observed and average predicted suspended sediment, respectively.

Results and Discussion
The  (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].

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: 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

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: 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][21][22][23][24][25][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 τ b = 0.   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.030d 1.686 and lower boundary of ω = 0.009d 1.686 for 50 mm ≤ d ≤ 3290 mm; O'Connor's [23] average of ω = 0.002d 1.71 and lower boundary ω = 30 × 1.00865d 0.1d for a particle size of 270 mm ≤ d ≤ 6240 mm; and Williams' [24] lower boundary of ω = 0.079d 1.3 for a particle size of 10 mm ≤ d ≤ 1500 mm ( Figure 6).  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 m 3 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).

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 Q T = Q + kσ n , where Q is the mean discharge, k is the frequency factor, and σ n is the standard deviation of the maximum instantaneous flow, respectively. The frequency factor is given by k = y t −y n s n , where y n is the mean and s n is the standard deviation of Gumbel's reduced variate; y t is given by y t = −ln ln T T−1 . The observed highest flood in 1975 was 3280 m 3 ·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.

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 M w 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. 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].

Hysteresis Curve and Hysteresis Index (HI mid ) Analysis
The relationship between the suspended sediment concentration and fluvial discharge can be studied by the nonlinear relationship between them known as hysteresis [44]. Generally, a clockwise hysteresis loop is formed due to an increasing concentration of sediment that forms more rapidly during rising limb, which suggests a source of sediment close to the monitoring point and sediment depletion in the channel system. Conversely, an anticlockwise hysteresis loop shows a long gap between the discharge and concentration peak, which suggests that the source is located far from the monitoring point or bank collapse [45,46].
Clockwise hysteresis loops were developed, increasing the suspended sediment load on the rising limb of hysteresis from December to July, leading to a maximum value of the suspended sediment load of 10,691 kg·s −1 for a fluvial discharge of 1053 m 3 ·s −1 on August 2009. The suspended sediment load decreased on the falling limb of hysteresis from July/September to November. Overall, these six years were characterized by distinct clockwise hysteresis patterns (Figure 10a).
The HI mid is a numerical indicator of hysteresis, which effectively shows the dynamic response of suspended sediment concentrations to flow changes during storm events [47].
The midpoint discharge was calculated by Lloyd [46] and Lawler [47]: where k is 0.5, Q max is the peak discharge, and Q min is the starting discharge of an event.
The hysteresis index value was calculated by Lloyd [46] and Lawler [47]: where Q sRL and Q sFL are the suspended sediment on the rising and falling limb, respectively.

Yearly Suspended Sediment Yield and Prediction by Different Models
A regression equation derived from the observed data (2006 to 2011) of the suspended sediment versus the discharge of the river shown in Figure 10b is given by: The total suspended sediment yield from the catchment is given by: where Y s is the total annual sediment yield from the catchment, C wi is the suspended sediment concentration in mg·L −1 , Q wi is the fluvial discharge in m 3 ·s −1 , dt is the time interval, t i and t i+1 are the preceding and succeeding time in seconds, respectively. This study showed that the median ASSL transported by KaliGandaki River in the hydropower reservoir was 0.003 Mt during winter, increased to 0.026 Mt during spring, was 41.405 Mt during the summer season, and decreased 0.175 Mt during the autumn season ( Figure 11a). Compared to the seasonal transport of suspended sediment, more than 96% of the suspended sediment was transported during the summer season. This depicts a wide seasonal variability of the suspended sediment caliber, which was nearly 14,000 times higher than the winter season (Figure 11a). The maximum observed ASSL transported by the river was 58.426 Mt in 2009, and after that it decreased (Figure 11b).
The HI mid ≈ 0 indicated a weak hysteresis loop whereas HI mid > 0 indicated a clockwise hysteresis loop, and HI mid < 0 an anticlockwise hysteresis loop. Moreover, the maximum HI mid developed was +2.64 in 2006, depicting the higher sediment transport rate in the rising limb but lower sediment transport rate in the falling limb (Figure 10a). The minimum HI mid developed was +0.53 in 2008, depicting the nearly same paths of the rising and falling limb and indicating a weak hysteresis loop (Figures 10a and 11b).
Different types of MLR, NLMR, general power, log transform linear, and ANNs models, including inputs of the fluvial discharge and average rainfall of the catchment, were developed to select the most suitable model and the results are shown in Tables 2-6, respectively. The performance parameters of MLR and NLMR were satisfactory but predicted negative sediment values for low fluvial discharges and low rainfall, thus these models are unacceptable.    The RMSE, PBIAS, RSR, R 2 , and NSE values of the general power model, log transform models, and ANNs are shown in Tables 4-6. In general, the model simulation can be judged as "satisfactory" if NSE > 0.50, and RSR ≤ 0.70, and if the PBIAS value is ±25% for the stream flow and the PBIAS value is ±55% for the sediment [35]. In this study, the predicted values from ANNs (4−10−1−1) showed an RMSE value of 1982 kg·s −1 , PBIAS value of +14.26, RSR value of 0.55, R 2 value of 0.71, and an NSE value of +0.70, which indicates that the ANNs model's performance was satisfactory. Figure 12a-d show the comparison between the model's predicted transport rates of the suspended sediment discharge in kg·s −1 of the SRC, log transform power model, log transform linear models, and ANNs and the observed suspended sediment values respectively.
Among the SRC, power, log transform, and ANN models, the best median ASSL predicted by the ANN model was 37.611 Mt for the period of 2006 to 2011, whereas the observed median ASSL was 41.678 Mt. The mean ASSL transported by the river to the hydropower reservoir was 40.904 ± 12.453 Mt for 2006 to 2011 and the ANNs' predicted mean value was 35.190 ± 7.018 Mt ( Figure 13). Struck [8] reported that the average annual suspended sediment transported by this river was 36.9 ± 10.6 Mt.  Comparison of different models' predicted and observed yearly total suspended sediment transport. Central lines indicate the median, and bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers, the '+' symbol represents outliers (1.5-fold interquartile range), and the circle shows the mean value.

Conclusions
Shear stress, specific stream power, and flow velocity are important key hydraulic parameters to describe sediment transport in river systems. The monsoon fluvial discharge and landslide dam outburst flood (LDOF) were responsible for boulder movements in Kali Gandaki River, Nepal. The lower boundary equation derived from a broad range of observed and calculated data sets estimated that a maximum particle size of 840 mm was transported by the monsoon fluvial discharge from 2003 to 2011. The ASSL transported by KaliGandaki River in the hydropower reservoir increased from winter to pre-monsoon to monsoon, respectively, and decreased in the post-monsoon period. It was estimated that 40.904 ± 12.453 Mt suspended sediment is lost annually from the higher Himalayas. Additionally, the ANN model provided satisfactory results for the prediction of the suspended sediments' transport rate in Kali Gandaki River, where the annual predicted mean ASSL value was 35.190 ± 7.018 Mt. These parameters are important for visualizing sediment loss from the higher Himalayas to the sea and also for monitoring the dead storage volume of reservoirs for hydroelectric power generation.