Handling non-stationary flood frequency analysis using TL-moments approach for estimation parameter

Climate change is a change in the form of long-term seasonal patterns of temperature, precipitation, humidity, wind, and seasons happening in a period of time up to decades or millions of years. The changes of extreme weathers or climate events phenomena (e.g., heat waves and heavy downpours) influence the experiences of climate change that increasingly face natural disasters such as drought and floods. A report from the Intergovernmental Panel on Climate Change (IPCC) predicted that climate change can lead to sea warming and glacial meltdown that increases the sea level between 7 and 23 inches by 2,100 and possibly causing the Pacific islands to be flooded (Trenberth et al. 2007). From the report, a warmer climate is also expected to increase the risk of drought and more flood events due to increased rainfall. Malaysia is one of the countries that have been affected, where 18% of coastal areas are affected by climate change phenomena (Ismail 2018). Due to the increasing amount of precipitation as an effect of climate change, it is expected that there will be a contribution to the non-stationary characteristics of flood data.

The non-stationary flood series means processes that are not stationary and that have statistical properties with deterministic functions of time. Environmental factors that cause floods, such as climate, land cover, and global warming will be changing over time (i.e., inconstant trend) (Khaliq et al. 2006). The existence of non-stationarity in extreme flood events makes flood frequency analysis more difficult in resources management. This is because the changes of time or covariates in the parameter model makes the probability distribution model for floods become more complicated (Gado & Nguyen 2016a). These consist of method suitability in the change analysis and the response of climate impacts on water resources.

Many researchers have studied flood frequency analysis using non-stationary hydrological data. Khaliq et al. (2006) developed parametric values of different times in distribution parameters and quantiles of hydrological data. Villarini et al. (2009, 2010) demonstrated the generalized additive models for location, scale, and shape parameters as a tool to overcome the problem of slowly varying the time trends and abrupt changes in the annual peak records. Some researchers' experiences with non-stationary probability distribution models involve parameters that vary over time. Researchers implemented trends into the first and second moments of flood data to perform non-stationary flood frequency analysis (NFFA) (see Gado & Nguyen 2016a, 2016b; Mat Jan et al. 2018).

An appropriate parameter estimation method should be established to deal with flood frequency analysis under the non-stationary context. In most recent studies, parameter estimation has been measured using various techniques in order to handle NFFA. There were many methods of handling NFFA that used by researchers. Some researchers (see Ouarda & El Adlouni 2011) used generalized maximum likelihood analysis in order to avoid the numerical problems that can arise in the maximum likelihood estimators (MLE) method during the short series data parameter estimation. Gado & Nguyen (2016a, 2016b) introduced L-moments in the non-stationary version (LM-NS) by removing linear trends from the flood series and estimating the parameter and quantile of the generalized extreme value (GEV) distribution using the L-moments method. In their study, the MLE method as a comparison method investigated the performance of stationary and non-stationary GEV and Gumbel models using Monte Carlo simulation (Gado & Nguyen 2016a) and the application of real flood data from 32 annual maximum flow series in Quebec, Canada (Gado & Nguyen 2016b). From both experiments, LM-NS outperformed the MLE and stationary method in the simulation analysis for most cases. This method is also recommended in the model parameter estimation for the non-stationary background.

Similar to this study, the L-moments method is used as a comparison method for the model parameter estimation of GEV distribution since previous research proved that the L-moments method was a better method than the MLE method (Gado & Nguyen 2016a, 2016b). This paper compares the different ways in which a Trimmed L-moments (TL-moments) method is proposed as a main approach in parameter estimating of the GEV model. The TL-moments method has been used in many investigational studies of flood frequency analysis. Since the development of TL-moments as a generalization of L-moments by Elamir & Seheult (2003), many studies have been conducted regarding the TL-moments method (see Shahzad & Ashgar 2013; Mat Jan & Shabri 2017). However, the TL-moments method has not yet been used in NFFA. Based on the capabilities of the four levels of TL-moments, TL(η,0) (η = 1, 2, 3, 4) in the flood frequency analysis (see Mat Jan & Shabri 2017), this study proposes the TL-moments method as an innovative method for NFFA.

The purpose of this study is to investigate the behavior of the proposed TL-moments method in the presence of NFFA and apply the TL-moments method along with the GEV distribution. The performances of this method are evaluated by comparing the results of performances of the L-moment method for the stationary and non-stationary GEV models. The annual maximum flood series from Peninsular Malaysia and Monte Carlo simulation are utilized. For this intention, the error computation criteria are used to calculate the individual performance of the probability distribution function by achieving the minimum error of the relative root mean square error (RRMSE).

