Non-Stationary Flood Frequency Analysis Using Cubic B-Spline-Based GAMLSS Model

: Under changing environments, the most widely used non-stationary flood frequency analysis (NFFA) method is the generalized additive models for location, scale and shape (GAMLSS) model. However, the model structure of the GAMLSS model is relatively complex due to the large number of statistical parameters, and the relationship between statistical parameters and covariates is assumed to be unchanged in future, which may be unreasonable. In recent years, nonparametric methods have received increasing attention in the field of NFFA. Among them, the linear quantile regression (QR-L) model and the non-linear quantile regression model of cubic B-spline (QR-CB) have been introduced into NFFA studies because they do not need to determine statistical parameters and consider the relationship between statistical parameters and covariates. However, these two quantile regression models have difficulties in estimating non-stationary design flood, since the trend of the established model must be extrapolated infinitely to estimate design flood. Besides, the number of available observations becomes scarcer when estimating design values corresponding to higher return periods, leading to unreasonable and inaccurate design values. In this study, we attempt to propose a cubic B-spline-based GAMLSS model (GAMLSS-CB) for NFFA. In the GAMLSS-CB model, the relationship between statistical parameters and covariates is fitted by the cubic B-spline under the GAMLSS model framework. We also compare the performance of different non-stationary models, namely the QR-L, QR-CB, and GAMLSS-CB models. Finally, based on the optimal non-stationary model, the non-stationary design flood values are estimated using the average design life level method (ADLL). The annual maximum flood series of four stations in the Weihe River basin and the Pearl River basin are taken as examples. The results show that the GAMLSS-CB model displays the best model performance compared with the QR-L and QR-CB models. Moreover, it is feasible to estimate design flood values based on the GAMLSS-CB model using the ADLL method, while the estimation of design flood based on the quantile regression model requires further studies.


