Frequency Analysis of Snowmelt Flood Based on GAMLSS Model in Manas River Basin, China

: With the acceleration of human economic activities and dramatic changes in climate, the validity of the stationarity assumption of ﬂood time series frequency analysis has been questioned. In this study, a framework for ﬂood frequency analysis is developed on the basis of a tool, namely, the Generalized Additive Models for Location, Scale, and Shape (GAMLSS). We introduced this model to construct a non-stationary model with time and climate factor as covariates for the 50-year snowmelt ﬂood time series in the Kenswat Reservoir control basin of the Manas River. The study shows that there are clear non-stationarities in the ﬂood regime, and the characteristic series of snowmelt ﬂood shows an increasing trend with the passing of time. The parameters of the ﬂood distributions are modelled as functions of climate indices (temperature and rainfall). The physical mechanism was incorporated into the study, and the simulation results are similar to the actual ﬂood conditions, which can better describe the dynamic process of snowmelt ﬂood characteristic series. Compared with the design ﬂood results of Kenswat Reservoir approved by the China Renewable Energy Engineering Institute in December 2008, the design value of the GAMLSS non-stationary model considers that the impact of climate factors create a design risk in dry years by underestimating the risk.


Introduction
In recent decades, with global warming and the potential influence of human activities on climate change or directly changing the hydrological cycle [1], the assumption of random independent co-distribution of flood time series has been greatly challenged. In fact, all water-related infrastructures were and are based on conventional stationary frequency analysis methods (assuming that there is no trend or sudden change in flood time series data) for design flood analysis, and the reliability of the results obtained is questioned [2,3]. Therefore, it is particularly necessary to consider non-stationary methods for flood frequency analysis. Several researchers have begun exploring the validity of this hypothesis for flood conditions in many parts of the world, taking into account the effects of global climate variability and human activities [4][5][6][7][8][9]. These studies show that the flood time series clearly violates the assumption of stationarity.
Previous studies have shown that the hydrological series in the Manas River Basin is mainly dominated by abrupt changes [10,11]. The variability shown by hydrological and climatic variables is the main reason for its abrupt change [12,13]. However, with the exception of Chen et al. [10,14], little research has identified the stationarity of flood status in the Manas River basin, which is located in the inland river in the arid region. These results show that the flood frequency analysis based on a stationary condition cannot meet the realistic flood control safety standard.
For flood frequency analysis in non-stationary environments, most hydrologists use indirect methods. The method mainly includes the rainfall-runoff relationship method of the watershed, the decomposition and composition method of time series, and the hydrological model method. In recent years, to restore the original hydrological time series, an increasing number of methods have been adopted to directly study the non-stationary hydrological series by using probability theory and mathematical statistics theory, such as conditional probability method, mixed distribution method and time-variant moment method [15][16][17][18][19]. In the literature, most of the studies on the frequency analysis of nonstationary floods assumed trends in time. However, trends may change in the short-and long-terms as a result of climate variability and the intensification of human activities, which are the true drivers, so it is not entirely correct to use a model that only depends on time to predict. Generalized Additive Models for Location, Scale, and Shape (GAMLSS) [20] is a kind of time-varying moment model, it has a more flexible modeling framework than traditional statistical models. At the same time, it can be added to consider climate change and water conservancy engineering factors to construct a non-stationarity model, so that the model can better describe changes in flood conditions over time [21][22][23][24].
Under the background of environmental changes and the study of the sudden change of the snowmelt flood time series of the Manas River in the study area, the prediction of a model that only depends on time is not completely correct. Therefore, the aim of this paper is to address non-stationary modelling of snowmelt floods in the Kenswat River in the arid area, and demonstrate that the incorporation of climate indices (P 1 , P 3 and T 78 ) may result in appropriate covariates to describe changes of floods. We have used the GAMLSS proposed by Rigby and Stasinopoulos [20] to incorporate external covariates. In addition, in order to analyze the importance of the non-stationary modeling in flood frequency analysis, we compared the results of the Kenswat Reservoir design flood simulated by the non-stationary modeling with the design flood results approved by the China Renewable Energy Engineering Institute in December 2008. The research results can provide scientific guidance for the flood control safety, flood control management and basin planning of the Manas River hydraulic engineering.

