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

Abstract

With the acceleration of human economic activities and dramatic changes in climate, the validity of the stationarity assumption of flood time series frequency analysis has been questioned. In this study, a framework for flood 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 flood time series in the Kenswat Reservoir control basin of the Manas River. The study shows that there are clear non-stationarities in the flood regime, and the characteristic series of snowmelt flood shows an increasing trend with the passing of time. The parameters of the flood 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 flood conditions, which can better describe the dynamic process of snowmelt flood characteristic series. Compared with the design flood 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.

1. 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 non-stationary 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₁, P₃ and T₇₈) 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.

2. 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² and annual average runoff is 12.21 × 10⁸ m³. 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², 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⁸ m³ [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.

3. Data and Methods

3.1. 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₁, P₃, P₅, P₇, P₁₅ and P₃₀) before the occurrence of the flood peak were collected at Kenswat Hydrological Station from 1957 to 2006 as basic data.

3.2. Methods

3.2.1. 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].

3.2.2. 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,\dots ,n$ obey the probability density function $f\left(y_i|{\theta }^i\right)$ where ${\theta }^i=\left({\theta }_{1i},{\theta }_{2i},{\theta }_{3i},{\theta }_{4i}\right)=\left({\mu }_i,{\sigma }_i,{\upsilon }_i,{\tau }_i\right)$ is the distribution (statistic) parameter vector corresponding to a certain moment, n is the number of observations, and p is the number of distribution (statistic) parameters. The parameters ${\mu }_i$ and ${\sigma }_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 ${\upsilon }_i$ and ${\tau }_i$, respectively. Let $g_k(•)$ denote ${\theta }_k$ a monotonic function relationship with the corresponding explanatory variable $X_k$ and the random effect term, generally expressed as:

$$g_k\left({\theta }_k\right)={\eta }_k=X_k{\beta }_k+{\displaystyle \sum _{j=1}^{J_k}{Z}_{jk}{\gamma }_{jk}},$$

(1)

where ${\theta }_k$ are vectors of length n; ${\beta }_k={\left({\beta }_{1k},{\beta }_{2k},\cdots ,{\beta }_{I_{k}k}\right)}^T$, which is a regression parameter of length $I_k$ Vector; $X_k$ is an explanatory variable matrix of $n\times I_k$; $Z_{jk}$ is a known fixed design matrix of $n\times q_{jk}$; ${\gamma }_{jk}$ is a normal distribution random variable vector of $q_{jk}$ dimension; $Z_{jk}{\gamma }_{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 θ₁ and scale parameter θ₂ with covariate is used.

In this paper, we considered as candidates four widely used two-parameter distribution functions in modelling streamflow data (

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
LOGNO	$\begin{array}{l}{f}_Y\left(y|\mu ,\sigma \right)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\frac{1}{y}\exp \left\{-\frac{{\left[\log \left(y\right)-\mu \right]}^2}{2{\sigma }^2}\right\}\\ y>0,\mu >0,\sigma >0\end{array}$	$\begin{array}{l}E\left(Y\right)={\omega }^{1/2}{e}^{\mu }\\ Var\left(Y\right)=\omega \left(\omega -1\right)e^{2\mu }\\ \omega =\exp \left({\sigma }^2\right)\end{array}$
GA	$\begin{array}{l}{f}_Y\left(y|\mu ,\sigma \right)=\frac{1}{{\left({\sigma }^2\mu \right)}^{1/{\sigma }^2}}\frac{y^{1/{\sigma }^2-1}{e}^{-y/\left({\sigma }^2\mu \right)}}{\Gamma \left(1/{\sigma }^2\right)}\\ y>0,\mu >0,\sigma >0\end{array}$	$\begin{array}{l}E\left(Y\right)=\mu \\ Var\left(Y\right)={\sigma }^2{\mu }^2\end{array}$
GU	$\begin{array}{l}{f}_Y\left(y|\mu ,\sigma \right)=\frac{1}{\sigma }\exp \left[\left(\frac{y-\mu }{\sigma }\right)-\exp \left(\frac{y-\mu }{\sigma }\right)\right]\\ -\infty <y<\infty ,-\infty <\mu <\infty ,\sigma >0\end{array}$	$\begin{array}{l}E\left(Y\right)\cong \mu -0.57722\sigma \\ Var\left(Y\right)={\pi }^2{\sigma }^2/6\end{array}$
PIG	$\begin{array}{l}{p}_Y\left(y|\mu ,\sigma \right)={\left(\frac{2\alpha }{\pi }\right)}^{\frac{1}{2}}\frac{{\mu }^y{e}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\sigma $}\right.}{K}_{y-\frac{1}{2}}\left(\alpha \right)}{{\left(\alpha \sigma \right)}^{y}y!}\\ y=0,1,2,\cdots ,\infty ,\mu >0,\sigma >0\end{array}$	$\begin{array}{l}{\alpha }^2=\frac{1}{{\sigma }^2}+\frac{2\mu }{\sigma }\\ K_{\lambda }\left(t\right)=\frac{1}{2}{\displaystyle {\int }_0^{\infty }{x}^{\lambda -1}\exp \left\{-\frac{1}{2}t\left(x+x^{-1}\right)\right\}}\end{array}$

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