Four series of annual maximum streamflow, located in Johor, Malaysia, are used to view other aspects of performances between the comparison methods, L-moments and TL-moments method. These stations have been chosen for a minimum of 15 years, ending in 2009, to evaluate the effect of flood behavior. The collected data were obtained from the Department of Irrigation and Drainage, Ministry of Natural Resources and Environment, Malaysia. Table 1 presents the list of selected streamflow with their characteristics.

This study is conducted based on several methods comprising detection of a possible change point and non-stationarity of the streamflow. Then, a flood frequency analysis is conducted to estimate the parameters and quantile of the underlying distribution. In this analysis, the non-stationary model is developed and estimated by selected methods. Following the performance comparison by using the simulation and bootstrap-based study, it is necessary to assess the method ability and the accuracy of the non-stationary model to predict rare flood among river stations. The final step of this study is the selection of the best model to highlight the capabilities of the proposed method in the selection of realistic models based on the experimental data. The detailed explanation of each analysis of the study will be discussed herein. Figure 1 presents an overview of the study.

A hydrologic time series is said to be non-stationary if there is a trend or sudden changes in the mean and variance (or called change point) of the probability distribution in the time series (Tan & Gan 2015). Hence, a change point analysis is implemented in order to detect the difference between the changes from one period of time (years) to another period with the status remaining the same until the next shift occurs (Khaliq et al. 2009).

In this study, the Pettitt test (Pettitt 1979), a non-parametric test, is applied. This test is a non-parametric test that allows changes in the mean to be detected when the point of change is unknown. This test has been selected by many researchers in the detection of change point in hydrological series (see Tan & Gan 2015; Li et al. 2016). Under the null hypothesis H₀, the variable T follows the distributions with no change on the location parameter, while the alternative hypothesis is the existence of the change point in the series. The series of random variables (X₁, X₂, … X_T) exist with the change point at time τ when X_t (for t = 1, 2 … τ) has the same distribution function F₁(x), and X_t (for t = τ + 1 … τ) has the same distribution function F2(x), where F1(x) ≠ F2(x) (Jha & Sigh 2013).

The change point of the series is located at K_T, provided that the statistics is significant.

The significant change point of the hydrological analysis is said to be existing in the series if the probability is equal to or greater than 0.95, where the approximate significance probability for a change point is P = 1 − p. When there is a change point in a streamflow series, the series is separated into two subseries at the position of the change point.

Considering trends in datasets for the non-stationary analysis of flood frequency and magnitude will increase the number of the estimated parameters from fixed-length time series. The development of new models could gain dynamic progress of the flood flow properties over time (Villarini et al. 2009). The new statistical models enable capturing dynamically the change of probability density functions over time, in order to obtain strong flood estimation. The non-stationary models used in this study are developed, as seen in Table 2.

The location and scale are modeled as a linear function of time (GEV1 and GEV2) where time, t, is expressed as a covariate. The exponential function has been introduced in the scale parameter of the GEV2 model to keep the positivity of α (Yilmaz & Perera 2013). A non-stationary model with the mean as a parabolic time of function is named as the GEV3 model. The shape parameter in GEV1, GEV2, and GEV3 remained constant since it is difficult to estimate with precision and it is unrealistic to model it as a smooth function of time (Hasan et al. 2013). The parameters are estimated using TL-moments as the proposed method and the competitor method called the L-moments method. The process of parameter estimation is explained in the next section.

The best model among the considered non-stationary models is the one with the minimum value of the selected criterion. For the selection of the best model, this study uses the Bayesian information criterion (BIC) (Schwarz 1978). The value of BIC is calculated by using formula BIC = −2 l + m × ln(n), where the total of negative log-likelihood (−l) and the number of model parameters (m) are multiplied with the n observations of data in order to select the best model among the deterministic function (Sober 2002). The model with the minimum BIC value is considered to be the right model.

The method of L-moments and TL-moments has been increasingly used for the parameter estimation in hydrology. The popularity of the L-moments method is because of its robustness with respect to outliers in the data. Otherwise, the TL-moments method assigns zero weight to extreme observations, which are easy to compute when estimating from a sample containing outliers (Noura 2010). The formation of L-moments is from linear combinations of probability weighted moments (PWM) (Hosking 1990).