Study Area
The Manas River originates in the Yilian Habir Mountain on the northern slope of the Tianshan Mountains in China. The total length of the main stream is 324 km and the basin is bounded by 43 • 27 N to 45 • 21 N latitude and 85 • 01 E to 86 • 32 E longitude. The study area (Figure 1) is the Kenswat hydrological station, with basin area of 4637 km 2 and annual average runoff is 12.21 × 10 8 m 3 . At an altitude of more than 3600 m, it is covered with snow all the year round, and there are modern glaciers with an area of 608.25 km 2 , which is the main source of supply for various rivers. The study area has a typical temperate continental climate with a dry climate. The average annual temperature is 5.9 • C, the annual average rainfall is 338.2 mm and the average annual evaporation is 1550.6 mm [25,26].
The Kenswat reservoir (43 • 58 N, 85 • 57 E) is located 2 km upstream of the Kenswat hydrological station. The normal storage level of the reservoir is 990 m, the maximum dam height is 129.4 m, the total storage capacity is 1.88 × 10 8 m 3 [27,28]. Geological exploration of the Kenswat Reservoir began in January 2003. After years of design planning, the reservoir design flood was approved by the China Renewable Energy Engineering Institute in December 2008, and a preliminary design report of the Kenswat hydraulic engineering was formed in December 2009 [29]. The Manas River is a river with frequent flood disasters. Floods mainly occur in July and August during the flood season, where the average rainfall can reach 43.56 mm and the average daily temperature is 22.4 • C, which fully shows the characteristics of rivers during the flood season in summer. Water 2021, 13, x FOR PEER REVIEW 3 of 20

Dataset
The measured annual maximum peak discharge (Qmax) and daily averaged discharge time series from the Kenswat Hydrometric Station on the Manas River over the period 1957-2006 are used ( Figure 1). From the daily averaged discharge time series, five records are created: annual maximum 1, 3, 7, 15 and 30-day flood volume (Wmax1, Wmax3, Wmax7, Wmax15, and Wmax30), by accumulating the daily discharge value. At the same time, the time series of mean temperature of July and August and the influence of rainfall of 1, 3, 5, 7, 15 and 30 days (P1, P3, P5, P7, P15 and P30) before the occurrence of the flood peak were collected at Kenswat Hydrological Station from 1957 to 2006 as basic data.

Correlation Analysis
In this study, the Pearson correlation coefficient was used to analyze the correlation between the temperature mean series in July and August and the snowmelt flood characteristic time series (satisfying the normal distribution). The Spearman rank correlation coefficient was used to analyze the time series of precipitation affected by the early stage of different days and the snowmelt flood characteristic time series. The statistical significance of Pearson and Spearman tests was set to 5% [30][31][32].

Dataset
The measured annual maximum peak discharge (Q max ) and daily averaged discharge time series from the Kenswat Hydrometric Station on the Manas River over the period 1957-2006 are used ( Figure 1). From the daily averaged discharge time series, five records are created: annual maximum 1, 3, 7, 15 and 30-day flood volume (W max1 , W max3 , W max7 , W max15 , and W max30 ), by accumulating the daily discharge value. At the same time, the time series of mean temperature of July and August and the influence of rainfall of 1, 3, 5, 7, 15 and 30 days (P 1 , P 3 , P 5 , P 7 , P 15 and P 30 ) before the occurrence of the flood peak were collected at Kenswat Hydrological Station from 1957 to 2006 as basic data.

Correlation Analysis
In this study, the Pearson correlation coefficient was used to analyze the correlation between the temperature mean series in July and August and the snowmelt flood characteristic time series (satisfying the normal distribution). The Spearman rank correlation coefficient was used to analyze the time series of precipitation affected by the early stage of different days and the snowmelt flood characteristic time series. The statistical significance of Pearson and Spearman tests was set to 5% [30][31][32].