The GD of GAMLSS model is defined as follows:

$$GD=-2l({\widehat{\theta }}_i),i=1,2,3,4$$

(2)

where $l\left({\widehat{\theta }}_i\right)$ 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:

$$GAIC=GD+\#df$$

(3)

where $df$ 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 \left(n\right)$ (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].

4. Results

4.1. Correlation Analysis of Temperature and Snowmelt Flood

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.

	Qmax	Wmax1	Wmax3	Wmax7	Wmax15	Wmax30
r	0.3673	0.3670	0.3818	0.4030	0.4186	0.4290
p-value	0.0087	0.0088	0.0062	0.0037	0.0025	0.0019

summarizes the Pearson correlation coefficients of the mean temperature series in July and August (T₇₈) 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.

4.2. Correlation Analysis between Precipitation in Early Stage and Snowmelt Flood

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₁₅ and P₃₀, there is a strong correlation between the previous-period 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₃, while the W_max1, W_max3, W_max7, W_max15 and W_max30 have the best correlation with the P₁, and all the correlations can be tested by a 0.05 significance level. The correlation coefficients between the P₁ 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₁ and the W_max7, W_max15 and W_max30, the P₃ 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₁ and P₃, P₇, P₁₅ and P₃₀ show a downward trend, among which the correlation coefficient between P₁ and P₃₀ is only 0.255, and the correlation coefficient between P₃ and P₃₀ is 0.610. Therefore, the P₁ and P₃ with the strongest correlation with the snowmelt flood characteristic series are selected as the influencing factors of the model.

4.3. 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. The Akaike information criterion (AIC) was used to determine the fit of Generalized Additive Models for Location, Scale, and Shape (GAMLSS) under different probability distributions of snowmelt flood time series. Lower AIC represent a better performance.

Distribution	Qmax	Wmax1	Wmax3	Wmax7	Wmax15	Wmax30
LOGNO	630.94	804.73	890.34	948.26	1010.42	1060.20
GA	638.12	811.39	895.91	952.86	1016.10	1064.42
GU	696.45	872.56	948.20	1001.48	1075.57	1115.66
PIG	633.01	1183.23	1307.98	1387.78	1534.77	1582.77

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. Residual distribution moments and Filliben coefficients for optimal fitting distribution of snowmelt flood time series (model 0).

Snowmelt Flood Characteristic Series	Distribution	Mean	Variance	Skewness	Kurtosis	Filliben Coefficient
Qmax	LOGNO	0	1.020	0.919	3.596	0.968
Wmax1	LOGNO	0	1.020	0.958	3.976	0.961
Wmax3	LOGNO	0	1.020	0.909	3.918	0.959
Wmax7	LOGNO	0	1.020	0.864	3.927	0.962
Wmax15	LOGNO	0	1.020	1.072	5.089	0.959
Wmax30	LOGNO	0	1.020	0.851	4.510	0.966

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 each model that the residuals of each model cannot meet the requirements of obeying the normal distribution.

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.