Equations (9)–(11) present the L-moments and TL-moments, TL(η,0) (η= 1, 2, 3, 4) of the GEV distribution. Thus, parameters ξ, α, and k for h= 0, 1, 2, 3, 4 are estimated by solving Equations (9)–(11) simultaneously. After the determination of the corresponding parameters of the GEV distribution for different levels (η= 1, 2, 3, 4) of TL-moments, the required return period is estimated from the quantile x(F) of the GEV distribution (Equation (7)).

Thus, in the context of the non-stationary case, the steps of the estimation parameter by TL-moments for the GEV2 model (location and scale modeled as a linear function of time) described after Strupczewski et al. (2016) are as follows:

• 1.

  The non-stationarity in the mean, ξ_t and standard deviation, α_t of the time series is estimated by the fitting linear model.

• 2.
  Transform the data to the standardized data so as to deprive of trends:

  

  (12)
• 3.

  Estimate the stationary parameter (ξ, α, k) and quantiles, x(F) of a chosen distribution by L-moments and TL-moments methods on the stationary sample.

• 4.
  Estimate the non-stationary quantile of the distribution by:

  

  (13)

  where F is the cumulative distribution. The parameters estimated by L-moments and TL-moments method for all models are summarized in Table 3.

The performance of different TL-moments levels of GEV distribution is evaluated in this section. In order to investigate the capabilities of the GEV distribution between L-moments and TL-moments method for non-stationary and stationary models (GEV0, GEV1, GEV2, and GEV3), the Monte Carlo simulation technique is used. The quantiles which correspond to the updated parameter have been estimated at the end of the time series (t = n).

In this Monte Carlo simulation, 1,000 samples are generated for the recurring intervals of 20 to 1,000 years, using different sample sizes of 15, 30, 50, and 100 years. Based on data from Peninsular Malaysia, the shape and mean (ξ₁) parameter of the GEV distribution is in the range of (−0.3, +0.3) to see the sensitivity of the methods to the values of different trends while the location parameter of the trend is set at ξ₀ = 0 for the GEV1 model. For the GEV2 model, ξ₀ = 0 and ξ₁ = −0.1 for time-dependent in mean. Time-dependent in standard deviations, α₀ = 0 and two cases of α₁ = ±0.1 are considered to see the effect on the positive and negative trends in the standard deviation of the method (Gado & Nguyen 2016a). For the GEV3 model, ξ₀ = 0, ξ₁ = 0.3, and ξ₂ = −0.005. For GEV1 and GEV3 models, the standard deviation, α, is equated to 1, respectively (Wang 1990), since it is the standard form of any distribution (NIST/SEMATECH 2012).

In addition, the residual bootstrap method is applied to test the accuracy of the comparison of GEV non-stationary models in rare floods, as studied by Gado & Nguyen (2016b) and Khaliq et al. (2006). This study generates a new sample from the original time series data using this bootstrap method. Thus, the parameters and quantiles of the non-stationary GEV model are estimated from the bootstrapped samples.

For a particular recurrence interval, N is the total number of samples, Qs is the simulated quantiles of the nth sample and Qc is the calculated quantiles from the observed data. The minimum achieved values are highlighted in order to better illustrate the result. In this study, boxplots are used as one of the tools for grouping the results based on their statistical properties; the aim is to select the most suitable method of GEV distribution based on the minimum RRMSE values achieved.

The four streamflow records detailed in the previous section are tested in two main properties. The Pettitt test previously stated is used to generate significance for the presence of a change point and the location of a possible change point. The results are presented in Table 4. As shown in the table, two streamflows (Johor and Linggui stations) recorded the significance of change point smaller than 0.95. Thus, these two streamflows show an insignificant change point while another two streamflows (Bekok and Kahang stations) have a change point bigger than 0.95, where Bekok station achieved 1.00 significance of change point.

The presence of the change point for the streamflow is also illustrated in Figure 2. From the plot, the significant change point stations (Bekok and Kahang) show the abrupt changes in the mean at 1989 for Bekok station and Kahang station at 1995. Johor station and Linggui station also showed the change point presence in the mean of the streamflow series. Since these stations are identified as not significant, only the figure of abrupt changes for Bekok station and Kahang station are presented in this study.

