Estimating the design flood under the influence of check dams by removing nonstationarity from the flood peak discharge series

Design flood estimation plays an important role in water project design, water resources planning, and flood risk control. It is usually obtained using frequency analysis which is a technique of fitting a probability distribution to a series of observations for defining the probabilities of future occurrences of some events of interest, e.g., an estimate of a flood magnitude corresponding to a chosen risk of failure (Khaliq et al. 2006). The conventional hydrologic frequency analysis method requires the hydrological series to be independent and identically distributed (Salas & Obeysekera 2014; Qin et al. 2015; Read & Vogel 2015). For the flood peak discharge Q with distribution ⁠, we generally assume that the series is stationary, and that the distribution parameters are time invariant. In northwestern China, climate and topographical conditions cause serious water and soil loss. As such, many check dams have been built to try to mitigate this phenomenon. In the context of the construction of thousands of check dams, the stationarity of the flood peak discharge series gradually disappears. The nonstationary flood frequency analysis has become a point of concern for engineers and local managers.

Under nonstationary conditions, in the absence of understanding the physical mechanisms, we usually assume that the distribution parameters of the hydrological series change with time; the time-varying distribution is determined by fitting the time-varying characteristic of the distribution parameters from a statistical perspective. Based on this logic, Strupczewski et al. (2001a, 2001b) and Strupczewski & Kaczmarek (2001) proposed the time-varying moments method. The Generalized Additive Models for Location, Scale, and Shape (GAMLSS), which was proposed by Rigby & Stasinopoulos (2005), makes the application of time-varying moments method more flexible. Usually, the model uses time as the covariate so it can describe the trend variation in the hydrological time series; however, this is not sufficient to reflect the driving force of the nonstationary change. Therefore, a number of studies have analyzed the influence of climate change and anthropogenic activity on the nonstationary hydrological time series by using them as the covariates (Villarini et al. 2009; López & Francés 2013; Xiong et al. 2014; Zhang et al. 2015). Using the traditional definition of the return period, where the return period is equal to the reciprocal of the exceedance probability, results in unrealistic design values that are not appropriate for real-world applications because the design value varies from one year to the next (Yan et al. 2017).

In order to cope with nonstationarity, some studies have extended the concept of return period to a more general interpretation as the expected waiting time (EWT), which is the expected time interval between two exceedance events (Wigley 1988; Olsen et al. 1998). Alternatively, other works have employed the return period as the time period during which the expected number of exceedance (ENE) events is equal to 1 (Parey et al. 2007, 2010). Cooley (2013) and Yan et al. (2017) pointed out that the EWT definition requires indefinite extrapolation along the trend, while the ENE definition only requires extrapolation for a finite number of years. Based on the EWT or ENE definitions, we can obtain the unique design value for a nonstationary series with a given return period. However, under nonstationary conditions, different initial and end time of calculation period results in different hydrological design values; the EWT and ENE methods do not have a clear definition of the calculation period. In practice, many studies take the project's design lifetime as the calculation period (Rootzén & Katz 2013; Yan et al. 2017).

Some studies illustrate that the return period does not reflect the risk of failure of a project, and recommend to replace it with the reliability, which is the probability of no exceedance events occurring over the lifetime of the project (Rootzén & Katz 2013; Salas & Obeysekera 2014; Read & Vogel 2015). Reliability is a more effective index than the return period because it summarizes the joint probability instead of the average probability (Serinaldi & Kilsby 2015). On this basis, some reliability-based methods such as the ‘design life level’ (Rootzén & Katz 2013), ‘equivalent reliability’ (Hu et al. 2017), and ‘average design life level’ (Read & Vogel 2015) have been proposed successively.

Whether the return period or the reliability, we usually construct its expression according to the time-varying probability distribution, and then use it to calculate the design value. However, it is impossible to accurately describe the characteristics of hydrological time series or to predict the influence degree of climate change or human activities in the future simply by pure statistical approaches unless you have a clear understanding of the physical mechanism. Even if meteorological or anthropogenic factors are chosen as the covariates, the nonstationary changes are usually characterized using data-based regression technique. This nonstationary change is usually fitted with linear or some simple nonlinear curves, but the actual conditions must be extremely complicated, and it is difficult to answer how this nonstationary change will evolve in the future. Besides, according to statistical theory, because the nonstationary hydrological series no longer has ergodicity, the population distribution cannot be inferred from it.