4.4. Based on the Results of the Non-Stationary Model with Time Variables: Model 1

4.4.1. Model Fitting Evaluation

The non-stationary model is constructed with time t as the explanatory variable, and the cumulative probability distribution parameters θ₁ and θ₂, 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. 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.

Distribution	θ1	θ2	GD	AIC	SBC	θ1	θ2	GD	AIC	SBC
	Qmax					Wmax1
LOGNO	t	t	620.22	628.22	635.86	t	t	791.88	799.88	807.53
GA	t	t	626.19	634.19	641.84	t	t	795.94	803.94	811.59
GU	t	t	677.45	685.45	693.10	t	ct	856.48	862.48	868.21
PIG	t	ct	626.37	632.37	638.10	t	ct	1175.83	1181.83	1187.57
	Wmax3					Wmax7
LOGNO	t	t	878.16	886.16	893.81	t	t	935.70	943.70	951.35
GA	t	t	881.42	889.42	897.07	t	t	938.05	946.05	953.70
GU	t	ct	931.93	937.93	943.67	t	ct	983.88	989.88	995.62
PIG	t	ct	1300.30	1306.30	1312.04	t	ct	1378.81	1384.81	1390.54
	Wmax15					Wmax30
LOGNO	t	t	996.86	1004.86	1012.51	t	t	1045.25	1053.25	1060.90
GA	t	t	999.76	1007.76	1015.41	t	t	1046.86	1054.86	1062.51
GU	t	ct	1058.44	1064.44	1070.18	t	ct	1097.66	1103.66	1109.39
PIG	t	ct	1525.68	1531.68	1537.42	t	ct	1574.07	1580.07	1585.81

. For the Q_max, the LOGNO distribution is the best fitting distribution. The distribution parameter θ₁ exhibits a linear dependence on time t in the four candidate distributions models. The distribution parameter θ₂ 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.

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 θ₁ exhibits a linear dependence on time t, and the distribution parameter θ₂ is constant and independent of the time trends.

Figure 4 and

Table 6. Residual distribution moments and Filliben coefficients for optimal fitting distribution of snowmelt flood time series (model 1).

Snowmelt Flood Characteristic Series	Distribution	Mean	Variance	Skewness	Kurtosis	Filliben Coefficient
Qmax	LOGNO	0	1.020	0.830	3.562	0.977
Wmax1	LOGNO	0	1.020	0.543	2.563	0.977
Wmax3	LOGNO	0	1.020	0.476	2.577	0.978
Wmax7	LOGNO	0	1.020	0.379	2.420	0.979
Wmax15	LOGNO	0	1.020	0.457	2.921	0.982
Wmax30	LOGNO	0	1.020	0.180	2.665	0.990

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.

4.4.2. 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.

4.5. 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₁, P₃ and T₇₈) 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. 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.

Distribution	θ1	θ2	GD	AIC	SBC	θ1	θ2	GD	AIC	SBC
	Qmax					Wmax1
LOGNO	T78 + P3	t	610.08	618.08	625.73	T78 + P1	P1	765.05	775.05	784.61
GA	T78 + P3	ct	612.61	620.61	628.26	T78 + P1	P1	764.72	774.72	784.28
GU	T78 + P3	ct	661.83	669.83	677.48	T78 + P1	P1	779.19	789.19	798.75
PIG	P3	ct	621.21	627.21	632.94	T78	P1	1164.72	1172.72	1180.37
	Wmax3					Wmax7
LOGNO	T78 + P1	P1	856.04	866.04	875.60	T78 + P1	P1	919.13	929.13	938.69
GA	T78 + P1	P1	855.85	865.85	875.41	T78 + P1	P1	919.56	929.56	939.12
GU	T78 + P1	T78 + P1	865.91	877.91	889.38	T78 + P1	P1	933.82	943.82	953.38
PIG	T78 + P1	ct	1278.51	1286.51	1294.15	P1	ct	1371.10	1377.10	1382.84
	Wmax15					Wmax30