Generalized Additive Models for Location, Scale, and Shape (GAMLSS) Theory
For non-stationary hydrological time series analysis, the parameters of the selected distribution in the modeling framework can vary as a function of explanatory variables. In this paper, we use the (semi-)parametric regression model proposed by Rigby and Stasinopoulos in 2005, GAMLSS. The model can add a variety of explanatory variables, and fit the linear and non-linear functional relationship between the statistical parameters of the response variable series and the explanatory variables. In this study, the model response random variable Y (Q max , W max1 , W max3 , W max7 , W max15 and W max30 in this paper) has a parametric cumulative distribution function, and the parameters can be modeled as a function of selected covariates, namely time or climate index. Therefore, the stationary model (model 0) was established, in which the distribution parameters do not depend on covariates; the time-varying model (model 1), in which the distribution parameters vary as function of time only; and the model that incorporates covariates (model 2), in which the distribution parameters can vary as a function of climate.
The GAMLSS model assumes that the independent random variable observations y i , i = 1, 2, . . . , n obey the probability density function f y parameter vector corresponding to a certain moment, n is the number of observations, and p is the number of distribution (statistic) parameters. The parameters µ i and σ i are generally defined as the location parameter and the scale parameter vector, corresponding to mean vector and mean squared (or coefficient of variation) vector, respectively, represented as random variables. The other parameters in the distribution are collectively referred to as shape parameters. The general shape parameters are only two at most, and skewness vector and kurtosis vector of the random variable series are represented by υ i and τ i , respectively. Let g k (•) denote θ k a monotonic function relationship with the corresponding explanatory variable X k and the random effect term, generally expressed as: where θ k are vectors of length n; β k = β 1k , β 2k , · · · , β I k k T , which is a regression parameter of length I k Vector; X k is an explanatory variable matrix of n × I k ; Z jk is a known fixed design matrix of n × q jk ; γ jk is a normal distribution random variable vector of q jk dimension; Z jk γ jk is the j-term random effect term; q jk is the random influence factor dimension of the j-term random effects. The distribution characteristics of hydrological time series are often described by statistical parameter mean, mean square error and skewness coefficient. The model usually adopts a two-parameter or three-parameter probability distribution function. In this paper, a two-parameter model with linear variation of location parameter θ 1 and scale parameter θ 2 with covariate is used.
In this paper, we considered as candidates four widely used two-parameter distribution functions in modelling streamflow data (Table 1): Log normal (LOGNO), Gamma (GA), Gumbel (GU), Poisson inverse Gaussian (PIG) [7]. It can be seen from these alternative distribution functions that they have the characteristics of both exponential distribution and power function, which accord with the actual situation of hydrological time series. When the model incorporates the selected covariates, the relation between the distribution parameters and the selected covariates will increase the complexity of the model. In order to avoid model overfitting, the Akaike information criterion (AIC), the global fit deviation (GD) and the Schwarz Bayesian criterion (SBC) [20,33] were used to optimize the model selection. Table 1. Summary of the four two-parameter distributions considered in this study to model the streamflow data. LOGNO means Log normal distribution; GA means Gamma distribution; GU means Gumbel distribution; PIG means Poisson inverse Gaussian distribution.