For a long time, we have tried to reformulate and extend the conventional hydrological frequency analysis method to adapt the nonstationary conditions. Since there is insufficient knowledge of the physical mechanism, for nonstationary frequency analysis, new challenges and doubts are encountered as one problem after another is solved. However, under the stationary conditions, the conventional frequency analysis method has been widely proved to be scientific and reasonable. The conventional method has been mastered by engineers, and there are a large number of guides and criteria to ensure its reliability. Therefore, the conventional frequency analysis method should not be discarded easily. If the stationarity of a hydrological series can be reconstructed, then the conventional method can still be used to calculate the design value.

For a basin that has a high number of check dams, we analyze the influence of the check dams on the nonstationary changes in the flood peak discharge series according to its construction and operation characteristics. According to this influence law, we try to remove the nonstationary influence of the check dams from the flood peak discharge series and transform it into a new stationary series. The de-nonstationarity idea is not new in hydrology. Previous studies are based on the theory of stochastic hydrology to divide the hydrological series into a deterministic component (trend, periodicity, etc.) and a stochastic component (Kottegoda 1980; Xie et al. 2005). After removing the deterministic component, the stationary stochastic component can be obtained. However, this simplification is not accurate enough for the complex changes of hydrological series. The mechanism-based de-nonstationarity method proposed in this paper does not consider the type of nonstationarity, no matter what type of nonstationary changes of hydrological variable is due to the corresponding nonstationary change of its influence factors. The de-nonstationarity process is, indeed, the process of removing the nonstationary effects of these factors. Moreover, through theoretical derivation and practical application, we prove that the de-nonstationarity method proposed in this study can achieve stationarity in both the first and second moments of the flood peak discharge series. The resulting uncertainty of design flood, determined using the bootstrap method (Efron 1992), is small.

The check dams constructed in northwestern China have largely alleviated the water and soil loss. There are more than 20,000 check dams in Yulin in the Shaanxi Province alone (Jiang & Shuai 2011). This kind of high-density water project leads to significant nonstationary changes in the flood peak discharge series, which poses a great challenge to the water project construction and flood risk management. We take the Mahuyu River Basin (MRB) in Yulin City, Shaanxi Province, China as an example to calculate the design flood in the presence of check dams. The Mahuyu River is located in the hinterland of the Loess Plateau, sourced from Shiyaogou in Hengshan County, and flows from west to east into Wuding River. The MRB lies in the area bounded by 37°44′–37°55′N, 109°37′–110°02′E (Figure 1). The control station for the MRB is the Mahuyu hydrological station, which controls an area of 372 km². The MRB flood peak discharge series recorded at the Mahuyu station from 1962 to 2010 is shown in Figure 2. We employed the Mann–Kendall (M–K) test (Mann 1945; Kendall 1975) and the Breusch–Pagan (B–P) (Breusch & Pagan 1979) test to analyze the trends in the first and second moments, respectively (Table 1). Under the influence of check dams, both the first and second moments of the MRB flood peak discharge series exhibit nonstationary changes, where the nonstationarity in the second moment is more significant. Except for the trend, we do not consider other types of nonstationarity in this study.

A check dam is generally composed of two parts: the dam and the horizontal pipe drainage structure. The slow drainage (about 1.0–1.5 m³/s) of the horizontal pipes in the check dam causes the sediment in the flood flow to be deposited in front of the dam, which helps mitigate the soil loss. The runoff in the control area of the check dams cannot be discharged to the outlet of the basin in time. That is to say, for a flood event, the runoff generated in the control area of the check dam provides little contribution to the flood at the outlet of the basin. Therefore, we infer that the nonstationary changes in the flood peak discharge series are essentially caused by the change in the control area of the check dams. During the flooding process, areas in the basin that are not controlled by the check dams still generate runoff effectively, so this part of the basin area is defined as the effective runoff generation area. Figure 3 shows the interannual variation of the effective runoff generation area in the MRB. Dozens of check dams in the MRB were built year by year from 1962 to 2010, and the effective runoff generation area also changed gradually.

The presence of the check dams resulted in a reduced effective runoff generation area, which leads to nonstationary changes in the flood peak discharge series. After removing the influence of the effective runoff generation area from the original nonstationary series, the remaining part of the flood peak discharge series should be statistically stationary.