Based on the comparison between the two tests (change point detection and non-stationary test), the ADF test showed that Bekok station, Kahang station, and Johor station had significant non-stationary streamflow. Meanwhile, Linggui station presented stationary streamflow from the test. This shows the exact relationship between the change point analysis and non-stationary analysis, where a significant change point streamflow indicated a non-stationary result, and vice versa. However, it was different for Johor station, which produced conflicting results between the two tests.

As stated in the previous section, the non-stationary models with different values of shape parameter (k) and mean slope (ξ₁) at both ranges from −0.3 to 0.3 were used in the five methods (L-moments and TL-moments, TL(η,0), (η = 1, 2, 3, 4)). The standard deviation (α) and ξ₀ were fixed to 1 and 0, respectively. In this study, the simulation result with the selected sample size (n = 50 and 100) is presented in Figure 3.

The boxplot is used as an illustration tool in order to evaluate the strength of the L-moment and TL-moments method. The criteria for selecting the most suitable methods are based on the minimum value of the median level where the ratio of the estimated and real quantiles is close to 1 or both ends of the boxplot. Figure 3 displays the boxplots of the RRMSE result using the methods considered (L-moment and TL-moments method) for (a) GEV1 model, (b) GEV2 model, and (c) GEV3 model. Figure 3(a) shows that all boxplots illustrating the median value of RRMSE of the GEV1 model are in the range (1.1,1.0) which show minor differences between L-moment and TL-moments method. However, the TL-moments method is more robust (close to 1) than the L-moment method for the GEV1 model.

The boxplot for the GEV2 model (Figure 3(b)) shows that the minimum dispersion in the median RRMSE values indicated by both ends of the boxplot is also the criterion for selecting the most suitable method. The L-moment and TL(1,0) method give a better performance that is shown by minimum dispersion compared to other TL-moments, TL(η,0), (η = 2, 3, 4) method. The GEV3 model (Figure 3(c)) at positive k and negative ξ₂ presents the TL(η,0), (η = 2, 3, 4) method, which is more robust than the L-moment method based on the smallest interquartile series of the boxplots.

Thus, TL-moments method tends to a positive trend parameter in mean and standard deviation for the GEV1, GEV2, and GEV3 models in a simulation study. The TL-moments method can produce better estimation results compared to L-moments since the former eliminates some of the lowest values and gives more weight to the largest value, which usually carries important information.

The bootstrap resampling technique is applied to choose the suitable model in the prediction of flood for 200 years and 1,000 years by seeing the performances from the comparison methods, L-moments and TL-moments methods. Table 5 presents 0.995 and 0.999 return periods of minimum value of the RRMSE results in the flood series calculated by L-moments and TL-moments for three selected models: GEV1, GEV2, and GEV3 models. The smaller value of the RRMSE results calculated by selected method are in bold to highlight the result in Table 5.

For station 2130422 Bekok and 1737451 Johor, TL-moments, TL(η,0), (η = 1, 2, 3, 4) methods produce small values of RRMSE results in all selected models, GEV1, GEV2, and GEV3 for all quantiles. However, smaller RRMSE values are also obtained by stationary and non-stationary L-moments method in both GEV2 and GEV3 models where each model represents station 2235401 Kahang and 1836401 Linggui. Thus, these findings suggest a role for the proposed TL-moments method as the estimating method in the non-stationary parameter estimation field.

From the simulation and bootstrap study, GEV1, GEV2, and GEV3 models show good results in simulation and real data application when estimated by TL-moments method. Thus, the TL-moments method gives zero weight to extreme values and makes it easier to compute and is more robust than the L-moments method.

Table 6 shows the results of the BIC test of the GEV0, GEV1, GEV2, and GEV3 models for the four selected stations. The best values of BIC values are highlighted in the table. The test is used to find the best model for the non-stationary flood estimation.

For station 2130422 Bekok, this flood series is shown as significant non-stationary streamflow from the Pettitt test and ADF test (Table 4). The GEV3 model is found to be the best model for Bekok station based on the BIC test. This outcome is also highlighted by the results of the quantile estimated (Table 5) that illustrates the TL-moments method (i.e., TL(1,0) and TL(4,0)) as the chosen estimated method in the non-stationary flood series of 2130422 Bekok station for the non-stationary model.

Stations 2235401 Kahang and 1737451 Johor also presented as non-stationary flood peak stations, which indicated the GEV3 and GEV2 models were the best fit models for the data. From the quantile with estimated RRMSE result (Table 5), the 2235401 Kahang station chooses the L-moments method as the estimated method for the GEV3 model. Otherwise, the TL-moments method is selected by station 1737451 Johor for all non-stationary models.