Distribution Probability Density Function Distribution Moments
The GD of GAMLSS model is defined as follows: where l θ i is the logarithmic likelihood function corresponding to the estimated value of the regression parameter. Meanwhile, the generalized Akaike information criterion (GAIC) is introduced for judgment, and its expression is: where d f is the overall degree of freedom of the model, # is the penalty factor. If the penalty factor # = 2, it is called the AIC; if # = log(n) (n is the sample size of the explanatory variable), it is called the SBC. The AIC and the SBC are two special cases of the GAIC. The model with the smallest GAIC value is taken as the optimal model. Because the value of the maximum likelihood does not provide information about the quality of the fitting [34]. Therefore, we examine the first four statistical moments of the residuals and the Filliben correlation coefficients [35] ensure that the selected models can adequately describe the systematic part. For the remaining independent and identically distributed random noise, we use the residual diagnostic plot (residuals vs. response, qq-plots and worm plots) for visual inspection [36,37]. Table 2 summarizes the Pearson correlation coefficients of the mean temperature series in July and August (T 78 ) and the snowmelt flood characteristic time series. It can be seen from the table that the Pearson correlation coefficient between the mean temperature series and the annual maximum peak discharge (Q max ) and annual maximum flood volume series (W max ) is between 0.3-0.5 (statistically there is a medium correlation). For all correlation coefficients the null hypothesis that the correlation coefficient is equal to 0 can be rejected at a 0.05 significance level. Not all of the time series non-stationarity of snowmelt flood is caused by temperature changes, there are also by rainfall, underlying surface and other factors. Therefore, there is a moderate correlation between the mean temperature series and the characteristic series of snowmelt flood, and the calculated results are reasonable and reliable. Table 2. Summary for the correlation between the mean temperature series in July and August and the snowmelt flood characteristic series: r represents the Pearson correlation coefficient; p-value means the significance value; Q max means the measured annual maximum peak discharge; W max1 , W max3 , W max7 , W max15 , and W max30 mean the annual maximum 1, 3, 7, 15 and 30-day flood volume, respectively.  Figure 2 shows the Spearman rank correlation coefficients of rainfall series of different days before flood and snowmelt flood characteristic time series. It can be found from Figure 2 that, except for the P 15 and P 30 , there is a strong correlation between the previousperiod rainfall series of other days and the snowmelt flood characteristic series. The prior-period rainfall series with the largest correlation coefficient is selected as the optimal correlation series. The Q max has the strongest correlation with the P 3 , while the W max1 , W max3 , W max7 , W max15 and W max30 have the best correlation with the P 1 , and all the correlations can be tested by a 0.05 significance level. The correlation coefficients between the P 1 and the W max1 and W max3 are all between 0.5-1.0, showing a strong correlation; while the correlation coefficients between the P 1 and the W max7 , W max15 and W max30 , the P 3 and the Q max are all between 0.3-0.5, showing a medium correlation. Meanwhile, the Spearman rank correlation coefficients of the W max1 and W max3 , W max7 , W max15 and W max30 are calculated respectively. The study found that the correlation is gradually decreased but not significant, and W max1 and W max30 is the lowest value (0.879), but there is still a strong correlation. For the rainfall series of different days before a flood, the correlation coefficients between P 1 and P 3 , P 7 , P 15 and P 30 show a downward trend, among which the correlation coefficient between P 1 and P 30 is only 0.255, and the correlation coefficient between P 3 and P 30 is 0.610. Therefore, the P 1 and P 3 with the strongest correlation with the snowmelt flood characteristic series are selected as the influencing factors of the model. tively. The study found that the correlation is gradually decreased but not signif Wmax1 and Wmax30 is the lowest value (0.879), but there is still a strong correlatio rainfall series of different days before a flood, the correlation coefficients betwe P3, P7, P15 and P30 show a downward trend, among which the correlation coeff tween P1 and P30 is only 0.255, and the correlation coefficient between P3 and P3 Therefore, the P1 and P3 with the strongest correlation with the snowmelt flood c istic series are selected as the influencing factors of the model.