For hydrological frequency analysis, using sample series to estimate the population distribution will inevitably produce uncertainty. However, because the de-nonstationarity method combines the conventional frequency analysis with the deterministic influence of the check dams, the design flood will have a smaller uncertainty than those calculated using purely statistical approaches. On the other hand, because the nonstationary hydrological series no longer has ergodicity, we cannot actually estimate the population distribution according to the nonstationary sample series. The de-nonstationarity method transforms the nonstationary series into a stationary reconstructed series with ergodicity. However, with potential observation errors and a low sample size, it is still necessary to assess the design value uncertainty.

The bootstrap method proposed by Efron (1992) has been widely used in the quantitative analysis of uncertainty. Serinaldi & Kilsby (2015) pointed out that the bootstrap method strictly depends on the available information without any asymptotic assumptions, and that it does not rely on specific parameter estimation methods. Moreover, this method is relatively easy to implement, regardless of the complexity of the model. We use the bootstrap method to estimate the confidence intervals of the design flood. The procedure of the bootstrap method is as follows:

1. After removing the nonstationarity from the original flood peak discharge series Q(t), we apply the conventional frequency analysis method to stationary reconstructed series RS(t) and calculate the design value for a certain return period.

2. The design value of the original nonstationary flood peak discharge series is calculated with Equation (6).

3. We then resample the reconstructed series RS(t) with replacements to obtain a new sample series ⁠.

4. Based on the new sample series ⁠, we calculate the new design value for a specific return period.

5. Using Equation (6), we determine the design value of the nonstationary flood peak discharge series for a specific return period.

6. We then repeat the resampling process N times to obtain N design values: ⁠.

7. For a certain significance level α, the confidence interval of the design value is ⁠.

After removing the nonstationary influence of the effective runoff generation area, we obtained the stationary reconstructed series RS(t) (Figure 4). The hypothesis test results show that both the first and second moments of the reconstructed series of the MRB flood peak discharge are stationary (Table 2). Because the reconstructed series is stationary, we can use the conventional hydrological frequency analysis method to calculate the design value. The presence of the check dams drives changes in both the magnitude of the flood peak discharge and the rate of change relative to the mean, which reflects the nonlinear characteristics of the hydrological system. Because the statistics of B–P test decrease more significantly after removing the nonstationary effects of the effective runoff generation area from Q(t), we infer that the changes in the effective runoff generation area have a more significant impact on the second moment of Q(t) than they do on the first moment. Therefore, in the hydrological frequency analysis, we should account for nonstationarity in both the first moment and the second moment.

The effect of removing the nonstationary influence of the effective runoff generation area by formula (3) is similar to a kind of unitization, that is, the reconstructed series RS(t) is a rough approximation of the flood peak discharge per unit area; when the exponent in that expression is equal to 1, RS(t) represents the flood peak modulus series. The construction of check dams reduces the actual runoff generation area. However, in regions that are beyond the control area of the check dams, the runoff generation process for a flood event does not change, resulting in the flood peak modulus within the effective runoff generation area that is unaffected by the presence or absence of the check dams. Therefore, it is reasonable that the reconstructed series RS(t), which approximates the flood peak modulus, is stationary after the influence of the effective runoff generation area has been removed.

The nonstationary change of the flood peak discharge series and the effective runoff generation area in MRB can also be characterized as abrupt changes that occurred in the 1970s. The existing nonstationary hydrological frequency analysis methods, which follow the time-varying characteristics of the hydrological series, are always limited to a simple type of nonstationary change in fitting these nonstationary changes, but the actual conditions must be extremely complicated. In fact, nonstationary changes in a hydrological series are caused by the corresponding nonstationary changes in certain influencing factors. No matter what type of nonstationarity, it can be removed from the hydrological series using the de-nonstationarity method. To cope with nonstationarity, it is necessary to identify the nonstationary influence factors and their corresponding physical mechanisms.

In this paper, the frequency analysis of the reconstructed series of flood peak discharge (1962–2010) of Mahuyu hydrological station at the outlet of the MRB was carried out. Four distributions, Weibull (WEI), Gumble (GU), General extreme value (GEV), and Pearson type III (P-III), are selected as the alternative distributions in this study (Table 3).