However, there are dramatic differences in the results found for station 1836401 Linggui where the non-stationary GEV model (GEV2) is found to be the best model for this station, while it presented as the stationary data flood series based on the ADF test. This station chooses the L-moments method for estimate quantiles of GEV2 model as presented in Table 5.

Table 7 presents the model parameters of the performance criteria for all stations. The estimated model parameters are chosen based on the results of model selection (Table 6). From the table, Kahang and Johor station show a high value of location and scale parameters. The scale parameter refers to the characteristic of flood risk for a station (known as the gradient of the extreme discharge) (Hounkpe et al. 2015). Kahang station shows the highest value for this parameter. The shape parameter links to a greater extreme discharge for a greater absolute value of the shape parameter. This refers to Kahang station, which also showed the presence of the change point for the streamflow (Figure 2).

According to the location parameter, μ₁ obtained from Table 7, the values are positive, and estimated mostly from TL-moments method. Similarly, the TL-moments method is the best in the case of positive trends in the mean, shown from the simulation study. Thus, this shows a connection between the selection parameter in the simulation study and the real data of streamflow stations.

The purpose of the current study is to decide upon the stationarity in the annual maximum flow series and to investigate the best model of four selected stations in Johor estimated by the comparison method between L-moments method and the proposed method (TL-moments method). The estimated parameters of the GEV distribution in the context of stationary (GEV0) and non-stationary were applied to the original flood data series. Three non-stationary models presented a linear trend of the mean (GEV1), the linear trend in mean and standard deviation (GEV2), and parabolic trend in mean (GEV3).

The Pettitt test and ADF test were used to check the abrupt changes and stationarity in the flood series. Two stations were identified with significant change point streamflow even when all stations detected abrupt changes in the mean of the AMS series. However, three stations were identified as the non-stationary streamflow series where the non-stationarity of the stations in Johor was due to the presence of the change point in the streamflow series. It should also be noted that the detected change point may not be statistically robust enough due to the limited length of time series.

By determining the efficiency of the TL-moments approach for the non-stationary case, the ability of each method was carried out through Monte Carlo simulation. The boxplots applied to the better illustration of the RRMSE result. As well, the performance of the proposed method in selecting the best model for accurately predicting rare floods was studied through residual bootstrap.

From the simulation findings, the TL-moments method was more appropriate for estimating the parameters of the non-stationary GEV model than the L-moments method for GEV1 and GEV3 models. However, the TL-moments method competed closely with the L-moments method for cases in the GEV2 model, in which L-moment and TL(1,0) methods were found to be the most accurate methods in this model. Thus, TL-moments method can produce better estimation results compared to L-moments; TL-moments can assign zero weight to extreme observations (i.e., easy to compute) and more robust than L-moments method.

By using the residual bootstrap and likelihood test, the performances of flood for 200 years and 1,000 years mostly observed similar results from these studies. This can be clearly seen by the non-stationary stations (Bekok, Kahang, and Johor) where the good results obtained from their selected model in the residual bootstrap had shown a good fit to the similar model based on the BIC test. However, significant differences from both results was shown in this study too. This was seen in Linggui station, which had a good fit with the GEV2 model based on the BIC test. Nevertheless, the TL-moments method estimated poor results of RRMSE in the GEV2 model according to the residual bootstrap for this station.

The research also showed the parameters of the selected model estimated from the successful methods in the simulation study. The location parameter estimated showed the same trend parameter in the standard deviation as the parameters used in the simulation. Therefore, it is likely that such connections exist between simulation parameters used with the actual data.

Also, the proposed method, namely, the TL-moments method, can be used to estimate the non-stationary parameter model by improving the quality of flood estimation in the research of flood frequency analysis in a non-stationary context. The additional parameter in the non-stationary model for extreme flood frequency analyses by linking it with time also plays an important role in this study. Also, the connection between non-stationarity and climate change should be taken into account in future studies.

Special thanks to the Ministry of Higher Education (MOHE) and Ministry of Science, Technology and Innovation (MOSTI) Malaysia for funding this research under Vot 4F681, 18H55, and 4F875. We thank the Department of Irrigation and Drainage, Ministry of Natural Resources and Environment, Malaysia for providing the data for this study. Last but not least, we would like to express our utmost gratitude to Universiti Teknologi Malaysia.