Results with Stationary Approaches: Models 0
This section presents the fitted stationary model (models 0) for the Kenswat hydrological station control basin of the Manas River. Table 3 summarizes the optimal fitting distribution of the characteristic series of snowmelt flood under the stability model using AIC. It can be seen from the table that for the Q max the LOGNO distribution and the PIG distribution have similar fitting effects; for the W max , the LOGNO distribution and the GA distribution have similar fitting effects. In the four candidate distributions models, the AIC value of the LOGNO distribution is the smallest, which is the optimal fitting distribution.  Table 4 and Figure 3 summarize the fitting quality of the optimal fitting distribution in model 0, which are based on the residual plots and the estimation of the first four moments of the residuals. When the snowmelt flood sample size of 50, the Filliben's coefficient is greater than or equal to 0.977 to pass the 5% significance level. The results show that the residual Filliben's coefficients of the optimal fitting distribution model for the snowmelt flood characteristic series are all less than 0.977, which fails the significance test. Moreover, as can be seen from Figure 3, part of the standard residual points of the optimal model of the annual maximum flood time series are outside the 95% confidence interval, which does not meet the requirements. It can also be seen from the residual distribution moments of    Based on the above analysis, model 0 has low fitting accuracy and poor fitting effect. Therefore, the use of a stationary model under the two-parameter distributions can no longer be satisfied, and it is necessary to consider the frequency analysis and calculation of the non-stationary model for the snowmelt floods characteristic series under the mutation condition. The non-stationary model is constructed with time t as the explanatory variable, and the cumulative probability distribution parameters θ 1 and θ 2 , namely the mean and variance of the corresponding characteristic series of snowmelt flood are considered. We use GD, AIC and SBC discriminant criteria to determine the optimal fitting distribution, the optimal covariate of the distribution parameters, and the functional relationship between the distribution parameters and the optimal covariates under the non-stationary model, as shown in Table 5. For the Q max , the LOGNO distribution is the best fitting distribution. The distribution parameter θ 1 exhibits a linear dependence on time t in the four candidate distributions models. The distribution parameter θ 2 exhibits a linear dependence on time t in the LOGNO, GA and GU distribution models, but is independent of the time trends in the PIG distribution model. Table 5. Summary for the fitted models type 1 and the type of dependence between time and the distribution parameters: t means linear dependence; ct refers to a parameter that is constant. Global fit deviation (GD) denotes the evaluation value of the GD. AIC denotes the evaluation value of the AIC criterion. Schwarz Bayesian criterion (SBC) denotes the evaluation value of the SBC criterion. Lower AIC and SBC represent a better performance. For the W max , the LOGNO distribution and the GA distribution have similar fitting effects. However, among the four parameter distributions, the LOGNO distribution has the smallest AIC value, which is the best fitting distribution. In the LOGNO and GA distribution models, the distribution parameters exhibit a linear dependence on time t. In the GU and PIG distribution models, the distribution parameter θ 1 exhibits a linear dependence on time t, and the distribution parameter θ 2 is constant and independent of the time trends. Figure 4 and Table 6 summarize the fitting quality of the optimal distribution of Model 1, which are based on the residual plots and the estimates of the first four moments of the residuals. The results show that the Filliben's coefficient of the fitted residuals are all greater than or equal to 0.977, and all pass the significance test. At the same time, it can be known from the residual distribution moments of each model that the residuals of each model obey the normal distribution well. Therefore, the models fit the data adequately.

Analysis of Optimal Model Fitting Results
It can be seen from Figure 5 that most of the measured points of the snowmelt flood characteristic series are within the interval of 5th to 95th percentile gray scale, indicating that the LOGNO distribution parametric model can capture well the variability exhibited by the data. The GAMLSS non-stationarity model with time as the covariate can be used for trend analysis of snowmelt flood feature series. The characteristic series of snowmelt flood increased with time, and the fitting effect is close. The larger the quantile, the more obvious the trend of the increase in the quantile curve. Among them, the upward trend of the quantiles after 1993 have increased significantly. Taking the 95th percentile as an example, the specific performance was that the growth rate of the Qmax after 1993 increased by 39.03% compared with that before 1993, and the average growth rate of the Wmax after 1993 increased by 30.85% compared with that before 1993. However, with the passage of time, the variation trend of the snowmelt flood characteristic series is not an infinite increase. In actual conditions, the variation of snowmelt flood is affected by many factors such as climate change and human activities. The characteristic series of snowmelt flood should fluctuate. Therefore, model 1 can only describe the trend of the series over time, so it cannot fully describe the dynamic change process of the characteristic series under the influence of many factors.

