*2.3. Method*

#### 2.3.1. Reservoir Impact Assessment on LR Discharge

The current study intended to estimate the impact of reservoir functioning by using the modest, graphical, and useful method of double-mass curve (DMC) analysis for the consistent and long-term trend examination of hydro-meteorological data. The concept of the double-mass curve is that a plot of the two cumulative quantities during the same period displays a straight line as long as the proportionality between the two remains unchanged, and the slope of the line represents the proportionality. The advantages of this method are that it can smooth a time series and eliminate random components in the series, thus showing the main trends of the time series. In last 30 years, Chinese researchers have explored the effects of soil and water conservation measures and land use/cover changes on runoff and sediment using the DMC method, resulting in some very good outcomes [40]. For the current study, the double-mass analysis of annual LR discharge and LRB precipitation for the long time span of 1956–2016 and the chosen study time period from 2000 to 2016 were done separately to ensure the accuracy and verification of the change points, if any, in the hydrological time series.

The coefficient of variation (CV) was used to determine the variability of climatic and hydrological changes in the LRB by using the long-term available hydro-meteorological time series. CV is defined as:

$$CV = \frac{\sigma}{\mu} \ast 100\% \tag{1}$$

where '*σ*' is the standard deviation and '*μ*' is the mean.

The construction and operation of water conservancy projects (dams, channel modifications, drainage works, etc.) have transformed the seasonal distribution of and caused abrupt changes in streamflow at the basin scales [41–43]. How to attribute the physical causes of hydrological variability and how to correctly identify the human-caused signals from natural hydrological variability are therefore important questions [44–46]. Understanding hydrological variability is initially needed to solve these queries, as well as for hydrological simulation and forecasting, water resource management, control of water disasters, and many other water activities [47]. However, the correct detection and attribution of complex variability in hydrological processes at multi-time scales are still challenging tasks [48,49], and the difficulty has not been resolved, although many methods are presently used [50,51]. There have been a large number of methods such as the moving T-test (Student's T-test) [52], moving F-test (also known Fisher-Snedecor distribution) [53], Mann–Kendall test [54], and Pettitt test [55]. The evaluation of a trend in time series of hydro-meteorological phenomena has been done using the non-parametric Mann–Kendall test (MK) [56–58] in the MS Excel software supported by XLSTAT 2014 macro. This test is extensively used and can deal with missing and distant data. The test has two parameters that are substantial for trend detection: a significance level (*p*) that represents the power of the test and a slope magnitude estimate (MK-S) that represents the direction and volume of the trend. The trends in time series were completed by a calculation of Kendall coefficient 'τ' [59–61]. In the current study, the MK trend at a significance level of 5% (*p* < 0.05) was applied to the hydro-meteorological time series. Specific attention was paid to the streamflow exposure of the LR to the impact of reservoir operations by analyzing the MK trend and changes in it with time.

#### 2.3.2. SWAT Modelling of LR Streamflow

The SWAT model is among the most extensively applied open source, semi-distributed watershed-scale hydrologic models to simulate the water quantity, surface runoff, and quality of streamflow in river channels [62]. According to the working principle of the SWAT model, a watershed is initially divided into sub-basins, and each sub-basin is subdivided into hydrologic response units (HRUs) based on land use, topography, soil, and slope maps. The hydrologic cycle for each HRU is simulated based on the water balance, including precipitation, interception, surface runoff, evapotranspiration, percolation, lateral flow from the soil profile, and return flow from shallow aquifers. In this study, ArcSWAT2012 running on an ArcGIS 10.2 platform was used for watershed delineation and sub-basin discretization, resulting in 21 sub-basins that were further categorized to 149 HRUs. Srinivasan et al. (2010) [63] stated that since the accuracy of simulated streamflow may be reduced by not considering a reservoir or dam in a watershed, the calibration of the reservoir component is needed to improve the accuracy of simulated streamflow. The SWAT model provides four different reservoir outflow estimation methods: measured daily flow, measured monthly flow, average annual release rate for uncontrolled reservoir, and controlled outflow with target release. The selection of the method depends on available data regarding the reservoir. Therefore, the SWAT model was forced to simulate LR streamflow under the reservoir influence by using the default reservoir module of SWAT. The reservoir details including surface area, reservoir water volume, reservoir operational year, and monthly LR discharge data for the time span of 2000–2016 were poured into the model to simulate LR discharge under the hydraulic interventions of the Zhikong and Pangduo dams. The model was calibrated for the years 2005–2010 and validated for 2011–2016 with 500 simulations each using SWAT-CUP (SWAT Calibration and Uncertainty Procedures) algorithm. The global sensitivity method was used to rank the selected sensitive parameters.

The sensitivity and significance degree of each parameter were analyzed by the t-Stat and *p*-value, as well as the p-factor and r-factor; the higher the value of t-Stat, the greater its sensitivity is, and the lower the *p*-value, the greater the sensitivity is. The p-factor is the percentage of data that is enclosed by the 95PPU (95 Percent Prediction Uncertainty) band (ranging from 0 to 1, where 1 shows that all the prediction are within the 95PPU band), while the r-factor is the average width of the 95PPU band divided by the standard deviation of the measured variable (from 0 to <sup>∞</sup>, with 0 showing perfect match). During the calibration and validation periods, the calculated monthly streamflow was compared with the observed data from the Lhasa hydrometric station using the Nash–Sutcliffe coefficient (NSE) [64], the coefficient of determination (R2) [65], and the percent bias (PBIAS, %) [66]. Additionally, the observed and simulated discharge records were statistically tested for the correspondence between them using Pearson correlation [67], Spearman's correlation [67], and Kendall's rank correlation [68] for the credibility verification of SWAT modelling under the reservoir operations for LR streamflow simulation. Additionally, the forecasted river discharges using the observed and simulated hydrological time series were correlated using the same correlation tests.

#### 2.3.3. ARIMA Forecasting of LR Discharge Time Series

An ARIMA time series model, which was pioneered by Box and Jenkins (1970), was employed in this study [28]. In the ARIMA (*p*, *d*, and *q*) model, where, *p* is autoregressive (AR), *d* is differencing, and *q* is moving average (MA). ARIMA models have two common forms: one is non-seasonal ARIMA (*p*, *d*, and *q*) and the other is seasonal ARIMA (*p*, *d*, and *q*) ( *P*, *D*, and *Q*) *S*, where *P*, *D,* and *Q* represent seasonal parts and *p*, *d*, and *q* are non-seasonal parts of the model. SARIMA was applied in the current study using the observation-based LR discharge recorded at the Lhasa hydrometric station and the SWATsimulation-based LR discharge to predict the future hydrological time series for the LRB. The SARIMA was used for predicting the LR discharge using both the observed and simulated hydrological time series in R. The SARIMA models can be used for stationary time series data, which was ensured through decomposition for non-stationarity and log transformation for de-seasonality of the SWAT-simulated and observed hydrological time series. The SARIMA model developed for the observation-based and simulation-based LR hydrological time series was trained for the time span of 2005–2012 and validated for the years 2013–2016 under the reservoir influence to minimize the ambiguity in the prediction of LR streamflow, and LR discharge forecasting was carried out from 2017 to 2025.

We constructed the model and depicted the autocorrelation function (ACF) and partial autocorrelation function (PACF) of model residuals to confirm autoregressive and moving average parameters. The final automatic ARIMA model selection was carried out in the R environment. The ACF and PACF of residuals were determined to evaluate the goodness of fit. The two most commonly used ARIMA model selection criteria, the Akaike's information criterion (AIC) and the Bayesian information criterion (BIC), were examined and compared for ARIMA model selection. The AIC was used for the purpose of selecting an optimal model fit to given data. The model that had the minimum AIC was selected as a parsimonious model [69–72]. The BIC was also utilized for the identification of the best fit model for LR flow prediction. The model with the least BIC was suitable for time series prediction [73].

AIC, in general case, is:

$$AIC = 2k = n1 \text{n} (SSE/\text{n}) \tag{2}$$

where *k* is the number of parameters in the statistical model, *n* is the number of observations, and *SSE* is square sum of error given by

$$SSE = \sum\_{i=1}^{n} \varepsilon\_i^2 \tag{3}$$

BIC, in general, is given by

$$BI\mathbb{C} = k\mathbf{1}\mathbf{n}(n) + n\mathbf{1}\mathbf{n}(SSE/n) \tag{4}$$

R2, root mean square error (RMSE), and mean absolute percentage error (MAPE) were selected to assess the ability of SARIMA in forecasting the LR discharge. A Ljung–Box test [74] at *p* = 0.05 was also performed to make sure best fit of the SARIMA model for both time series for LR discharge. Finally, we applied the best fitted models to forecast the monthly LR discharge using the observation-based and simulation-based hydrological time series. The schematic representation of study design is shown in Figure 2.

**Figure 2.** Schematic representation of reservoir impact evaluation, seasonal AutoRegressive Integrated Moving Average (SARIMA) and Soil and Water Assessment Tool (SWAT) model setup for the current study. LR: Lhasa River; DMC: doublemass curve; MK-S: non-parametric Mann–Kendall test; CV: coefficient of variability; PACF: partial autocorrelation function; AIC: Akaike's information criterion; BIC: Bayesian information criterion; RMSE: root mean square error; MAPE: mean absolute percentage error; Q: Discharge (m3/s).
