Abstract

Non-stationary flood frequency analysis (NFFA) plays an important role in addressing the issue of the stationary assumption (independent and identically distributed flood series) that is no longer valid in infrastructure-designed methods. This confirms the necessity of developing new statistical models in order to identify the change of probability functions over time and obtain a consistent flood estimation method in NFFA. The method of Trimmed L-moments (TL-moments) with time covariate is confronted with the L-moment method for the stationary and non-stationary generalized extreme value (GEV) models. The aims of the study are to investigate the behavior of the proposed TL-moments method in the presence of NFFA and applying the method along with GEV distribution. Comparisons of the methods are made by Monte Carlo simulations and bootstrap-based method. The simulation study showed the better performance of most levels of TL-moments method, which is TL(η,0), (η = 2, 3, 4) than the L-moment method for all models (GEV1, GEV2, and GEV3). The TL-moment method provides more efficient quantile estimates than other methods in flood quantiles estimated at higher return periods. Thus, the TL-moments method can produce better estimation results since the L-moment eliminates lowest value and gives more weight to the largest value which provides important information.

INTRODUCTION

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

HYDROLOGICAL DATA

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.

Table 1

Characteristics of the selected stations

River nameArea (km²)Number of observationsLatitudeLongitude
2130422 Bekok 350 33 02 07 15 103 02 30 
2235401 Kahang 587 32 02 15 05 103 35 15 
1737451 Johor 1,130 45 01 46 50 103 44 45 
1836401 Linggui 209 15 01 53 45 103 41 30 
River nameArea (km²)Number of observationsLatitudeLongitude
2130422 Bekok 350 33 02 07 15 103 02 30 
2235401 Kahang 587 32 02 15 05 103 35 15 
1737451 Johor 1,130 45 01 46 50 103 44 45 
1836401 Linggui 209 15 01 53 45 103 41 30 

METHODOLOGY

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.

Figure 1

Overview of the study.

Figure 1

Overview of the study.

Change point and non-stationary analysis

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 Pettit 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 H0, 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 (X1, X2, … XT) exist with the change point at time T when Xt (for t = 1, 2 … T) has the same distribution function F1(x), and Xt (for t = τ + 1 … T) has the same distribution function F2(x), where F1(x) ≠ F2(x) (Jha & Sigh 2013).

The non-parametric statistic is defined as: 
formula
(1)
where 
formula
(2)

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

The confidence level associated with KT is determined approximately by: 
formula
(3)

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.

The non-stationary analysis is used in order to calculate whether or not a series of mean and variance values vary with time (Van Gelder et al. 2007). A statistical test is adopted to investigate non-stationarity in extreme streamflow in time series data. The Augmented Dickey–Fuller (ADF) test, introduced by Dickey & Fuller (1979) is applied to test the presence of unit roots in the series (difference stationarity). This test is conducted through the ordinary least square estimation of regression models incorporating either an intercept or a linear trend. The AR(1) model is considered as: 
formula
(4)
where x0 = 0, |ρ| ≤ 1, and ɛt is a real-valued sequence of independent random variables with mean zero and variance. If ρ = 1, the process {xt} is non-stationary, where it is known as a random walk process. In contrast, the process {xt} is stationary if |ρ| < 1. The MLE of ρ is the least squares estimator: 
formula
(5)
The null hypothesis of the ADF test is non-stationarity of the time series data. Whenever the p-value (probability) of the test statistics is lower than 0.05 confidence level, the null hypothesis is rejected.

Non-stationary GEV model candidates and model selection

The GEV distribution is a common distribution that is often applied in NFFA (Gado & Nguyen 2016a, 2016b). The cumulative distribution (cdf), f(x) and probability distribution (pdf), F(x) of GEV distribution are: 
formula
(6)
 
formula
(7)
where the parameters of the standard GEV distribution are location (ξ), scale (α), and shape (k) parameter.

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.

Table 2

Parameters for the GEV non-stationarity