LOGNO	T78 + P1	P1	983.80	993.80	1003.36	T78 + P1	T78	1035.10	1045.10	1054.66
GA	T78 + P1	P1	984.95	994.95	1004.51	T78 + P1	T78	1035.91	1045.91	1055.47
GU	T78 + P1	P1	1005.48	1015.48	1025.04	T78 + P1	P1	1056.75	1066.75	1076.31
PIG	T78	ct	1521.19	1527.19	1532.99	ct	ct	1578.77	1582.77	1586.60

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.

For the Q_max, the LOGNO distribution is the best fitting distribution. Among the three candidate indicators of T₇₈, P₁, and P₃, P₁ 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₇₈ and P₃ indicators tend to describe the linear variation of the distribution parameter θ₁, 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₃ 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₃. The T₇₈ 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₇₈. The distribution parameter of all distribution models is constant, indicating that the variance of the Q_max has little to do with climatic factors.

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₇₈, P₁, and P₃, P₃ 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₇₈ indicator of the W_max1, W_max7 and W_max15 are mainly expressed as the linear dependence of the distribution parameter θ₁, while the P₁ indicator is mainly expressed as the linear dependence of the distribution parameters θ₁ and θ₂. 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₇₈ indicator of the W_max3 and W_max30 are all expressed as the linear dependence between the distribution parameters θ₁ and θ₂, and the P₁ indicator is also expressed as the linear dependence between the distribution parameters θ₁ 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. Residual distribution moments and Filliben coefficients for optimal fitting distribution of snowmelt flood time series (model 2).

Snowmelt Flood Characteristic Series	Distribution	Mean	Variance	Skewness	Kurtosis	Filliben Coefficient
Qmax	LOGNO	0	1.020	0.352	3.425	0.985
Wmax1	GA	0	1.011	0.007	1.896	0.986
Wmax3	GA	0	1.019	0.015	1.863	0.987
Wmax7	LOGNO	0	1.020	0.082	1.795	0.985
Wmax15	LOGNO	0	1.020	0.246	1.916	0.983
Wmax30	LOGNO	0	1.020	0.108	2.194	0.992

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.

4.5.2. Analysis of Optimal Model Fitting Results

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. Summary of the results of modeling the snowmelt flood time series with model 2 under non-stationary conditions. The results show estimates of the 90th, 95th and 98th percentiles. The design flood values for the 50-, 20-, and 10-year return periods of Kenswat Reservoir (approved by China Renewable Energy Engineering Institute in December 2008).

Snowmelt Flood Characteristic Series	Extremum	Year	Quantile of Snowmelt Flood Time Series			Design Standard Value
			98%	95%	90%	50-Year	20-Year	10-Year
Annual maximum peak discharge (m3/s)	maximum	1996	1459	1172	966	1249	856	600
	minimum	1972	351	341	328
Annual maximum 1-day flood volume (105 m3)	maximum	1996	11,636	10,141	8937	7406	5206	3756
	minimum	1972	2013	2215	2077
Annual maximum 3-day flood volume (105 m3)	maximum	1996	26,029	22,623	19,998	15,920	12,090	9425
	minimum	1972	6424	5984	5578
Annual maximum 7-day flood volume (105 m3)	maximum	1996	41,058	35,625	31,400	29,430	23,120	18,620
	minimum	1972	13,721	12,744	11,939
Annual maximum15-day flood volume (105 m3)	maximum	1996	69,712	61,205	54,513	56,830	44,230	35,340
	minimum	1972	26,262	24,281	22,759
Annual maximum 30-day flood volume (105 m3)	maximum	1996	89,702	80,799	74,637	93,383	74,599	61,042
	minimum	1972	44,883	42,446	40,431

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³ s⁻¹, 341~1172 m³ s⁻¹, 328~966 m³ s⁻¹, 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.

5. 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₇₈ and P₃ indicators. The location parameter of the annual maximum flood volume series depended on the T₇₈ and P₁ indicators. Moreover, the scale parameter can be related only to the P₁ 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.