Analysis of Optimal Model Fitting Results
It can be seen from Figure 5 that most of the measured points of the snowmelt flood characteristic series are within the interval of 5th to 95th percentile gray scale, indicating that the LOGNO distribution parametric model can capture well the variability exhibited by the data. The GAMLSS non-stationarity model with time as the covariate can be used for trend analysis of snowmelt flood feature series. The characteristic series of snowmelt flood increased with time, and the fitting effect is close. The larger the quantile, the more obvious the trend of the increase in the quantile curve. Among them, the upward trend of the quantiles after 1993 have increased significantly. Taking the 95th percentile as an example, the specific performance was that the growth rate of the Q max after 1993 increased by 39.03% compared with that before 1993, and the average growth rate of the W max after 1993 increased by 30.85% compared with that before 1993. However, with the passage of time, the variation trend of the snowmelt flood characteristic series is not an infinite increase. In actual conditions, the variation of snowmelt flood is affected by many factors such as climate change and human activities. The characteristic series of snowmelt flood should fluctuate. Therefore, model 1 can only describe the trend of the series over time, so it cannot fully describe the dynamic change process of the characteristic series under the influence of many factors. by 39.03% compared with that before 1993, and the average growth rate of the Wmax after 1993 increased by 30.85% compared with that before 1993. However, with the passage of time, the variation trend of the snowmelt flood characteristic series is not an infinite increase. In actual conditions, the variation of snowmelt flood is affected by many factors such as climate change and human activities. The characteristic series of snowmelt flood should fluctuate. Therefore, model 1 can only describe the trend of the series over time, so it cannot fully describe the dynamic change process of the characteristic series under the influence of many factors.  .The red points indicate the actual measured points of the snowmelt flood feature series; the light gray region represents the area between 5th and 25th quantile percentile; the dark gray region in the middle represents the area between 25th and 75th percentile; the upper light gray region represents the area between the 75th and 95th percentile; the middle black curve represents the median (50th percentile).

Based on the Results of the Non-Stationary Model with Climatic Factors: Model 2 4.5.1. Model Fitting Evaluation
On the basis of Model 1, the climate factors (P 1 , P 3 and T 78 ) are used to replace the distribution parameter explanatory variable time t as a new interpretation, and the GAMLSS non-stationary model with climate factors as covariates (model 2). Table 7 summarizes the selected distributions as well as the optimal explanatory variable of the distribution parameters and the functional relationship between the distribution parameters and the optimal explanatory variable for Model 2. Table 7. Summary for the fitted models type 2 and the type of dependence between Climate factor and the distribution parameters: ct refers to a parameter that is constant. GD denotes the evaluation value of the GD. AIC denotes the evaluation value of the AIC criterion. SBC denotes the evaluation value of the SBC criterion. Lower AIC and SBC represent a better performance. For the Q max , the LOGNO distribution is the best fitting distribution. Among the three candidate indicators of T 78 , P 1 , and P 3 , P 1 does not pass the screening, indicating that the rainfall of the day before the occurrence of the flood peak is not suitable for describing the non-stationary changes of the Q max . The T 78 and P 3 indicators tend to describe the linear variation of the distribution parameter θ 1 , that is, the mean (location parameter) change of the Q max , which are more susceptible to the influence of the temperature in July and August and the impact of the rainfall in the 3 days before the flood peak appears. Among them, the P 3 indicator has an influence on the Q max of all distribution models, which is only reflected in the linear dependence between the Q max and the P 3 . The T 78 indicator also affects the Q max of models other than the PIG distribution, which is also only reflected in the linear dependence between the Q max and the T 78 . The distribution parameter of all distribution models is constant, indicating that the variance of the Q max has little to do with climatic factors.