ModelLocation parameter, ξScale parameter, αShape parameter, k
GEV1 ξ(t) = ξ0 + ξ1t α = constant k = constant 
GEV2 ξ(t) = ξ0 + ξ1t α(t) = exp(α0 + α1tk = constant 
GEV3 ξ(t) = ξ0 + ξ1t+ξ2t2 α = constant k = constant 
ModelLocation parameter, ξScale parameter, αShape parameter, k
GEV1 ξ(t) = ξ0 + ξ1t α = constant k = constant 
GEV2 ξ(t) = ξ0 + ξ1t α(t) = exp(α0 + α1tk = constant 
GEV3 ξ(t) = ξ0 + ξ1t+ξ2t2 α = constant k = constant 

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.

Parameter estimation method

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

The corresponding PWM equation for the GEV distribution is given by Equation (8) (Hosking et al. 1985): 
formula
(8)
where ξ is a location parameter, α is a scale or lower bound parameter, and k is the shape parameter of the GEV distribution. The PWM relationships (r=h, h+ 1, and h+ 2) of Equation (8) are shown in the following equations: 
formula
(9)
 
formula
(10)
 
formula
(11)

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:

  • 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: 
    formula
    (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: 
    formula
    (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.

Simulation and bootstrap method study

Table 3

Parameter estimation of GEV1, GEV2, and GEV3 model by the L-moments and TL-moments (η= 1, 2, 3, 4) methods

(a)
Method
L-momentsTL-moments (1,0)TL-moments (2,0)
Moments Mean    
Standard deviation    
Parameter estimation GEV1 model    
   
   
   
GEV2 model    
   
   
   
GEV3 model    
   
   
   
(b) 
Method   TL-moments (3,0) TL-moments (4,0) 
Moments Mean    
Standard deviation    
Parameter estimation GEV1 model    
   
   
GEV2 model    
   
   
GEV3 model    
   
   
(a)
Method
L-momentsTL-moments (1,0)TL-moments (2,0)
Moments Mean    
Standard deviation    
Parameter estimation GEV1 model    
   
   
   
GEV2 model    
   
   
   
GEV3 model    
   
   
   
(b) 
Method   TL-moments (3,0) TL-moments (4,0) 
Moments Mean    
Standard deviation    
Parameter estimation GEV1 model    
   
   
GEV2 model    
   
   
GEV3 model    
   
   

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 (ξ1) 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 = 0 for the GEV1 model. For the GEV2 model, ξ0 = 0 and ξ1 = −0.1 for time-dependent in mean. Time-dependent in standard deviations, α0 = 0 and two cases of α1 = ±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, ξ1 = 0.3, and ξ2 = −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.

The criteria for the selection of the suitable method and the model for each case would be based on the minimum calculated RRMSE when simulated, and calculated quantiles are compared for the 0.95–0.999 quantiles. The RRMSE is represented by: 
formula
(14)

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.

RESULT AND ANALYSIS

Change point and non-stationary analysis

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.

Table 4

Results of the Pettitt test and ADF test for the annual maximum streamflows

River namePettitt test
ADF test
tau/yearKT1-p
2130422 Bekok 19/1989 320 0.9999 Non-stationary 
2235401 Kahang 18/1995 162 0.9811 Non-stationary 
1737451 Johor 13/1977 148 0.5121 Non-stationary 
1836401 Linggui 6/1983 26 0.3518 Stationary 
River namePettitt test
ADF test
tau/yearKT1-p
2130422 Bekok 19/1989 320 0.9999 Non-stationary 
2235401 Kahang 18/1995 162 0.9811 Non-stationary 
1737451 Johor 13/1977 148 0.5121 Non-stationary 
1836401 Linggui 6/1983 26 0.3518 Stationary 

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.

Figure 2

The change point result from the Pettitt test for annual maximum streamflow of Bekok station and Kahang station.

Figure 2

The change point result from the Pettitt test for annual maximum streamflow of Bekok station and Kahang station.

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.

Simulation and bootstrap study

As stated in the previous section, the non-stationary models with different values of shape parameter (k) and mean slope (ξ1) 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 ξ0 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.

Figure 3

Boxplots of RRMSE values for (a) GEV1 model (n = 50, ξ0 = 0, ξ1 = −0.2, α = 1, and k = 0.2), (b) GEV2 model (n = 100, ξ0 = 0, ξ1 = −0.1, α0 = 1, and α1 = 0.02), and (c) GEV3 model (n = 50, ξ0 = 0, ξ1 = 0.3, ξ2 = −0.005, α = 1, and k = 0.2) with L-moments and TL-moments method estimated.

Figure 3

Boxplots of RRMSE values for (a) GEV1 model (n = 50, ξ0 = 0, ξ1 = −0.2, α = 1, and k = 0.2), (b) GEV2 model (n = 100, ξ0 = 0, ξ1 = −0.1, α0 = 1, and α1 = 0.02), and (c) GEV3 model (n = 50, ξ0 = 0, ξ1 = 0.3, ξ2 = −0.005, α = 1, and k = 0.2) with L-moments and TL-moments method estimated.

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 ξ2 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.

Table 5

The quantile estimated results of RRMSE values for GEV1, GEV2, and GEV3 models at 0.995 and 0.999 return periods estimated by L-moments and TL-moments method based on the residual bootstrap

No.River nameMethodGEV1
GEV2
GEV3
0.9950.9990.9950.9990.9950.999
2130422 Bekok Stationary 0.670 0.600 0.543 0.487 0.731 0.675 
L-mom 0.766 0.685 0.637 0.571 0.748 0.691 
TL(1,0) 0.557 0.498 0.342 0.306 0.506 0.467 
TL(2,0) 0.545 0.487 0.277 0.248 0.464 0.429 
TL(3,0) 0.567 0.507 0.234 0.209 0.476 0.440 
TL(4,0) 0.570 0.510 0.209 0.187 0.506 0.467 
2235401 Kahang Stationary 0.767 0.714 0.224 0.208 0.124 0.117 
L-mom 0.763 0.711 0.347 0.322 0.056 0.053 
TL(1,0) 0.632 0.589 0.209 0.194 0.065 0.062 
TL(2,0) 0.546 0.508 0.088 0.082 0.066 0.062 
TL(3,0) 0.532 0.495 0.102 0.095 0.094 0.088 
TL(4,0) 0.596 0.555 0.145 0.135 0.135 0.127 
1737451 Johor Stationary 0.649 0.600 0.748 0.709 0.347 0.327 
L-mom 0.606 0.561 0.571 0.541 0.332 0.312 
TL(1,0) 0.543 0.503 0.291 0.276 0.291 0.274 
TL(2,0) 0.520 0.481 0.119 0.113 0.278 0.262 
TL(3,0) 0.501 0.463 0.051 0.049 0.280 0.264 
TL(4,0) 0.486 0.449 0.142 0.134 0.291 0.274 
1836401 Linggui Stationary 0.085 0.084 0.174 0.152 0.087 0.087 
L-mom 0.017 0.017 1.016 0.887 0.079 0.079 
TL(1,0) 0.213 0.212 1.024 0.894 0.057 0.057 
TL(2,0) 0.256 0.256 0.935 0.816 0.040 0.040 
TL(3,0) 0.149 0.149 1.073 0.937 0.161 0.161 
TL(4,0) 0.009 0.009 1.042 0.910 0.071 0.071 
No.River nameMethodGEV1
GEV2
GEV3
0.9950.9990.9950.9990.9950.999
2130422 Bekok Stationary 0.670 0.600 0.543 0.487 0.731 0.675 
L-mom 0.766 0.685 0.637 0.571 0.748 0.691 
TL(1,0) 0.557 0.498 0.342 0.306 0.506 0.467 
TL(2,0) 0.545 0.487 0.277 0.248 0.464 0.429 
TL(3,0) 0.567 0.507 0.234 0.209 0.476 0.440 
TL(4,0) 0.570 0.510 0.209 0.187 0.506 0.467 
2235401 Kahang Stationary 0.767 0.714 0.224 0.208 0.124 0.117 
L-mom 0.763 0.711 0.347 0.322 0.056 0.053 
TL(1,0) 0.632 0.589 0.209 0.194 0.065 0.062 
TL(2,0) 0.546 0.508 0.088 0.082 0.066 0.062 
TL(3,0) 0.532 0.495 0.102 0.095 0.094 0.088 
TL(4,0) 0.596 0.555 0.145 0.135 0.135 0.127 
1737451 Johor Stationary 0.649 0.600 0.748 0.709 0.347 0.327 
L-mom 0.606 0.561 0.571 0.541 0.332 0.312 
TL(1,0) 0.543 0.503 0.291 0.276 0.291 0.274 
TL(2,0) 0.520 0.481 0.119 0.113 0.278 0.262 
TL(3,0) 0.501 0.463 0.051 0.049 0.280 0.264 
TL(4,0) 0.486 0.449 0.142 0.134 0.291 0.274 
1836401 Linggui Stationary 0.085 0.084 0.174 0.152 0.087 0.087 
L-mom 0.017 0.017 1.016 0.887 0.079 0.079 
TL(1,0) 0.213 0.212 1.024 0.894 0.057 0.057 
TL(2,0) 0.256 0.256 0.935 0.816 0.040 0.040 
TL(3,0) 0.149 0.149 1.073 0.937 0.161 0.161 
TL(4,0) 0.009 0.009 1.042 0.910 0.071 0.071 

Values in bold indicate that the smaller value of RRMSE.

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.

Model selection

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.

Table 6

The BIC test of the GEV0, GEV1, GEV2, and GEV3 models for the four selected stations

River nameBIC test
GEV0GEV1GEV2GEV3
2130422 Bekok − 339.1416 − 321.594 − 321.591 −320.831 
2235401 Kahang − 410.4218 − 405.384 − 404.168 −402.436 
1737451 Johor − 551.6774 − 551.412 −550.31 − 545.534 
1836401 Linggui − 138.1357 − 131.579 −57.5333 − 105.163 
River nameBIC test
GEV0GEV1GEV2GEV3
2130422 Bekok − 339.1416 − 321.594 − 321.591 −320.831 
2235401 Kahang − 410.4218 − 405.384 − 404.168 −402.436 
1737451 Johor − 551.6774 − 551.412 −550.31 − 545.534 
1836401 Linggui − 138.1357 − 131.579 −57.5333 − 105.163 

Values in bold indicate the smaller value of BIC test for each station.

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

Table 7

Parameters of the GEV model

StationModelMethodLocation, ξ
Scale, α
Shape, k
μ0μ1μ2α0α1
Bekok GEV3 TL(4,0) 40.029 0.536 −0.020 12.653 – −0.150 
Kahang GEV3 L-mom 216.046 11.693 −0.282 129.132 – −0.218 
Johor GEV2 TL(3,0) 182.354 1.492 – 0.154 147.007 −0.038 
Linggui GEV2 TL(4,0) 19.975 2.368 – 3.960 −9.629 −0.101 
StationModelMethodLocation, ξ
Scale, α
Shape, k
μ0μ1μ2α0α1
Bekok GEV3 TL(4,0) 40.029 0.536 −0.020 12.653 – −0.150 
Kahang GEV3 L-mom 216.046 11.693 −0.282 129.132 – −0.218 
Johor GEV2 TL(3,0) 182.354 1.492 – 0.154 147.007 −0.038 
Linggui GEV2 TL(4,0) 19.975 2.368 – 3.960 −9.629 −0.101 

According to the location parameter, μ1 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.

CONCLUSION

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.

ACKNOWLEDGEMENTS

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.

REFERENCES

REFERENCES
Dickey
D. A.
&
Fuller
W. A.
1979
Distribution of the estimators for autoregressive time series with a unit root
.
Journal of the American Statistical Association
74
(
366
),
423
431
.
doi:10.2307/2286348
.
Elamir
E. A. H.
&
Seheult
A. H.
2003
Trimmed L-moments
.
Computational Statistics & Data Analysis
43
(
3
),
299
314
.
doi:10.1016/S0167-9473(02)00250-5
.
Gado
T. A.
&
Nguyen
V.-T.
2016a
An at-site flood estimation method in the context of non-stationarity I. A simulation study
.
Journal of Hydrology
535
,
710
721
.
doi:10.1016/j.jhydrol.2015.12.063
.
Gado
T. A.
&
Nguyen
V.-T.
2016b
An at-site flood estimation method in the context of non-stationarity II. Statistical analysis of floods in Quebec
.
Journal of Hydrology
535
,
722
736
.
doi:10.1016/j.jhydrol.2015.12.064
.
Hasan
H.
,
Salam
N.
&
Adam
M. B.
2013
Modelling extreme temperature in Malaysia using Generalized Extreme Value distribution
.
International Journal of Mathematical and Computational Sciences
7
(
6
),
983
989
.
doi:10.1063/1.4801267
.
Hosking
J. R. M.
1990
L-Moments: analysis and estimation of distributions using linear combinations of order statistics
.
Journal of Royal Statistical Society B
52
,
105
124
. .
Hounkpè
J.
,
Diekkrüger
B.
,
Badou
D.
&
Afouda
A.
2015
Non-stationary flood frequency analysis in the Ouémé River Basin
.
Benin Republic. Hydrology
2
(
4
),
210
229
.
doi:10.3390/hydrology2040210
.
Ismail
A.
2018
Perubahan iklim, pemanasan global (Climate change, global warming), Utusan Online. Available from: http://www.utusan.com.my/sains-teknologi/sains/perubahan-iklim-pemanasan-global-1.747481 (accessed 31 January 2019)
.
Jha
M. K.
&
Singh
A. K.
2013
Trend analysis of extreme runoff events in major river basins of Peninsular Malaysia
.
International Journal of Water
7
(
1
),
142
158
.
doi:10.1504/IJW.2013.051995
.
Khaliq
M. N.
,
Ouarda
T. B. M. J.
,
Ondo
J. C.
,
Gachon
P.
&
Bobée
B.
2006
Frequency analysis of a sequence of dependent and/or non-stationary hydro-meteorological observations: a review
.
Journal of Hydrology
329
(
3–4
),
534
552
.
doi:10.1016/j.jhydrol.2006.03.004
.
Khaliq
M. N.
,
Ouarda
T. B. M. J.
,
Gachon
P.
,
Sushama
L.
&
St-Hilaire
A.
2009
Identification of hydrological trends in the presence of serial and cross correlations: a review of selected methods and their application to annual flow regimes of Canadian rivers
.
Journal of Hydrology
368
(
1–4
),
117
130
.
doi:10.1016/j. jhydrol.2009.01.035
.
Li
G.
,
Xiang
X.
&
Guo
C.
2016
Analysis of non-stationary change of annual maximum level records in the Yangtze River Estuary
.
Advances in Meteorology
2016
,
1
14
.
doi:10.1155/2016/7205723
.
Mat Jan
N. A.
&
Shabri
A.
2017
Estimating distribution parameters of annual maximum streamflows in Johor, Malaysia using TL-moments approach
.
Theoretical and Applied Climatology
127
,
213
227
.
doi:10.1007/s00704-015-1623-7
.
Mat Jan
N. A.
,
Shabri
A.
,
Hounkpè
J.
&
Badyalina
B.
2018
Modelling non-stationary extreme streamflow in Peninsular Malaysia
.
International Journal of Water
12
(
2
),
116
140
.
doi:10.1504/IJW.2018.10012408
.
NIST/SEMATECH
2012
e-Handbook of Statistical Methods
.
NIST/SEMATECH
.
Available from: http://www.itl.nist.gov/div898/handbook/ (accessed 12 March 2015)
.
Noura
A. T. A. E.
2010
TL-moments of the exponentiated generalized extreme value distribution
.
Journal of Advanced Research
1
(
4
),
351
359
.
doi:10.1016/j.jare.2010.06.003
.
Ouarda
T. B. M. J.
&
El-Adlouni
S.
2011
Bayesian non-stationary frequency analysis of hydrological variables
.
Journal of the American Water Resources Association
47
(
3
),
496
505
.
doi:10.1111/j.1752-1688.2011.00544.x
.
Pettitt
A. N.
1979
A non-parametric approach to the change point problem
.
Applied Statistics
28
(
2
),
126
135
.
doi:10.2307/2346729
.
Schwarz
G.
1978
Estimating the dimension of a model
.
Annals of Statistics
6
(
2
),
461
464
. .
Shahzad
M. N.
&
Ashgar
Z.
2013
Parameter estimation of Sigh Maddala distribution by moments
.
International Journal of Advanced Statistics and Probability
1
(
3
),
121
131
.
doi:10.14419/ijasp.v1i3.1206
.
Sober
E.
2002
Instrumentalism, parsimony, and the Akaike framework
.
Philosophy of Science
69
,
S112
S123
.
doi:10.1086/341839
.
Strupczewski
W. G.
,
Kochanek
K.
,
Bogdanowicz
E.
,
Markiewicz
I.
&
Feluch
W.
2016
Comparison of two nonstationary flood frequency analysis methods within the context of the variable regime in the representative polish rivers
.
Acta Geophysica
64
(
1
),
206
236
.
doi:10.1515/acgeo-2015-0070
.
Tan
X.
&
Gan
T.
2015
Non-stationary analysis of annual maximum streamflow of Canada
.
Journal of Climate
28
(
5
),
1788
1805
.
doi:10.1175/JCLI-D-14-00538.1
.
Trenberth
K. E.
,
Jones
P. D.
,
Ambenje
P.
,
Bojariu
R.
,
Easterling
D.
,
Klein Tank
A.
,
Parker
D.
,
Rahimzadeh
F.
,
Renwick
J. A.
,
Rusticucci
M.
,
Soden
B.
,
Zhai
P.
2007
Observations: Surface and atmospheric climate change
. In:
IPCC Fourth Assessment Report: Climate Change 2007. Working Group I: The Physical Science Basis
(
Solomon
S.
,
Qin
D.
,
Manning
M.
,
Chen
Z.
,
Marquis
M.
,
Averyt
K. B.
,
Tignor
M.
&
Miller
H. L.
, eds).
Cambridge University Press
,
Cambridge
,
UK
, pp.
235
336
. .
Van Gelder
P. H. A. J. M.
,
Wang
W.
,
Vrijling
J. K.
2007
Statistical estimation methods for extreme hydrological events
. In:
Extreme Hydrological Events: New Concepts for Security. NATO Science Series
(
Vasiliev
O. F.
,
van Gelder
P. H. A. J. M.
,
Plate
E. J.
&
Bolgov
M. V.
, eds).
Springer
,
Dordrecht
,
The Netherlands
, pp.
199
252
.
doi:10.1007/978-1-4020-5741-0_15
.
Villarini
G.
,
Serinaldi
F.
,
Smith
J. A.
&
Krajewski
W. F.
2009
On the stationarity of annual flood peaks in the continental United States during the 20th century
.
Water Resources Research
45
(
8
),
W08417
.
doi:10.1029/2008WR007645
.
Villarini
G.
,
Smith
J. A.
&
Napolitano
F.
2010
Non-stationary modelling of a long record of rainfall and temperature over Rome
.
Advances in Water Resources
33
(
10
),
1256
1267
.
doi:10.1016/j.advwatres.2010.03.013
.
Wang
Q. J.
1990
Estimation of the GEV distribution from censored samples by method of partial probability weighted moments
.
Journal of Hydrology
120
(
1–4
),
103
114
.
doi:10.1016/0022-1694(90)90144-M
.
Yilmaz
A. G.
&
Perera
B. J. C.
2013
Extreme rainfall non-stationarity investigation and Intensity–Frequency–Duration relationship
.
Journal of Hydrologic Engineering
19
(
6
),
1160
1172
.
doi:10.1061/(ASCE)HE.1943-5584.0000878
.