Variability and Trend Detection in the Sediment Load of the Upper Indus River

Abstract: Water reservoirs planned or constructed to meet the burgeoning energy and irrigation demands in Pakistan face a significant loss of storage capacity due to heavy sediment load from the upper Indus basin (UIB). Given their importance and the huge investment, assessments of current UIB sediment load and possible future changes are crucial for informed decisions on planning of optimal dams’ operation and ensuring their prolonged lifespan. In this regard, the daily suspended sediment loads (SSLs) and their changes are analyzed for the meltwater-dominated zone up to the Partab Bridge and the whole UIB up to Besham Qila, which is additionally influenced by monsoonal rainfall. The gaps between intermittent suspended sediment concentration (SSC) samples are filled by wavelet neural networks (WA-ANNs) using discharges for each site. The temporal dynamics of SSLs and discharges are analyzed using a suite of three non-parametric trend tests while the slope is identified using Sen’s slope estimator. We found disproportional spatio-temporal trends between SSLs and discharges caused primarily by intra-annual shifts in flows, which can lead to increased trap efficiency in planned reservoirs, especially upstream of Besham Qila. Moreover, a discernible increase in SSLs recorded at Partab Bridge during summer is being deposited downstream in the river channel. This is due to a decrease in river transport capacity in the monsoonal zone. These findings will not only help to identify these morphological problems, but also accurately anticipate the spatio-temporal changes in the sediment budget of the upper Indus River. Our results will help improve reservoir operational rules and sediment management strategies for existing and 30,000-MW planned dams in the UIB.


Introduction
Estimation of the suspended sediment loads (SSLs) is important in the design and operation of water structures and in the planning of sediment management (yield reduction, routing and removal) to preserve their live storage capacities [1][2][3][4][5].The temporal variations and changes in SSLs are also an important indicator of the effectiveness of existing watershed management practices or tectonic and land-sliding activities in the catchment area.Being a water stressed country amongst the top ten most climate-affected countries [6,7], Pakistan has a total water storage capacity of only 30 days (equal to 10% of the annual available water), which has been depleting due to heavy sedimentation transported through the Indus River system from the young Hindukush-Karakoram-Himalaya (HKH) ranges [8].For example, amongst three big reservoirs, the Tarbela dam has lost 35% of its storage capacity since 1974 due to trapping of approximately 8 km 3 of sediment in the reservoir ponding area [9].The Warsak dam constructed on Kabul River has filled with 60 Mt of SSL annually in the 30 years after its construction, and no structural or non-structural remedies can reverse its depleting storage capacity [10].Mangla dam, the second largest Pakistani water storing facility, had an initial storage of 7.1 billion m 3 (BCM), which was reduced to 5.6 BCM in 2005 due to sedimentation.In 2009, an additional 9 m rise of the dam increased the storage to 9.1 BCM, which cost one billion USD over five years.However, the rise created technical problems such as an increase of seepage through the dam embankment in addition to the displacement of 45,000 people living in the vicinity of the dam [11].In view of the transboundary nature of the source of water, such a decrease in water storage capacity in Pakistan exacerbates the instability and geopolitical tensions of the region [12].Hence, the assessment of prevalent sediment patterns and their projected changes are vital for the optimization of sediment management processes to ensure the water and food security in the country and to regulate the transboundary water availability pressures.
Although there are many studies assessing the climate-induced adverse impacts on the UIB river flow patterns [13][14][15][16][17][18][19], few have investigated the impact of flow pattern changes on the sediment load capacity [20,21].Furthermore, the studies conducted in this regard differ widely in their suggested estimates.For instance, the SSL to Tarbela Dam (the country's largest) or at the immediately upstream Besham Qila discharge gauge is reported to range from 200 Mt year −1 -675 Mt year −1 over the past 50 years (Table 1).Such uncertainty leads to poor design quality of the operating rules for existing dams and those under construction.Moreover, the studies have generally estimated the SSL by using empirically-developed sediment rating curves (SRCs), whose accuracy is limited as they oversimplify the relationship between the suspended sediment concentration (SSC) samples and the observed discharges [22][23][24].
Table 1.Estimates published on the suspended sediment load (SSL) of the upper Indus River.−1 )  Estimated by 480
The accuracy of SRCs is also limited since it does not model the complex sediment transport processes related to hysteresis phenomena and marked hydrological variations within the UIB, such as: (a) the fluvial erosion and transport processes, which interact with other sediment-producing processes; (b) temporary sediment storage in the main river channel [34]; (c) aggradation and degradation phases of landslides [35]; (d) on average, 5-10 high flow waves of an average 10-12-day duration during the monsoon period; (e) different transit times of discharge and sediment and their different lag times from several sources to the gauge stations.Given that SRCs are employed in the estimation process, there may be a marked compromise in the design quality of reservoir sedimentation prevention measures, as apparent from the current sedimentation rate of the Tarbela and Mangla dams.Since the assessment of the SSL patterns is important for the management of water-related structures, watershed management practices and the sediment budget of the Indus, it is necessary to detect the temporal changes in sediment transport, which are influenced by the river discharge responses and hysteresis phenomena, requiring frequent discharge and SSC observations.As opposed to the discharge time-series typically available on a daily resolution, the SSCs are intermittently sampled, which can affect the trend outcome needed to reconstruct the non-observed days.However, to deal with the non-linear nature of the time series, the wavelet transform coupled artificial neural networks (WA-ANNs) outperform SRCs, since they are able to model theoretically any kind of relationship between the dependent and independent variables without having to know their physical relationship [36][37][38][39][40][41].The wavelet transform decomposes the data time series up-to J levels in the time, space and frequency domains and reveals the information from a given data scenario [42].The temporal scale of the decomposition provides information on temporary storage, aggradation and degradation phases, high flow waves and transit time of the sediment load in the detail coefficients.Given these details, i.e., the detail coefficients along with the approximation coefficient, ANN precisely models the hysteresis phenomena.WA-ANNs have been used successfully over the last decade for reconstructing the missing data by adjusting the hysteresis phenomena in sediment load processes [43][44][45][46].
In assessing temporal dynamics of SSLs and discharges, non-parametric tests are assumed to be more robust as compared to their parametric counterparts, in view of the fact that the sediment load data are not normally distributed, owing to the highly nonlinear nature of the sediment transport processes.However, several non-parametric tests may also result in distinct estimates, which requires employing a suite of successful non-parametric methods and then quantifying their associated uncertainty to build more confidence in the results.
Analyzing discharges and SSCs at two different sites over the past 50 years, this study for the first time shows how changes in the flow patterns are affecting the sediment transport capacity of the UIB for the meltwater-dominated zone (up to Partab Bridge site) and for the whole UIB (up to Besham Qila), which is additionally influenced by the summer monsoonal rainfall regime.The gaps between intermittently sampled SSCs are filled using the wavelet transforms coupled with the artificial neural networks (WA-ANNs).The temporal discharge and SSL dynamics are robustly assessed using a suite of three widely-used non-parametric approaches, including: (1) the innovative trend analysis (ITA), which can analyze the trends in low, medium and high annual SSLs without requiring any assumptions, such as serial correlation, non-normality, sample numbers and others [47]; (2) the Mann-Kendall (MK) and the seasonal Kendall (SK) tests together with Sen's slope method; the MK test detects a trend in a time series without requiring normally-distributed input data [48,49]; Sen's slope method estimates its true slope, while SK analyzes annual trends by removing the seasonal cycles in a time series; (3) a change point detection test, which reveals the changing tendency in the SSL series on monthly and annual scales [50,51]; (4) mean monthly variations, which detect monthly changes based on differences from the (a) first and last decades and (b) monthly regression equations of the analyzed records.