In order to select the best-fitting distribution, we evaluated the overall simulation accuracy of the alternative distributions with the Kolmogorov–Smirnov (K–S) test (Kolmogorov-Smirnov et al. 1933), the Nash–Sutcliffe efficiency (⁠⁠) (Nash & Sutcliffe 1970) between the theoretical and empirical distributions (quantiles), and the root mean square error (RMSE). The goodness of fit of the four alternative distributions of the reconstructed MRB flood peak discharge series is shown in Table 4.

For stationary series RS(t), the K–S test statistic and show that the GEV distribution is optimal, while and RMSE indicate that the WEI distribution is a better fit. In general, the difference between the empirical probability and the theoretical probability for the same sample point at the tail of the distribution is small (both close to 1), but the difference between the empirical and theoretical quantile is relatively larger. Therefore, is more sensitive to the tail than ⁠. We conclude that the WEI distribution is the optimal distribution for the reconstructed series of the MRB flood peak discharge. The Q − Q plot in Figure 5 also shows the good performance of the optimal distribution.

We estimated the distribution parameters using the L-moments method (Hosking 1990). The L-moment is a linear combination of the probability weighted moments. Compared with the regular moment estimation method, the L-moments method is less affected by individual data points. After removing the influence of the effective runoff generation area from the flood peak discharge series, we estimated the distribution parameters of the reconstructed stationary series (Table 5) and calculated the design quantiles for a certain return period (Figure 6).

In order to obtain the design value for a certain return period of the original nonstationary flood peak discharge series, the design quantile for the same return period of the reconstructed series RS(t) should be multiplied by the influence of the effective runoff generation area under design conditions according to Equation (6). The changes in the effective runoff generation area mainly occurred in the 1970s. In recent years, the construction of check dams has been saturated, so the effective runoff generation area will not change drastically in the immediate future. By this logic, we can use the effective runoff generation area data from 2010 as the effective runoff generation area under the design conditions. According to Equation (6), the design quantiles for nonstationary flood peak discharge series are shown in Figure 7. The thin lines represent the 95% confidence intervals (CIs) for the design floods. As the return period increases, the CIs get wider, and the uncertainty in the design flood grows.

Because the existing nonstationary frequency analysis techniques are far from mature, most hydrological projects are still designed using the conventional hydrological frequency analysis method, regardless of whether the sample series is stationary. Due to the significant decreasing trend, if we calculate the design flood using the nonstationary MRB flood peak discharge series (the dashed line in Figure 7), the design flood for the same return period is larger than that of the value derived from the series we transformed with the de-nonstationarity method. The design values of the conventional method are nearly three times that of the de-nonstationarity method. The 100-year flood peak discharge is 1,829 m³/s, which is generated within the effective runoff generation area of only 56.5 km². This apparently larger, even somewhat impractical, design flood is unreliable. Because of the nonstationarity of the flood peak discharge series, it is obviously impossible to accurately assess the risk of flood by using conventional methods. Besides, the uncertainty of the design flood calculated by the original nonostationary flood peak discharge series is larger than that calculated by the de-nonstationary method. When the return period is more than 50 years, the width of uncertainty CIs of the design value by the conventional method has exceeded the design value itself. According to the physical mechanism of check dams on flood peak discharge, the de-nonstationarity method removes the nonstationarity and greatly reduces the uncertainty of the design value.

Some reliability-based nonstationary hydrological frequency analysis methods were introduced in the Introduction, including the expected number of exceedance (ENE), design life level (DLL), equivalent reliability (ER), and average design life level (ADLL). In this paper, we also use these four methods to calculate the design flood of MRB. With time and effective runoff generation area as the covariates, the GAMLSS model was used to model the nonstationary annual maximum flood peak discharge series. The alternative distributions were also WEI, GU, GEV, and P-III. To determine the best-fitting model (Table 6), the Akaike information criterion (AIC) (Akaike 1974) was employed.

Following the time-varying characteristics of the annual maximum flood peak discharge series, the four reliability-based nonstationary hydrologic frequency analysis methods were used to calculate the design value. It was assumed that the hydrologic structures of the MRB are planned to be in service from 1962 to 2061, with a lifespan of 100 years. Uncertainties were also calculated with the bootstrap method. The design values of different return periods are listed in Table 7 (other details are not shown).