Distribution
It can be seen in table that the LOGNO and GA distributions offer the best overall results in modelling the annual maximum flood time series frequency. The GA is used to better fit the W max1 and W max3 ; LOGNO distribution is used to better fit the W max7 , W max15 and W max30 . Among the three candidate indicators of T 78 , P 1 , and P 3 , P 3 does not pass the screening, indicating that the rainfall of the day before the occurrence of the flood peak is not suitable for describing the non-stationary changes of the W max .
The T 78 indicator of the W max1 , W max7 and W max15 are mainly expressed as the linear dependence of the distribution parameter θ 1 , while the P 1 indicator is mainly expressed as the linear dependence of the distribution parameters θ 1 and θ 2 . Therefore, the mean value is mainly affected by the temperature in July and August and the rainfall of the day before the flood peak, while the variance is mainly affected by the rainfall of the day before the flood peak. The T 78 indicator of the W max3 and W max30 are all expressed as the linear dependence between the distribution parameters θ 1 and θ 2 , and the P 1 indicator is also expressed as the linear dependence between the distribution parameters θ 1 and θ 2. Therefore, the mean value is affected by the temperature in July and August and the rainfall of the day before the flood peak. However, the variance of the W max3 is mainly affected by the rainfall of the day before the peak, and the temperature is less affected; while the variance of the W max30 is mainly affected by the temperature in July and August, and the influence of the rainfall of the day before the peak is small.
The results of Figure 6 and Table 8 are similar to the previous analysis results in Section 4.4.1 (Figure 4 and Table 6). This result supports the inference that the models fit the data adequately.   13 °C compared with the average temperature in 1956-1995; the contribution rate of precipitation increase to runoff increase is 59.64% [11,12]. The variation of the underlying surface also has an effect on the series characteristics, but the effect is not significant, because the underlying surface in the study area is less affected by human activities.    Figure 7 summarizes the quantile gray-scale map of the optimal fitting distribution model with the climate factor as covariable for the snowmelt flood characteristic series. This parametric model (model 2) can capture the dynamic change process of the snowmelt flood characteristic series well under the influence of environmental change, which is the most obvious process in the description of the large flood events in 1996 and 1999. In the case of ignoring individual outliers in the figure, the main cause of the dynamic change in the snowmelt flood characteristic series in the figure is climate change, which mainly includes change in temperature and rainfall. According to statistics, the temperature and rainfall in the Manas River basin have shown an increasing trend in recent years. The average temperature in the basin during 1996-2014 has increased by 2.13 • C compared with the average temperature in 1956-1995; the contribution rate of precipitation increase to runoff increase is 59.64% [11,12]. The variation of the underlying surface also has an effect on the series characteristics, but the effect is not significant, because the underlying surface in the study area is less affected by human activities. Figure 8 and Table 9 summarize the 90th, 95th and 98th percentile curves, corresponding extreme values and occurrence years of the Q max and W max under the optimal fitting distribution model with climate factor as the covariate.
Model 2 indicate the existence of periods in which flood frequency experienced significant variability (decreases and increases). The maximum values of the 90th, 95th and 98th percentiles curves of the snowmelt flood characteristic series all appeared in 1996, while the minimum values all occurred in 1972. Due to the impact of climate change, the dynamic range of the Q max at the 98th, 95th and 90th percentiles are 351~1459 m 3 s −1 , 341~1172 m 3 s −1 , 328~966 m 3 s −1 , respectively. The dynamic change process of the W max1 , W max3 , W max7 , W max15 and W max30 is basically similar to that of the Q max , the dynamic range are shown in Table 9. In the non-stationary snowmelt flood frequency analysis, the 98%, 95%, and 90% quantiles represent the flood events with the probability of exceeding 0.02, 0.05, 0.1 (i.e., return period of 50, 20 and 10-year) respectively. Table 9 also shows the design flood results of Kenswat Reservoir approved by the China Renewable Energy Engineering Institute in December 2008. Comparing the design flood results of Kenswat Reservoir with the snowmelt flood quantile values, it can be seen that the snowmelt flood value should show a dynamic change process under the combined influence of climate change and human activities, that is, it should have a dynamic range of change. The design flood value is a static behavior which is used to measure the snowmelt flood value under unstable conditions and can lead to two possible major problems: In dry years, it may appear conservative, while in wet years, especially in years when major floods occur, there may be certain risks. the snowmelt flood characteristic series in the figure is climate change, which mainly in-cludes change in temperature and rainfall. According to statistics, the temperature and rainfall in the Manas River basin have shown an increasing trend in recent years. The average temperature in the basin during 1996-2014 has increased by 2.13 °C compared with the average temperature in 1956-1995; the contribution rate of precipitation increase to runoff increase is 59.64% [11,12]. The variation of the underlying surface also has an effect on the series characteristics, but the effect is not significant, because the underlying surface in the study area is less affected by human activities.
Water 2021, 13, x FOR PEER REVIEW 16 of 20 (e) Wmax15 ( f) Wmax30 Figure 7. The series of snowmelt flood is modeled using the optimal distribution, with parameters depending on climate factor (Model 2).The red points indicate the actual measured points of the snowmelt flood feature series; the light gray region represents the area between 5th and 25th quantile percentile; the dark gray region in the middle represents the area between 25th and 75th percentile; the upper light gray region represents the area between the 75th and 95th percentile; the middle black curve represents the median (50th percentile). Figure 8 and Table 9 summarize the 90th, 95th and 98th percentile curves, corresponding extreme values and occurrence years of the Qmax and Wmax under the optimal fitting distribution model with climate factor as the covariate. The series of snowmelt flood is modeled using the optimal distribution, with parameters depending on climate factor (Model 2).The red points indicate the actual measured points of the snowmelt flood feature series; the light gray region represents the area between 5th and 25th quantile percentile; the dark gray region in the middle represents the area between 25th and 75th percentile; the upper light gray region represents the area between the 75th and 95th percentile; the middle black curve represents the median (50th percentile). (e) Wmax15 ( f) Wmax30 Figure 8. The 90%, 95% and 98% quantile plots of the optimal fitting distributions for Model 1 and Model 2. Figure 8. The 90%, 95% and 98% quantile plots of the optimal fitting distributions for Model 1 and Model 2.