Study Area and Data Description
With a total length of 2880 km and a drainage area of 912,000 km 2 , the Indus River is one of the largest in south Asia.It starts from China and then travels through India and across the whole of Pakistan, finally draining into the Arabian Sea.The drainage of the Indus River is divided into upper and lower parts, typically at the Besham Qila discharge gauge or around 65 km downstream at, so far, its only reservoir, Tarbela, which is one of the largest earth-filled dams in the world (Figure 1).The Besham Qila site located at an elevation of 580 masl drains the mostly high-altitude area of 165,515 km 2 , 12% of which is covered with the Hindukush-Karakoram-Himalaya (HKH) glaciers and permanent ice, while the seasonal snow cover varies between 3 and 67% [18,52,53].Mean annual discharge of the UIB at Besham Qila is 2405 m 3 /s, which constitutes roughly half of the annual surface water availability in Pakistan [18,53].More than 70-80% of such discharge is generated from the melting of snow and glaciers, making the Indus River amongst the most melt-water-dependent rivers in the world [12].The second study site at Partab Bridge is located at an elevation of 1250 masl about 300 km upstream of Besham Qila, draining around 95% of the cryospheric region and contributing around 75% of the Besham Qila flows.The rest of the Besham Qila flows are mostly received from the monsoonal rainfall from July-September.This 300-km river reach, Partab Bridge to Besham Qila, has gained in importance due to the many planned hydraulic structures.For example, the tenders for three major hydropower projects, Bunji 7100 MW (190 m high), Bhasha dam 4500 MW (272 m) and Dasu 5400 MW (242 m), have been completed for construction located downstream of the Partab Bridge gauge [54].In addition, the river reach contains huge sediment deposits due to landslides and tectonic activities.
Since both gauges feature large drainage areas, overall variations in their discharges and SSCs are not as abrupt as in the small catchments, but such variations are still large (Table 2), indicating the occurrence of frequent hydrological events within the UIB.For instance, 1973, 1988 and 1994 were the exceptional flow years at Besham Qila with a total volume of 98.95, 95.31 and 94.88 billion m 3 (BCM), respectively (Figure 2).The year with the highest peak flow was 1984 with a volume of 89.33 BCM.In contrast, only a 61.54 BCM flow volume was observed in 1982, distinguishing it as an extremely low flow year.For Partab Bridge, the exceptional flow years were 1973, 1994 and 1990 with a total volume of 76.5, 69.7 and 69.6 BCM, respectively.On the other hand, 1965 and 1982, with a volume of 42.16 and 46.8 BCM, respectively, were the extremely low flow years.Based on flow patterns, the UIB can be classified into a low flow cycle of 1974-1977 (dry period) and a high flow cycle of 1988-1992 (wet period) with their annual average volume being 71 and 85 BCM, respectively.In drought years, with wet winter and dry summer, the share of glacier melt increases, maintaining the water supply to the Indus River [12].[2] 21,218 ---Q -100 at Dasu damsite [2] 15,078 ---  The specific suspended sediment load (SSL) from the drainage area of the Indus at Besham Qila is estimated to be 1197 Mt km −2 year −1 , more than 90% of which reaches the Partab Bridge and Besham Qila during the high flow period that spans May-September.Such a heavy sediment load is mainly due to glacial bedrock erosion from a large number of small, but steep catchments that directly discharge into the Indus [55].Generally, the peak SSL correlates well with the peak discharge with a short time lag, particularly for Besham Qila during the monsoon season when discharge varies significantly within a short span of a few days, accompanied by an immediate and large increase in the sediment concentration [56].The SSC average grain size distribution for the UIB is about 45.7% sand, 39.9% silt and 14.4% clay [2].
The daily discharges and the discontinuous suspended sediment concentration (SSC) samples were collected for Partab Bridge over the period 1962-2008 and for Besham Qila over the period 1969-2008.Following the U.S. Geological Survey (USGS) guidelines, discharges at these gauges are measured using AAcurrent meter, while the SSC samples are taken once a week in winter and twice a week in summer, depending on the availability of labor and sampling feasibility [57][58][59].The total SSC samples within the collection periods on record were 3213 and 2117, representing around 22% and 14% of the daily time series for the Besham Qila and the Partab Bridge sites respectively.Due to low sampling frequency at Partab Bridge, we decided to use all available data samples.The long length of these samples improves the learning of the WA-ANN model, which in turn leads to better reconstruction of missing SSLs.The outliers in sediment data samples were excluded by examining the general behavior of the river and river catchment.More details on data collection, data quality and period of records for the Indus River can be found in [33].