As shown in Table 7, the design flood calculated by the DLL method is the largest, followed by the ENE method. The design values of these two methods even exceed that of the conventional method, which is obviously unreasonable. The design flood of the ENE method increased rapidly with the increase of the return period, and almost did not change when the return period was more than 20 years. This change characteristic of the design value of the ENE method is also unrealistic. The design value of ER and ADLL methods are relatively small and close to that of the de-nonstationarity method. However, the uncertainty of these two methods increases rapidly as the recurrence period increases, and even exceeds the magnitude of the design value. On the other hand, these methods all take the effective runoff generation area as the covariate to consider its influence. In MRB, although the effective runoff generation area of a flood event is significantly reduced, the runoff generation process outside of the control area of check dams should not change, that is, the flood peak modulus should not change. The flood peak modulus of the design flood obtained by these reliability-based methods and the conventional method is obviously increased. The de-nonstationarity method proposed in this paper has certain advantages in terms of the magnitude of the calculated design value, the uncertainty, and the consideration of the physical mechanism.

Meantime, the de-nonstationarity method also avoids some problems that are difficult to solve in the current methods. First, the time-varying probability distribution is the basis of the existing methods, but it is extremely difficult to accurately fit the time-varying characteristics of the flood peak discharge series. However, no matter what kind of nonstationarity (trend, abrupt change, etc.) caused by the influence factors, the de-nonstationarity method can remove the influence according to the physical mechanism. Second, in the absence of ergodicity, we cannot use the nonstationary sample series to estimate the population. Third, under nonstationary conditions, the different initial and end time of the calculation period result in different design values. The de-nonstationarity method transforms the nonstationary series into the stationary reconstructed series, thus avoiding this problem.

Check dams, as a kind of engineering measure for water and soil conservation, lead to the nonstationary changes in the flood peak discharges series. This nonstationarity creates unreliable design flood values when they are calculated with the conventional hydrological frequency analysis method. Previous works have attempted to explore the time-varying characteristics of nonstationary hydrological series using purely statistical approaches. Due to the influence of climate change and human activity, it is difficult to accurately describe the characteristics of hydrological time series or predict how the series will evolve in the future using statistical approaches alone. If we cannot guarantee the reliability of the time-varying probability distribution, then the design value under nonstationary conditions is also unreliable.

The conventional method, which requires hydrological series to be stationary, is a mature and robust scientific strategy. In the presence of check dams, nonstationarity in the flood peak discharge series is caused by changes in the effective runoff generation area. By quantifying the power function relationship between the flood peak discharge and the basin area, we can remove the nonstationarity from the flood peak discharge series. Once the reconstructed series is stationary, we can calculate the design quantile via the conventional method, and then obtain the design value of the original nonstationary series according to the effective runoff generation area under design conditions.

We applied this de-nonstationarity method to the nonstationary MRB peak flood discharge series. After removing the nonstationary influence of the effective runoff generation area, the nonstationary flood peak discharge series was transformed into the stationary reconstructed series, and the stationarity is simultaneously achieved in the first and second moments. Because the nonstationarity in the effective runoff generation area has a stronger impact on the second moment of the flood peak discharge series than it does on the first moment, we cannot ignore nonstationarity in higher order moments in a changing hydrological environment. The design flood calculated using the de-nonstationarity method is obviously smaller than that calculated according to the original nonstationary flood peak discharge series. By removing the nonstationarity from the flood peak discharge series, this method can greatly reduce the uncertainty of the design flood. Compared to the four reliability-based methods, the de-nonstationarity method has certain advantages in terms of the magnitude of the calculated design value, the uncertainty, and the consideration of the physical mechanism.

The nonstationary problems caused by climate change or human activities cannot be solved through statistical approaches alone. It is necessary to explore the physical mechanism of the nonstationary change and establish a new method of combining physical mechanism with statistical theory. For other cases, it should also be possible to eliminate the nonstationary effects of influence factors on the hydrological series according to the influence law. Therefore, the exploration of the physical mechanism of the interaction between hydrological elements should be the focus of future research, which can also promote our understanding of hydrological processes.

This research is supported by the National Natural Science Foundation of China (Grant NO.51679184). We are grateful to the editor and the anonymous reviewers for their detailed and insightful comments, which helped us to improve the quality of this paper. We thank LetPub (www.letpub.com) for its linguistic assistance during the preparation of this manuscript.

All relevant data are included in the paper or its Supplementary Information.