Conclusions
The flood frequency under non-stationary condition between the annual maximum peak discharge and the annual maximum flood volume series in Kenswat Reservoir of Manas River covering the period 1957-2006 is analyzed. Use GAMLSS theory to construct a traditional stationarity model (model 0), and two non-stationarity models based on time as a covariate (model 1) and based on climate factors as a covariate (model 2). The main findings of this work can be summarized as follows: Departures from the traditional assumption of stationarity in the snowmelt flood series in the Manas River are clear.
In the modelling of time-varying parameters, GAMLSS provides a flexible modeling framework to represent the non-stationarity in snowmelt flood distribution. The study found that the characteristic series of snowmelt flood showed an increasing trend over time in the Kenswat Reservoir control basin.
The climate change-related covariables were incorporated into GAMLSS framework to model snowmelt flood, where the location parameter of the annual maximum flood peak series depended on the T 78 and P 3 indicators. The location parameter of the annual maximum flood volume series depended on the T 78 and P 1 indicators. Moreover, the scale parameter can be related only to the P 1 indicator, and does not show any significant dependence on temperature. The covariate model that incorporates the effects of rainfall and temperature can better describe non-stationarities in the frequency and magnitude of the snowmelt flood in Kenswat Reservoir.
Comparing the results obtained in Model 2 with the Kenswat Reservoir design flood results approved by the China Renewable Energy Engineering Institute in December 2008, the snowmelt flood value has a dynamic range, while the design flood value is a static behavior. For snowmelt flood time series with 0.02, 0.05, and 0.1 annual exceeding probabilities (corresponding to 50, 20, and 10-year regression periods under stationary conditions), the variations obtained are dramatic, with extended periods in which the flood quantile values are much higher than the existing design flood value. Therefore, using the existing design value to measure the snowmelt flood value under the combined influence of climate change and human activities will appear conservative in dry years, while in wet years, especially in the years of major floods, there may be greater risks than expected.