Introduction
Flood frequency analysis is very important for the construction of hydrological projects. The stationary assumption has served as the basic assumption in flood frequency analysis for decades. However, due to climate change and human activities, the spatial and temporal distribution of rainfall and the catchment conditions have been changed. The stationary assumption has been challenged by many researchers [1][2][3][4][5][6][7][8][9][10][11][12][13] and the rationality of the design results obtained by traditional stationary flood frequency analysis has been questioned [6,14]. Therefore, the nonstationary frequency analysis of flood series under changing environments is of great significance to ensure the rationality of flood design results [1,[11][12][13][14][15]. There are several well-written articles summarizing the existing methods for non-stationary flood frequency analysis (NFFA) [8,[16][17][18][19]. NFFA has become one of the research hotspots in the field of flood frequency analysis under changing environments.
The generalized additive models for location, scale and shape (GAMLSS) model is currently the most commonly used method in NFFA [5,6,[20][21][22][23][24][25]. However, due to the large number of statistical parameters and the complexity of the relationship between statistical parameters and covariates, the model structure of GAMLSS is relatively complex. At the same time, the GAMLSS model assumes that the relationship between statistical parameters and covariates will remain unchanged in future, which may be unreasonable. In recent years, nonparametric methods have received more and more attention in the field of hydrological analysis and calculation. The simple linear quantile regression (QR-L) model in particular has been employed in NFFA. This is mainly because the QR-L model does not need to determine statistical parameters and consider the relationship between statistical parameters and covariates, which can definitely simplify the process of model construction. Koenker and Basset [26] first proposed the quantile regression method, and then it was introduced into the field of hydrological frequency analysis. Barbosa [27] used quantile regression to analyze Baltic Sea level change and found that the slope at the maximum value was more significant. Mazvimavi [28] used the quantile regression method to analyze changes in the rainfall series in Zimbabwe, and then found that climate change effects were not yet statistically significant within the time series of total seasonal and annual rainfall in Zimbabwe. Wang et al. [29] analyzed the possible changes in monthly precipitation in the southern United States using quantile regression. Feng et al. [30] used the quantile regression method to analyze the variation characteristics of the annual precipitation and runoff series in the Luanhe River Basin, and found that the annual runoff series in the sub-basin was no longer stationary.
However, in the field of hydrological frequency analysis, the dependence between the covariate and the independent variable is complex so it is unreasonable to simply use linearity, and more complex statistical relations between covariates and independent variables need to be considered [31]. The non-linear quantile regression model of cubic B-spline (QR-CB) was recommended in NFFA by Nasri et al. [31]. Compared with the QR-L model, the QR-CB model is more reasonable. The construction of the cubic B-spline function is only related to the number and position of the nodes and the degree of freedom of the function, and is not affected by the variables, so the model is more robust. Hendricks and Koenker [32] proposed spline parameterization using conditional quantile functions to estimate household electricity demand in metropolitan areas of Chicago. Nasri et al. [33] established a Generalized Extreme Value (GEV) model with a cubic B-spline curve under the Bayesian framework, and suggested that in future we can focus on the comparison of the extreme value model with the regression quantiles in order to use different covariates in quantile estimation. Nasri et al. [31] used cubic B-spline quantile regression to perform a non-stationary hydrological frequency analysis of the annual maximum and minimum flow records for Ontario, Canada.
However, it is difficult to estimate flood design values based on either the QR-L or QR-CB models. Currently, there is no accurate method for estimating the design flood values based on the quantile regression model. On the contrary, researchers have developed several non-stationary design methods when using the GAMLSS model [34][35][36][37][38][39]. Yan et al. [38] compared different design methods and recommended both equivalent reliability (ER) and average design life level (ADLL) for practical use because the design floods estimated by these two methods are linked with the design life of projects and possess reasonable confidence intervals.
Some researchers have tried to build a non-linear GAMLSS model for NFFA, and the results have shown that compared with the stationary model, variation types such as cubic spline function and parabolic function possess a better performance [22,40,41]. In this study, we attempted to develop a cubic B-spline-based GAMLSS model (GAMLSS-CB) by combining the GAMLSS model with the cubic B-spline. In the GAMLSS-CB model, the relationship between statistical parameters and covariates was fitted by the cubic B-spline under the GAMLSS model framework. We then compared its difference from the QR-L and QR-CB models. In this paper, the annual maximum flood series of four representative stations in the Weihe River Basin and the Pearl River Basin were selected as the research objects. The flood series of these stations are representative because they are located in different climate regions of China and exhibited either increasing or decreasing trends. Finally, based on the optimal non-stationary model, we also estimated the non-stationary design flood values using the ADLL method proposed by Yan et al. [25,38].
The paper is organized as follows: the second section introduces the methodologies, the third section gives the study areas and data, the fourth section gives the results, the fifth section is the discussion, and the sixth section is the conclusions.

Mann-Kendall Trend Test Method
The Mann-Kendall nonparametric trend test method was used to detect the long-term change trend of the precipitation and runoff series. This method has no need for the sample series to obey a specific distribution, and is not disturbed by a few outliers. It has been widely used in the trend analysis of hydrological and meteorological data. The Mann-Kendall test method [30,42] is as follows: Let i Y ( 1, 2, ..., ) i n = be a random variable for hypothesis testing; n represents the observed length of the sample, and the standardized test statistics are:

The Linear Quantile Regression (QR-L) Model
The QR-L model is based on the conditional quantile of the dependent variable Y required to regress the independent variable X and thus obtain the regression model under all quantiles. Linear quantile regression is related to linear least square regression and can be used to study the linear relationship between dependent variables and one or more independent variables. The accuracy of the parameter estimation can be independent of the distribution of the sample data, and can provide a more comprehensive description of the data from different quantile points, which can accurately describe the independent variable X for the dependent variable Y of the variation range and the conditional distribution of the effect of the shape [30].
Assuming that the distribution function of the random variable Y is ( ) According to different quantiles τ , different quantile functions ω ω n n y y , the regression estimate for the τ -quantile is solved, where ( ) τ ρ ⋅ is an asymmetric loss function: The parameter estimates can be obtained by the following formula, and a set of ( ) β τ values is determined from a τ value:

The Non-Linear Quantile Regression Model of Cubic B-Spline (QR-CB) Model
The nonparametric quantile model allows the linear hypothesis to be relaxed and the optimal model can be determined based on the data distribution [31]. Currently the most popular nonparametric quantile method is spline regression, which can be smoothed by adjusting the number of nodes. This paper considers constructing a QR-CB model.
Assuming that the distribution function of the random variable Y is: where i P is the control vertice, ,3 ( ) i B t is the harmonic function (base function) of the cubic B-spline, and the general formula of the basis function , ( ) k n B t in the n-time B-spline is [42,43]: If n + m + 1 vertices i P ( 0, 1, 2, ..., ) = + i n m are given, a parameter curve of m + 1 segments n times can be defined. Therefore, the basis function ,3 ( ) i B t of this paper is specifically:

Model Definition
In the GAMLSS model, it is assumed that the observation value is the distribution statistical parameter vector corresponding to time t, f is the number of distribution parameters, and n is the number of observations [44]. Let ( ) k k g θ denote the function relationship between k θ and the corresponding covariate k Y , which is generally expressed as: where k η is the vector of length n, 1 2 , , , ( ) is the regression parameter vector of represents the random effect of the j th term [44], namely the functional dependence of the distribution parameters on explanatory variables jk γ . The dependence can be linear and also smooth [40]. Adding the smoothing term in Formula (10) can identify non-linear dependence when modeling the parameter distribution. In this study the smooth dependence is based on cubic B-spline functions.
The first two parameters 1 θ and 2 θ of model (10) are usually defined as the location parameter vector and the scale parameter vector. If there are other parameters in the distribution, they are defined as shape parameters [44]. If we do not consider the effect of random effects, then ( ) For the location parameter μ , the scale parameter σ and the shape parameter υ , the full-parameter model that takes the time t as a covariate without considering the random effect is: where , , t t t μ σ υ are time-varying location parameters, scale parameters, and shape parameters, respectively, which can reflect the variation characteristics of non-stationary series statistical parameters with time. The model parameters β are estimated by the RS method available in the gamlss package [45], and the parameters and independent variables are fitted by the cubic B-spline function.

Model Evaluation Criteria
This paper used the generalized Akaike information criterion (GAIC) as the GAMLSS model fitting evaluation index. The GAIC criterion is based on the concept of entropy, which can weigh the complexity of the model and the superiority of the model fitting effect, and is a commonly used model evaluation index. The GAIC calculation formula is: is the global fitting deviation of the GAMLSS model, f d is the overall degree of freedom of the model, the penalty factor taking # = 2 represents the Akaike information criterion (AIC) value, and the model with the smallest AIC value is taken as the optimal model. The residual distribution of the optimal model is analyzed by a normal quantile-quantile (QQ) graph. The normal QQ graph is drawn in the plane coordinate system with the empirical residual as the ordinate and the theoretical residual as the abscissa. The smaller the deviation of the data point from the 1:1 line, the better the performance of the model.

Model Probability Coverage Test
The performance of the QR-L, QR-CB, and GAMLSS-CB models was qualitatively analyzed according to the magnitude of the model probability coverage bias value. The model probability coverage first needs to calculate the ratio of sample points in the coverage of each quantile curve to the total number of sample points, and then determine the difference between this ratio and the quantile. The smaller the difference, the better the model performance.

Filliben Test
The method of determining the optimal distribution pattern of the flood series through the Filliben correlation coefficient is more convenient and reliable. The optimal fitting distribution of the series is determined by the size of the Filliben correlation coefficient. The larger the Filliben value, the better the model performance [38].
Assuming that the actual distribution of the normalized residual 1 2 , , ..., n r r r obeys the normal distribution, the ascending statistic is (1)

Design Flood Value
Estimating the return period can be done according to 1 / m p = under traditional stationary conditions, where p represents the exceedance probability of the cumulative probability distribution function, and the corresponding design flood value formula is represents the inverse function of the cumulative probability distribution function. For a given design value, the probability distribution function obeyed by the flood extremum series can be estimated from the historical observation sample points, the exceedance probability is determined by the function curve, and the design return period corresponding to the given flood event is estimated [38]. When using the GAMLSS-CB model to estimate the design value, for a given return period, there is a design value corresponding to each year, which is also difficult to apply in practical engineering. Therefore, this article uses ADLL estimates of the design flood values based on the GAMLSS-CB model. For projects to be built with a design period of 1 2 T T − ( 1 T is the project start year and 2 T is the project termination year), the annual average reliability The ADLL method considers that the annual average reliability of a design value under nonstationary conditions should be equal to the annual reliability

Study Areas and Data
Multi-year maximum flood series were selected from the Xianyang and Huaxian stations of the Weihe River and the Gaodao and Dahuangjiangkou stations in the Pearl River Basin, these being four important stations (Table 1). The flood series of these stations are representative because they are located in different climate regions of China and exhibited either increasing or decreasing trends. The basic overview of the four stations is shown in Table 1 and Figure 1 below.
The Weihe River Basin originates in Bird Rat mountain in Weiyuan County, Gansu Province, and flows through Gansu and Shaanxi provinces to the Yellow River in Tongguan County, Weinan City (Figure 1a). The river is 818 km long, with a basin area of 134,766 km 2 and geographical coordinates at 33°40′-37°26′ N and 103°57′-110°27′ E [47]. The interannual variation of runoff in the middle and lower reaches of the Weihe River is characterized by small southern and large northern runoff. The mainstream flow of the Weihe River is the largest in autumn, accounting for 38% to 40% of annual runoff, 32.8% to 34.2% in summer, 17.7% to 19.1% in spring, and 8.3% to 9.9% in winter. The rainfall in the middle and lower reaches of the Weihe River is concentrated in July, August and September, and there are many heavy rains and flood disasters.
The Pearl River originates in Maxiong Mountain, Qujing City, Yunnan Province, and flows through Yunnan, Guizhou, Guangxi, Guangdong, Hunan, Jiangxi, and northern Vietnam. It is injected into the South Sea from eight downstream estuaries, with a total length of 2320 km and a basin area of 45,369 km 2 (Figure 1b). The Pearl River Basin consists of four river systems including Xijiang, Beijiang, Dongjiang, and the Zhujiang Delta [48]. The average annual rainfall of the basin is 1200-2200 mm, and the annual runoff is more than 330 billion m 3 . The flood season accounts for 80% of the total annual runoff from April to September, and accounts for more than 50% of the annual runoff in summer (June-August).

Mann-Kendall Trend Analysis
A Mann-Kendall test was carried out on the historical data from each station using the trend package in the R language. From this, the P, |Zc|, and S values of the Huaxian, Gaodao, Dahuangjiangkou, and Xianyang stations were obtained. The P, |Zc|, and S values of the four stations are shown in Table 2 below. and Xianyang stations were significant at a 5% significance level, while the trend of Dahuangjiangkou station was significant at a 10% significance level. Moreover, the positive and negative relationships among the S values allow us to conclude that Gaodao and Dahuangjiangkou stations showed an increasing trend, while the other stations showed a decreasing trend. Huaxian and Xianyang stations showed a significantly decreasing trend, while Gaodao station showed a significantly increasing trend. Finally, Dahuangjiangkou station showed no significant upward trend. Figure 2 shows the linear trend line of the annual maximum flood series at each station.

Determination of Optimal GAMLSS-CB Model
In order to select the optimal probability distribution type, the annual maximum flood series of four stations in the Weihe River and the Pearl River Basin located in different climate regions of China was selected as the research object. The time t was the covariate, and the relationship between the statistical parameters and the covariate was the cubic B-spline function; AIC was the evaluation criterion, and the gamma distribution (two-parameter), lognormal distribution (two-parameter) and GEV distribution were compared respectively. The shape parameter of GEV is sensitive and difficult to estimate; thus, it is assumed to be constant in this study in keeping with other studies [49,50]. The normal QQ graph can be used to judge the performance of the GAMLSS-CB optimal model. For each station, a total of 12 non-stationary models were constructed considering the combination of distribution types and variation types for each location and scale parameter. The corresponding AIC values are shown in Table 3 below. Table 3. Akaike information criterion (AIC) values of the non-stationary models for each station. Note that letter "L" in the models' names represents location parameter and "S" represents scale parameters. Number "0" means the parameter is invariant while "t" means the parameter varies with time covariate. Besides, the AIC value in bold is the optimal model for each station. It can be concluded that the gamma distribution is the optimal distribution when using the GAMLSS-CB model. The non-stationary gamma distribution with both location parameters and scale parameters changing with time had the best performance for Huaxian and Xianyang stations: the AIC values were 1054.82 and 960.70. However, for Gaodao and Dahuangjiangkou stations, the optimal models were non-stationary gamma distribution with location parameters changing with time and the scale parameters remaining unchanged: the AIC values were 1060.90 and 1157.02. Figure  3 shows the QQ map of the optimal non-stationary model for each hydrological station. The results show that the optimal non-stationary model empirical residual and theoretical residual are stationary and distributed near the 1:1 line, indicating that the model has a good performance.

Qualitative Analysis of Model Performance
Qualitative analysis of the probabilistic coverage rate of models was conducted by calculating the quantile curve of each station using the QR-L, QR-CB, and GAMLSS-CB models; see the following Figures 4-6 for details. Table 4 below shows the probabilistic coverage rate of non-stationary models for each station. It can be concluded that in the QR-L model the probability coverage deviation of the Huaxian station model was 0.81%-4.84%, the deviation of Gaodao station was 0.08%-5.74%, the deviation of Dahuangjiangkou station was 0-0.36%, and the deviation of Xianyang station was 0-2.59% (Figure 4). The QR-CB model shows that the probability coverage deviation of Huaxian station was 0-2.42%, the deviation of the Gaodao station model was 0.08%-2.87%, the deviation of Dahuangjiangkou station was 0-1.79%, and the deviation of Xianyang station was 0.17%-5.17% ( Figure 5). The GAMLSS-CB model probability coverage deviation of Huaxian station was 0.16%-10.48%, the deviation of Gaodao station is 0.08%-6.97%, the deviation of Dahuangjiangkou station was 0.36%-5.36%, and the deviation of Xianyang station was 0.17%-4.31% ( Figure 6). Based on the overall data, the performance of the QR-L and QR-CB models was basically the same, while the performance of the GAMLSS-CB model was slightly weaker.

Quantitative Analysis of Model Performance
We quantitatively compared the performance of the model by calculating the Filliben correlation coefficient. The model Filliben correlation coefficient is shown in Table 5 below. According to the principle that when the Filleben correlation coefficient is larger, the performance is better, we found that the performance of each station model was good, especially the GAMLSS-CB model, which had the best model performance compared with the QR-L and QR-CB models, and that the performance of the QR-L model was better than the QR-CB model. Since the accuracy of qualitative analysis is affected by the rules of artificial counting, and the accuracy of quantitative analysis is more secure, in this study quantitative analysis was the main method and qualitative analysis was auxiliary. This study found that the GAMLSS-CB model had the best model performance compared with the QR-L and QR-CB models, based on qualitative analysis and quantitative analysis. Then the ADLL method was used to estimate the non-stationary design flood value of the GAMLSS-CB model.

Design Values of GAMLSS-CB Model
This study concluded that the GAMLSS-CB model had the best model performance compared with the QR-L and QR-CB models. Therefore, this study assumed the engineering design period to be 50 years, from 2015 to 2064, and used the GAMLSS-CB model to estimate the flood design value. Figure 7 below shows that the flood design values estimated by the ADLL method based on the GAMLSS-CB model were reasonable and reliable. It can be used for non-stationary engineering design due to its linkage with the design period of projects under changing environments.

Discussion
Currently, there are very few studies estimating design flood value based on the quantile regression model, largely because the design results are affected by the distribution of sample points when estimating design flood values based on the quantile regression model. The number of available observations in particular becomes lower when estimating design values corresponding to higher return periods, leading to unreasonable and inaccurate design values. For this reason, this study does not compare the design flood values estimated using the quantile regression model with those using other models. In future studies, more research efforts are needed to improve the accuracy of design values based on quantile regression. For example, to avoid the influence of rare sample points of extreme floods when estimating design values with higher return period, we can turn to use the "peak-over-threshold" (POT) sampling method to increase the sample size of extreme floods. Different from the annual maximum sampling method of selecting only one sample per year, the sample size can be expanded by selecting 2-3 high-quantile flood samples exceeding the threshold using the POT sampling method, thus solving the problem of scarce high-quantile flood samples in estimating design values based on quantile regression.
The GAMLSS-CB model has a good performance, and its method of estimating design flood value is stable and reliable. Therefore, the accuracy of the design flood value estimated by the model is guaranteed, but the currently constructed GAMLSS-CB model only considered time as a covariate, which has certain limitations. In later research, we should try to add population, climate, etc. as covariates to optimize the GAMLSS-CB model.
In this study, only the point estimation of design quantiles was presented, which lacks standard error estimation. Nasri et al. [33] developed the GEV-B-Spline model and evaluated the uncertainties of design quantiles using the bias (BIAS) and the root mean square error (RMSE). They found that the Bayesian estimation for the GEV-B-Spline model can achieve satisfactory results. Therefore, in future research, we will consider using Bayesian estimation or other methods to improve the accuracy of parameter estimation and explore the uncertainties of design flood values.

Conclusions
This paper constructed GAMLSS-CB, QR-CB and QR-L models based on the annual maximum flood series of four stations in the Weihe River and the Pearl River basin located in different climate regions of China, and compared the model performance of different non-stationary models using probability coverage and Filliben correlation coefficient. In addition, the design flood values under changing environments were also estimated using the optimal model. The main conclusions of this study are obtained as follows: (1) Through the Mann-Kendall trend test, it is concluded that both Huaxian station and Xianyang station showed a significantly decreasing trend, while Gaodao station showed a significantly increasing trend. In addition, Dahuangjiangkou station showed no significant upward trend. (2) The gamma distribution is the optimal distribution when using the GAMLSS-CB model. The non-stationary gamma distribution with both location parameters and scale parameters changing with time had the best performance for Huaxian and Xianyang stations, while for Gaodao and Dahuangjiangkou stations the optimal models were non-stationary gamma distribution with location parameters changing with time and the scale parameters remaining unchanged.