Methods
We analyze how changes in the flow patterns are affecting the sediment transport capacity of the UIB specifically for the meltwater-dominated zone (up to the Partab Bridge site) and for the whole UIB (up to the Besham Qila), which is additionally influenced by the summer monsoonal rainfall regime.In order to do this, we analyze the observed discharges and SSCs over the past 50 years.Since the SSCs are intermittently sampled and thus represent only a fraction of the daily discharge series of the study gauges, we reconstructed the SSCs for the non-observed days using the wavelet transforms coupled with the artificial neural networks (WA-ANNs).We then employ three non-parametric statistical approaches to analyze the monthly-to-annual scale temporal dynamics of the reconstructed SSLs and observed discharges.These methods include: (1) innovative trend analysis; (2) Mann-Kendall (and seasonal-Kendall) trend test and Sen's slope method; (3) the Pettitt change point test.We also analyzed temporal dynamics by performing decadal and regressional comparisons.

Wavelet Neural Network
Artificial neural networks are widely used in hydrology and water resources studies for data optimization, reconstruction of missing sediment load and prediction of sediment load trends.The ANN architecture acts as an information processing system containing several non-linear and interconnected elements in the form of layers connected via weights.The multi-layer perceptron (MLP) is a typical ANN, which consists of a number of nodes that are organized according to a particular arrangement.The layers process the information from the input layer to the hidden layer and further the hidden layer to the output layer for the generation of results.Generally, the hidden layers contain non-linear transfer functions to process the non-linear or linear information in order to build a relation between input and output variables.The output layer normally contains a linear transfer function to produce the output outside of the range of −1-1.Moreover, the hidden layers can also vary from single to multiple layers using different numbers of neurons.The size of a hidden layer and neurons within the hidden layer also affect the model performance.Less neurons in the hidden layer may affect the learning process, while more neurons in the hidden layer or the number of hidden layers restrict the efficiency in terms of computational time.The increase of neurons may also cause a network over fitting problem.The work in [60] suggested that the neurons for optimal generalization should be in the range from √ 2n + m to the value 2n + 1, where n and m represent the number of input and output nodes, respectively.
Wavelet transform (WT) decomposes signals into successive wavelet components corresponding to a time-domain signal within a frequency range.The original signal can be represented in terms of a wavelet expansion that utilizes the coefficients of the wavelet functions.Several wavelets can be constructed from a function Ψ known as a "mother wavelet", which is confined in a finite interval.That is, WT decomposes a given signal into frequency bands and then analyses them in time.WT are broadly classified into continuous wavelet transform (CWT) and discrete wavelet transform (DWT).CWT is defined as the sum over all time of the signal to be analyzed multiplied by the scaled and shifted versions of the transforming function Ψ.The CWT of a signal f (t) is defined as follows: where '*' denotes the complex conjugate.CWT searches for correlations between the signal and wavelet function.This calculation is done at different scales of a and locally around the time of b.
The result is a wavelet coefficient W a,b contour map.However, computing the wavelet coefficients at every possible scale (resolution level) necessitates a large amount of data and computation time.
DWT analyzes a given signal with different resolutions for different frequency ranges.This is done by decomposing the signal into coarse approximation and detail coefficients.For this, the scaling and wavelet functions are employed.Choosing the scales a and the positions b based on the powers of two (dyadic scales and positions), DWT for a discrete time series f i becomes: where i is integer time steps (i = 0, 1, 2, ..., N − 1 and N = 2 M ); m and n are integers that control, respectively, the scale and time; W m,n is the wavelet coefficient for the scale factor a = 2 m and the time factor b = 2 m n.The original signal can be reconstructed using the inverse discrete wavelet transform as follows: or in a simple form as: where A M,i is called an approximation sub-signal at level M and D m,i are detail sub-signals at levels m = 1, 2, ..., M. The approximation coefficient A M,i represents the high scale low frequency component of the signal, while the detail coefficients D m,i represent the low scale high frequency component of the signal.
There is a variety of mother wavelets such as the Haar wavelet, Daubechies wavelet, Coiflet wavelet and biorthogonal wavelet.Normally, the Daubechies wavelet, which also belongs to the Haar wavelet, has been performing better in sediment transport processes due to its ability to detect time localization information.Time localization information is useful in handling the seasonality and hysteresis phenomenon in flow discharge and sediment load processes.The Coiflet wavelet is more symmetrical than the Daubechies wavelet.Similarly, the biorthogonal wavelet has the property of a linear phase, which is needed for signal reconstruction [61].The selection of an appropriate mother wavelet depends on the type of application and data characteristics.
Before applying an ANN, the input discharge time series are decomposed using pre-selected wavelets.After data decomposition, a portion of the signal associated with certain frequency bands need to be eliminated if there is a poor correlation between the decomposed signal and the observation data.Only the decomposed signals that have significant correlation with the observation signal are used in the forecast model.Furthermore, on decomposed signals, the permutation of the logsig, tansig, radbas and purelin transfer functions was tested for the hidden and output layers.The Levenberg-Marquardt algorithm was used to train the networks due to its simplicity.The neurons in the hidden layer were selected based on the criteria described by [60].The stopping criteria of the models was a maximum of 1000 epochs.The final networks were saved for later use to reconstruct the missing SSCs in the daily time series.Due to the different data time series at both gauges, we developed two different WA-ANN models.Figure 3 shows the methodology of coupling WT with ANN for forecasting SSC in the study area.
next filter level The performance of the model was assessed employing the correlation coefficient (R), root mean square error (RMSE), mean absolute error (MAE) and the Nash-Sutcliffe efficiency (NSE).The correlation coefficient indicates a perfect fit at 1 and otherwise at 0. Similarly, RMSE and MAE indicate the best model performance when close to 0. The NSE ranges from −∞-1, where 1 represents a perfect match and 0 indicates that the model is no better than simply representing the mean value.The simulated results are normally considered 'good' when the NSE is higher than 0.75 and 'satisfactory' when it lies between 0.36 and 0.75 [62].

Innovative Trend Test
The innovative trend analysis (ITA) [47] divides a time series into two halves, where the latter half is plotted against the first, after being sorted in ascending order.Given both halves are identical to each other, the plot shows a scatter of points along a 1:1 (45°) line on the Cartesian plane.The scatter of points falling above (below) the 1:1 line indicates a monotonically-increasing (decreasing) trend.ITA does not require pre-whitening, a specific sample size, a serial correlation structure of the time series or a normal probability distribution.ITA can easily identify the variations and trends in the lower, medium or higher hydrological processes [63,64]

Mann-Kendall Test
The Mann-Kendall (MK; [65,66]) test can detect a trend in a time series without being affected by the outliers.With the use of normal approximation, the MK test statistic S is calculated as follows: where X i and X j are the adjacent data values, S is the sum of positive or negative signs, n is the number of observations and: The two important parameters of the MK test are the significance level and the slope.The former indicates the strength, while the latter indicates the magnitude and direction of a trend.If there are many tied data values, then the method specified for the number of data values greater than 40 is used ( [66], as reported by [67]).The variance of S (Equation ( 7)) takes into account these ties, where q is the number of tied groups and t p is the number of data in the p group.
After calculating S and its variance, an MK statistic Z is computed using Equation (8).A positive value of Z indicates an upward tend, whereas its negative value indicates a downward trend.If there is no detectable trend, then the MK statistic Z has a standard normal distribution.
To detect the season-wise monotonic trends, a slightly modified version of the MK test, namely seasonal Kendall (SK), is used, which runs the original MK test on each season (k) separately, where k can refer to any period of time within a year (e.g., months or four quarters of a year).The overall S statistic is then computed by adding each SK statistic (S k ) for m number of seasons, and the statistical significance of the trend can be assessed using the outcome of Equations ( 10) and (11) in Equation (8).(10) and: Based on sets of Monte Carlo simulations, [68] show that the presence of a positive serial correlation increases the variance of the distribution of S and thus increases the possibility of rejecting the null hypothesis of no trend; the same was also found by [69].By contrast, negative serial correlation diminishes the variance of the distribution and results in underestimation of the significant trend detection probability.To limit the influence of serial correlation, we applied a correction factor, described by [70] in Equation ( 8), as follows; Normally, the population serial correlation coefficient ρ j is replaced with the sample serial correlation coefficient r j ; 2  where j = 2, 3, ..., m − 1. ( The correction factor η k shrinks (expands) the MK statistics in the presence of positive (negative) serial correlation.
An estimate of trend magnitude, which is closely related to the MK test procedure, is known as Sen's slope estimator [71].The slope estimates of N k pairs of data of the k-season are first computed by: for all 1 ≤ i ≤ j ≤ n k and 1 ≤ l ≤ N k .The median of these P k,l values is Sen's slope estimator P k : Finally, P k is tested by a two-sided test at the (1 − α) × 100% confidence interval, and the true slope can be obtained.More details about the Mann-Kendall and Sen's slope tests can be found in [67,72].

Change Point Detection
We used the Pettitt test [50] to detect the qualitative and quantitative changes in SSL and discharge series.The Pettitt change point test is non-parametric and based on a version of the Man-n-Whitney statistics U j,n as follows: where j = 2, 3, ..., n whereas X i and X j are the adjacent data values, n is the number of observations and sgn can be quantified using Equation (6).The statistics K j and corresponding significance testing are given by: and: If p ≤ 0.05, a significant change point exist.

Decadal Analyses and Linear Regressions
Similar to the innovative trend method of [47], we divided the suspended sediment load (SSL) and discharge data into two time series of one decade each.The first time series consists of the initial decade of the dataset, and the second time series consists of the last decade of the dataset.To determine the mean annual and mean monthly changes, we compared the SSL and discharge shares of pre-selected spatial resolution for both gauges.
At the upper Indus River, the effect of high discharge events is influential; they transport a considerable amount of SSL [55].Therefore, we also explored the mean monthly changes in most effective discharges during the initial and last decades of the datasets.The work in [73] defined the most effective discharge as the midpoint of the range of flows, which over a certain period can transport a considerable proportion of the SSL.The effective discharge can be computed using sediment transport formulae and regional flow duration curves.In the present study, we used the effective discharge (Q/Q avg ) of 2.0-times the average discharge (Q avg ) for Besham Qila and 5.0-times the average discharge (Q avg ) for the Partab Bridge gauge as per the classifications formulated by [56].
To obtain the linear changes in each month during past 50 years, we developed linear regression equations of the reconstructed SSLs and observed discharges.Using these equations, we also quantitatively and qualitatively analyzed the changes.

Results
To analyze the trends in suspended sediment loads (SSLs) of the upper Indus River, we reconstructed the missing SSC data using wavelet neural networks (WA-ANNs) and then estimated corresponding SSLs using measured discharges, i.e., SSC × Q.The reconstructed daily data series in the form of monthly and annual SSLs were used in four different trend analysis techniques, namely: ( 1

Reconstruction of Daily Sediment Load Time Series
Based on several preliminary simulations for both gauges, we eventually trained both networks using 70% of the data for the training, 15% for testing and 15% for validation on a random basis.In a similar way, we also decomposed Q t , Q t−1 , Q t−2 , for Besham Qila and Q t , Q t−1 for Partab Bridge up to seven levels using the Daubechies (db1) wavelet.The best performing WA-ANN architectures reconstructed SSLs with a correlation coefficient R = 0.92 for both sites (Table 3).The RMSE and MAE for Partab Bridge were approximately two times more than Besham Qila; likewise, the standard deviation (SD) and mean in the actual SSC samples (Table 2).This difference shows the complexity in the transport process in the glacier influence zone of the upper Indus River at the Partab Bridge gauge.The NSE, which is used to analyze the model performance, was 0.85 for both stations.Therefore, we consider the WA-ANNs reconstructed suspended sediment load (SSL) series good as the NSE is higher than 0.75 [39,62,74].In addition, both WA-ANN models reconstructed SSLs with a cumulative difference of ±1% with the measurements.Thus, according to another comparison criterion, the models that led to an error between ±10% and ±15% are considered as accurate models [37].A comparison between the mass of suspended sediment sampled daily and computed results using WA-ANN models is also shown in Figure S1 of the Supplementary Material.
The reconstructed results showed a higher mean annual SSL of 171 Mt for Partab Bridge compared to 160 Mt at the downstream Besham Qila site (Figure 4).Moreover, the annual SSLs appear to have been rising at Partab Bridge since 1993 and causing the 10-year moving average to increase.In contrast, the annual SSLs have been decreasing at Besham Qila since 1993 (Figure 4).The similar changes in SSLs are also shown in linear and quadratic trends for both gauges (Figure 5).The statistical parameters of linear and quadratic line fittings are shown in Table S1 (supplementary part).

Innovative Trend Test for Annual Loads
The innovative trend test (ITA) shows a decreasing trend in low annual SSLs at Besham Qila against an increasing trend in high annual SSLs at the Partab Bridge site (Figure 6a,b).The frequencies have been increasing at both gauge sites.On the other hand, the overall annual flows at Partab Bridge show an increasing trend, while there are diverse trends at Besham Qila, where, apart from medium annual flows, the low and high flows have no discernible trend (Figure 6a,b).Contrary to Besham Qila, the increase in flows has also been causing an increase in SSLs at Partab Bridge.However, in the absence of hydraulic structures, urbanization or industrialization along the upper Indus River or within the UIB, this increase in annual SSLs noticed at the Partab Bridge did not appear at the downstream gauge, i.e., Besham Qila (Figure 6).

MK Test for Annual and Monthly Loads
The MK trend analyses show that the annual SSLs at Besham Qila have been decreasing at a rate of 0.634 Mt year −1 (Table 4).Calculating according to the same rate, this indicates a maximum decrease of 34 Mt from the estimate made by [28] (reported by [29]) for the Tarbela dam in 1982 (Table 1).Due to a negative slope of 0.634 Mt year −1 , it is possible that the estimates published in 1970s and 1980s show higher sediment loads compared to our estimate (160 Mt year −1 ) at Besham Qila.In contrast to the results of the MK test, the seasonal Kendall (SK) test showed an annual statistically-significant increasing trend at the Besham Qila (Z = 3.2) and Partab Bridge (Z = 4.1) gauges.This contrast in both tests results arises due to the addition of each season's trend in the SK test (Equations ( 10) and ( 11)).In addition, the results of the SK test are similar either using seasons as four quarters of a year (December-February, March-May, June-August and September-November) or each month as a season.
The monthly SSLs show a significant increasing trend in the winter months (November-February) at Besham Qila with a cumulative magnitude of 0.004 Mt year −1 .This is a slight cumulative magnitude, which is unbalanced by the decreasing trend of −0.24 Mt year −1 alone in August (Table 4).Surprisingly, sandwiching increasing trends, April at Besham Qila shows a declining trend only in SSLs.The monthly SSLs at Partab Bridge, in contrast to Besham Qila, show a declining trend of 0.33 Mt year −1 only in August.This trend is balanced by 0.36 Mt year −1 rise in June and September (Table 4).Despite the diversified trends at both gauges, May showed a statistically increasing and August a statistically decreasing trend, whereas in summer, only August at Besham Qila and June, August and September at Partab Bridge show any trends.However, their contribution (33% and 83%) is still higher than the magnitudes of the trends in the remaining months of the year.In summer (July-September), the mean SSL recorded at Partab Bridge is 141 Mt year −1 , while during the same period or even until October, only 125 Mt passed through the Besham Qila gauge; this apparently indicates a durable deposition of SSLs between both gauges in summer.

Change Point Detection Test
The test results show discernible change points in the monthly SSLs after 1982, whereas no change point was detected in annual SSLs at both gauges (Figure 7).Therefore, it might be possible that the peaks in annual SSLs recorded after 1993 at Partab Bridge gauge in Figure 4 in the absence of an increase in corresponding discharges may have been caused by degradation of landslides, which may have previously blocked the sediments [35].On the other hand, the interventions of landslides are marginal for river flow due to a mean discharge of about 2600 m 3 /s.Thus, the change points in monthly discharges are similar at both gauge stations (Table 5).As the Besham Qila site is located in a monsoon rainfall and snowmelt zone, no change in annual flows indicates a decrease in contribution from these sources (Table 5).
Interestingly, the magnitude of the increasing trend in SSLs over May at Partab Bridge was higher than Besham Qila's SSLs, which makes their mean loads approximately the same after 1997 (Figure 7).After 1997, there was no detectable increase in either parameter at either gauge station.Furthermore, September showed a significant increase of 60% no earlier than 1982, which is the highest magnitude or in the change in SSLs of the analyzed record.Compared to a noticeable increase in SSLs at Partab Bridge, surprisingly, the increasing loads are not being received at the downstream gauge (Figure 7).

Decadal Analyses and Linear Regressions
The decadal analyses only show decreasing trends in SSLs during peak summer (June and July) at Besham Qila and only over August at Partab Bridge (Figures 8ab, S2 and S3 (Supplementary Material)).The directions of changes in monthly SSLs are similar to their corresponding discharges except for July at Partab Bridge and August at Besham Qila.It might be possible that the high SSLs recorded in July at Partab Bridge have been causing the SSLs in the following month of August at the downstream gauge to increase, as shown in Figure 8. Similar deviations can be seen in effective discharges, where the SSL transport capacity of the river has been decreasing in summer (June, August and September) at Besham Qila and only over August at Partab Bridge (Figure 9).
The linear regressions also showed identical directions in the changes of the monthly SSLs and their corresponding discharges, except for April at Besham Qila (Table 6), where SSLs are decreasing against the increase in discharges.Nevertheless, there is a certain sensitive linear correlation between mean monthly SSLs and their corresponding discharges for the months in the effective discharges zone, depicted in Figure 10, where the axes represent the change in mean monthly discharges and SSLs (since 1969 and 1962) determined by linear regression equations (Figures S2 and S3 in the Supplementary Material).As can be seen from Figure 10, the change in SSLs is sensitive to the change in discharges, where a 1% change in discharges, on average, can cause a change of 3% in SSLs in the study area.However, compared to Besham Qila, the transport capacity of the river is more sensitive to the discharge change at Partab Bridge, where, for example, an 11% change in mean monthly discharges caused a 65% change in corresponding SSLs over September (Table 6).This may be due to the location of the major source of eroded sediments in the Karakoram parts of the basin that contributes SSLs disproportionate to its drainage area at Partab Bridge [59].On the other hand, the river slope is mild from Partab Bridge to Besham Qila, which causes substantial sediment storage of the incoming SSLs, particularly in summer.To gain an overall qualitative overview of the trends, we compared the results in Table 7.The comparison reveals that SSLs have been increasing in May and decreasing in August in the study area.Apart from that, they have been monotonically increasing during winter months from November-February and also March and May.Although the annual SSLs at both gauge sites showed minor trends, they are statistically insignificant (Table 7).
Table 7. Qualitative comparison of the trends in SSLs using different methods (blue triangles imply an upward trend, whereas red triangles imply a downward trend; a "-" represents statistically insignificant/no trend).ITA, innovative trend analysis.

Discussion
The WA-ANN models with a decomposition level of 7 (256 days) and a lag time of two and one day for Besham Qila and Partab Bridge, respectively, can precisely find the missing suspended sediment load for a given circumstance of hydrological data of the study area.Our findings show that the variation in flow patterns have been causing a significantly increasing trend in suspended sediment loads (SSLs) in May and a significantly decreasing trend in August at both Besham Qila and Partab Bridge gauges in the upper Indus River (Table 7).Contrary to the increase in high frequencies in low annual SSLs and river flows at Besham Qila (which is additionally influenced by monsoon rainfall), the frequencies in high SSLs and river flows are increasing at the Partab Bridge gauge, which is located just downstream of high elevation glacierized areas of the Karakoram and Himalayas (Figures 1 and 6).Even in the absence of hydraulic structures between both gauges, the high SSLs recorded at Partab Bridge during summer are not being transported to the downstream gauge.Furthermore, the mean monthly linear variations show that an average 1% change in monthly flows can cause a 3% change in SSLs (Figure 10).However, the sediment transport capacity of the river is more sensitive to discharge change from May-August at Besham Qila and in September at Partab Bridge.
The sediment transport processes at the upper Indus River are influenced by hysteresis phenomenon and alternative cycles of dray and wet seasons.Applying simple relationship between water discharge and sediment concentration in the modeling process cannot adjust and model these impact factors.Therefore, a temporal resolution of approximately one year with a lag time of one day in the glacier-influenced zone and two days in the rainfall-influenced zone can reduce the variations in sediment load reconstruction.The reconstruction variations, in particular, increase when for example in conventional methods (sediment rating curves and ANN), temporary sediment storage in the main river channel and different transit times of discharges and sediment from their sources to the gauges are not included.Therefore, the quality of hydraulic design and sediment load trends based on poor sediment load estimation ultimately can affect the accuracy of subsequent studies and the efficiency of the overall hydraulic structure and associated benefits.
Partab Bridge gauge is located just downstream of the snow-fed and glacial melt zone of the upper Indus River.Therefore, the results indicate two types of patterns at Partab Bridge: (1) snowmelt-and (2) glacial melt-dependent SSLs.The former (snowmelt dependent) SSLs have been shifting to the spring months (April, May and June) due to an increase in early snowmelts at low altitudes [15,17,75].Particularly in May, the significant increasing effect of early snowmelt has increased the SSL from 3.3 Mt year −1 -5.6 Mt year −1 , over the last 50 years (Figure 7).The effect of early snowmelt has also been noticed by [53], where they determined a 50 million m 3 increasing rate in May's flow at Partab Bridge.Interestingly, in comparison to Besham Qila (47 million m 3 ), the rate of increase in flow (50 million m 3 ) is higher at Partab Bridge and vice versa in SSLs (Tables 4 and 6).However, this increasing trend in flows extraordinarily increased SSLs at Partab Bridge, where after 1993, SSLs are identical to those of Besham Qila.The identical loads at both gauges point out either no increase in SSLs at Besham Qila's catchment or deposition downstream of Partab Bridge.
On the other hand, retrieval of glacial size depreciates the SSLs in August due to less water availability [14,19,56,76]; the SSLs have decreased to 34% (from 43%) over the past 50 years.It seems that glacial melt has shifted to July and September (Table 4).Although the increasing trend in both months is similar (Table 6), September's flow has remarkably increased the SSL from 9 Mt-15 Mt (similar to regressions where increase is 65%) at Partab Bridge (Figure 7).This significant increasing trend may be caused by the small increase in most effective discharge.It also shows the degree of sensitivity where only an 11% change in discharge caused a 66% change in SSL (Table 6).Furthermore, the remarkable increase in SSLs in September may reduce the reservoirs' life by increasing trap efficiency, where according to existing operation rules, the dams are normally filled to the maximum conservation level as late as 31 of August, such as at Tarbela dam.
Contrary to monotonically-increasing trends in SSL at Partab Bridge (except August), the Besham Qila gauge, located in the snow and rain-fed zone, has diversified mean monthly trends from winter to spring (Table 7).The rise in spring' SSLs at Besham Qila might be due to early snowmelt as at Partab Bridge [53].However, the most surprising trend outcome is the decrease in SSLs during April in contrast to the increase in discharges revealed by regressions (Table 6).In the MK test, April's SSLs also showed a decreasing trend, despite an increasing trends in proceeding and immediately succeeding months (Table 4).In April, half of the flow volume recorded at Besham Qila comes from Partab Bridge [77]; however, corresponding to a 36% linear mean monthly increase at Partab Bridge (Table 6), the increase in flow at Besham Qila is only 7%.Therefore, the corresponding increase in SSLs recorded at Partab Bridge during April may be deposited (due to less SSL transport capacity of the river) between Besham Qila and Partab Bridge, causing high SSLs during the succeeding month (May) at Besham Qila, when the river flows show a significant increasing trend at both gauges.Over August, on the contrary, the declining trend in SSLs at Besham Qila is statistically insignificant and seems to be associated with the decrease in the contribution of SSLs (Table 4) and flow volume (from which 84% of flow comes) from Partab Bridge.
As can be seen in Figure 2, over the past 40 years, at Besham Qila, the average annual volume of water was about 76 billion m 3 (BCMs), while the same average was 56 BCMs at Partab Bridge.That means the catchment at Partab Bridge (denoted by Zone 1) contributes 74.2% of the annual flow volume at Besham Qila.The remaining 25.8% in annual flow volume is contributed from the catchment between Partab Bridge and Besham Qila (denoted by Zone 2).The flow volume in Zone 2 is mostly generated from rainfall and snowmelt [17].The linear trend from 1969-2008 in Figure 5 shows an increase in flow volume at Besham Qila of around 3.90% (denoted by ∆Q), while the same increase at Partab Bridge is around 13.50% (denoted by ∆Q 1 ).The variation of water availability in the area between Partab Bridge and Besham Qila (denoted by ∆Q 2 ) can be approximated using the following mass balance equation: From this equation, we obtain the variation in flow in Zone 2 ∆Q 2 = −38%.As Zone 2 is influenced by rainfall and snowmelt, it seems that the negative variation is attributable to trends of these parameter.These parameters (snowmelt and rainfall) have further been causing a decrease in water availability (between both gauges) required to transport the increased SSLs coming from Partab Bridge.Thus, the annual SSL trends at Besham Qila have shown a decreasing tendency since 1969 (Figure 5).Similar results have also been noted by [12], where the decrease in rainfall in the study area has been buffered by the increase in glacier melts.Additionally, the rise in glacier melt or precipitation over the western region of the upper Indus Basin noted by [16] might have been the cause of the 60% increase in SSLs during September at Partab Bridge.However, this increase has not been received at the downstream gauge, possibly due to a statistically insignificant increase in discharge downstream of the same gauge till Besham Qila (Table 4).
In the future at the upper Indus River, the overall increase in flow volume is expected to reach 7-12% [15].This increase will mostly increase the flow share for the pre and post summer months, which will not be enough (it will be less than the most effective discharges) to transport an additional sediment load.Consequently, the annual SSLs will remain the same or will decrease slightly at Besham Qila.Therefore, in keeping with the current trends, the published sediment load estimates indicate an ongoing decline at Besham Qila (Table 1), since 1970 to the present.Regardless of the increase or decrease in the flow volume, the researchers agree on the shift in flow patterns at the upper Indus River [14,33,78].Since there are neither hydraulic structures at the upper Indus River/basin, nor land use changes that might have affected the situation, in contrast to [79,80] studies for the East or Thames River, the temporal changes in SSLs can only be due to climate change factors.In addition, the statistically-significant monthly SSL trends contradict the previous reservoir sedimentation studies, which simply used the past SSCs without modification to the future predictions, particularly for the hydraulic structures planned upstream of Besham Qila [2,3,29,81].Thus, using modified boundary conditions for reservoir sedimentation studies in the presence of trends can improve the overall quality of hydraulic designs and reservoirs' lives in the study area.
Nevertheless, the variations in SSLs, overall, may have serious implications for water storage, as well as the management of peak supply, peak demand and dam safety, which will require certain changes in the existing reservoirs' operational rules.These changes may include the use of additional (increased) water for power generation during low flows (winter months) and for irrigation or flushing operations in May when more water is available.Flushing in May when crops are at a mature stage and do not require irrigation will also provide the opportunity to re-fill the reservoirs in the succeeding high flow months (June-July).Although the overall flow volume at Besham Qila has been increasing slightly, the flow contribution of the catchment between Partab Bridge and Besham Qila (Zone 2) has been decreasing and causing substantial sediment deposition and an overall decrease in the SSLs received at Besham Qila.Despite the fact that we did not include the impact of landslides on sediment deposition, the current findings are of crucial importance for 143 existing or planned dams and other construction projects in the upper Indus River, especially upstream of the Partab Bridge, which has a glacier-fed catchment and is sensitive to change in river discharges.

Conclusions
Reconstructed suspended sediment load (SSL) time series using wavelet neural networks (WA-ANNs) along with the innovative trend test, the Mann-Kendall test, Sen's slope estimator, the change point detection test and linear regressions have shown a shifting trend from the summer (June, July and August) to the spring and winter months due to a change in water availability at the upper Indus River over the past 50 years.The spatio-temporal trends between discharges and SSLs are disproportionate.This disproportional behavior and the significant trends strongly disconfirm the hypothesis that future inflows and SSLs are similar to the previous ones for reservoir sedimentation studies for the upper Indus River.In addition, the SSLs recorded at Partab Bridge are depositing in the river channel between both gauges.This deposition process has led to a long-term decrease in SSLs, in contrast to a long-term increase in flow volume at the Tarbela dam.For future water and food security along the Indus River command area, it is necessary to estimate the impact of long-term SSL variations on the existing and planned water storage capacities of the reservoirs.Moreover, the impact of planned construction activities along the upper Indus River, which contains enormous sediment deposits, should be evaluated.

Figure 1 .
Figure 1.Locations of study gauges in the study area.Modified from [2].

Figure 3 .
Figure 3. Schematic diagram of a wavelet transform coupled to an artificial neural network (WA-ANN).SSC, suspended sediment concentration.

Figure 6 .
Figure 6.Results of innovative trend test showing a decreasing trend in low and high annual SSLs and flows at the Besham Qila and an increasing trend in high annual SSLs at Partab Bridge, along with an increase in all flows (legends for (a) also apply for (b)).

Figure 7 .
Figure 7. Significant change points in monthly SSLs determined using Pettitt test; black denotes Partab Bridge, and blue denotes Besham Qila.

SSLFigure 8 .Figure 9 .yFigure 10 .
Figure 8. Monthly share of SSL and flow volume in the first and last decade of the analyzed record.
Figure S2: Mean monthly linear and quadratic trends in SSLs and discharges at Besham Qila site from 1969-2008; Figure S3: Mean monthly linear and quadratic trends in SSLs and discharges at Partab Bridge site from 1962-2008.Table S1: Statistical parameters of annual linear and quadratic trends of reconstructed SSLs and observed discharges for the Besham Qila and the Partab Bridge sites.Note: Q s is annual SSL in Mt, Q is annual flow volume in BCM for Besham Qila (1969 ≤ y ≤ 2008) and Partab Bridge (1962 ≤ y ≤ 2008).

Table 2 .
Hydrological and sedimentological characteristics at the Besham Qila and the Partab Bridge gauges.

Table 3 .
Statistics of the best performing WA-ANN architectures for the Besham Qila and the Partab Bridge sites.

Table 4 .
Mann-Kendall's (MK) annual and monthly SSL and discharge trends for the Besham Qila and the Partab Bridge sites.The minus symbol for the MK statistics indicates a downward trend, whereas the (-) symbol without numbers means no trend.

Table 5 .
Significant change points in river flows determined using the Pettitt test.

Table 6 .
Mean monthly linear variations in SSLs and discharges (flows) at both gauges (each month's regression plots are presented in FiguresS2 and S3in the Supplementary Material).