Abstract

Quantifying forecast uncertainty is of great importance for reservoir operation and flood control. However, deterministic hydrological forecasts do not consider forecast uncertainty. This study develops a conditional probability model based on copulas to quantify forecast uncertainty. Three updating models, namely auto-regressive (AR) model, AR exogenous input model, and adaptive neuro fuzzy inference system model, are applied to update raw deterministic inflow forecasts of the Three Gorges Reservoir on the Yangtze River, China with lead times of 1d, 2d, and 3d. Results show that the conditional probability model provides a reasonable and reliable forecast interval. The updating models both enhance the forecast accuracy and improve the reliability of probabilistic forecasts. The conditional probability model based on copula functions is a useful tool to describe and quantify forecast uncertainty, and using an updating model is an effective measure to improve the accuracy and reliability of probabilistic forecast.

INTRODUCTION

Flooding is one of the most frequent and severe natural disasters in China; thus, flood forecasting plays a critical role in flood control, disaster reduction, and water resources management (Shen et al. 2015; Wu et al. 2017). Due to the extreme complexity of hydrological processes and the limitation of human knowledge, there inevitably remains uncertainty in the hydrological model output, which mainly results from the errors of model input, model structure, and model parameter (Liu & Gupta 2007). However, traditional and deterministic hydrological forecasting does not consider forecast uncertainty. Decision-makers cannot make a quantitative assessment of the possible risk of flood control and reservoir operation in the absence of uncertainty information. Incorrect decisions may result in casualties and economic losses (Palmer 2002). Therefore, it is essential to quantitatively assess the inherent uncertainty of hydrological forecasts.

To overcome the shortage of deterministic forecasts, many researchers have concentrated on the study of probabilistic forecast in recent years (Krzysztofowicz 1999; Duan et al. 2007; Montanari & Grossi 2008; Todini 2008; Pianosi & Raso 2012; Madadgar & Moradkhani 2014; Khajehei & Moradkhani 2017; He et al. 2018; Liu et al. 2018). The current hydrological uncertainty study mainly focuses on two approaches. The first one is the total factor coupling approach, which quantifies the uncertainties of each link in the rainfall–runoff process, such as rainfall input uncertainty, model structure uncertainty, and model parameter uncertainty. One example of this approach is ensemble forecast that generates a set of forecast values through different inputs, different models, and different parameter optimization methods to express forecast uncertainty (Palmer et al. 2005). Many statistical post-processing methods such as Bayesian model averaging, model conditional processor, and ensemble model output statistics are proposed to accommodate the under-dispersion of ensemble forecasts (Gneiting et al. 2005; Raftery et al. 2005; Zappa et al. 2008; Coccia & Todini 2011). The other one is the total error analysis approach, which directly processes the deterministic forecasts to realize the probabilistic forecasts, such as hydrological uncertainty processor, quantile regression, and Bayesian joint probability model (Krzysztofowicz & Kelly 2000; Weerts et al. 2010; Zhao et al. 2015). The core of the probabilistic forecast literature mentioned above is to establish multivariate joint distribution, where the multivariate Normal distribution is used commonly (Biondi & Todini 2018). Considering that hydrological variables often do not follow the Normal distribution, many data transformation methods such as Box-Cox transformation (Box & Cox 1964) or Normal quantile transform (Kelly & Krzysztofowicz 1997) have been proposed to handle any marginal distributions.

Most previous probabilistic flood forecasting models only consider the uncertainty of the hydrologic model structure and parameters; this study will consider the uncertainty of hydrologic model structure and parameters, as well as model input. So, a conditional probability model based on copula is developed to post-process daily streamflow forecasts and estimate uncertainty in this study. For the purpose of improving the accuracy and reliability of daily streamflow forecasts, three updating models, including auto-regressive (AR), AR exogenous input (ARXM), and adaptive neuro fuzzy inference system (ANFIS), are applied to correct forecasting errors. The effect of three updating models on the accuracy and reliability of probabilistic forecast is also explored and discussed. The remainder of this paper is organized as follows: the conditional probability model, three updating models, and evaluation criteria are explained in the next section. The study area and the hydrological forecasting model are presented in the section ‘Case studies’. In the following section, the results and discussion are given. Finally, the main conclusions are summarized in the last section.

METHODOLOGY

Copula function

The Copula function, introduced by Sklar (1959), is capable of constructing the multivariate joint distributions through marginal distribution and correlation structure, which is simple and flexible. It can construct multi-dimensional joint distributions through marginal distributions of several variables and their correlations. Suppose that the marginal distributions of n-dimensional random variables are , where n is the number of random variables and is the value of random variables . Let be the joint distribution function of random variables . According to Sklar's theorem (Sklar 1959), the multivariate cumulative distribution function (CDF) can be written as follows: 
formula
(1)
where C is the copula function and is the parameter of C. If are continuous, then C is unique.

There are different families of copulas that have been employed in a variety of applications (Nelsen 2006). Of all the copula families, Archimedean copulas find a wide range of applications in hydrology, since they can be easily constructed and possess many good properties like simplicity and adaptability. Gumbel-Hougaard, Clayton, and Frank copulas are three commonly used one-parameter Archimedean copulas, which are widely applied in hydrology by many authors (Zhang & Singh 2006; Duan et al. 2016; Zhong et al. 2017; Yin et al. 2018a, 2018b), and chosen in this study to establish the joint distribution between forecasts and observations. The expression of these copulas and the correlation between their parameters and Kendall correlation coefficients are listed in Table 1. The parameters of the two-dimensional Archimedean copulas can be estimated by Kendall correlation coefficients.

Table 1

Three commonly used two-dimensional Archimedean copula functions

Copula  θ ∈  
Gumbel-Hougaard    
Clayton    
Frank    
Copula  θ ∈  
Gumbel-Hougaard    
Clayton    
Frank    
To obtain the optimal copula function, the root mean square error (RMSE) is used to measure the goodness of fit between the empirical and the theoretical joint distribution (Zhang & Singh 2007), and defined as follows: 
formula
(2)
where and are the empirical and theoretical joint CDF values, respectively. The smaller the value of RMSE is, the better the copula function.

Conditional probability

Let predictand H be the observed discharge whose realization h is being forecast. Let estimator S be the output discharge generated by a corresponding deterministic forecast model whose realization s constitutes a point estimate of H. Let random variable H0 represent the observed discharge at the time k = 0 when the forecast is prepared; then Hk (k = 1, 2, …, N) is the observed discharge and Sk (k = 1, 2, …, N) is the corresponding deterministic forecast discharge at lead time k. What the rational decision-maker then needs is not a single number sk, but the distribution function of predictand Hk conditional on Sk=sk (Liu et al. 2018).

Let and be the marginal CDFs, respectively. The corresponding probability density functions (PDFs) are and . According to Sklar's theorem, the joint CDF of the H and S can be expressed as follows: 
formula
(3)
For a given S=s, the conditional CDF of observation h can be derived as follows: 
formula
(4)
The corresponding conditional PDF of h can be expressed as follows: 
formula
(5)
where is the PDF of the twodimensional copula function.

The solution of conditional probability model based on copula (Equation (5)) is given by the following steps:

  • 1.

    In the model calibration period, the parameters of marginal distribution and copula function are estimated by the data series of observations (H) and forecast values (S), respectively.

  • 2.
    In the real-time flood forecasting, the observations are unknown for the lead times of 1d, 2d, and 3d. We discretize and determine the distribution range [hmin, hmax] for the observation H, according to its real flood probability distribution characteristics. Then the discrete h can be expressed as follows: 
    formula
    (6)
     
    formula
    (7)
    where hmin and hmax are the minimum and maximum observation, n is the number of discretization, and is the discretization spacing.
  • 3.

    Calculate the density function of copula and the density function value , according to the forecast values (sk) and the discrete value h.

  • 4.

    Calculate the conditional density function value by Equation (5).

The median value of the random variable H is used as a point estimate result. The median value is the 50% quantile of the random variable H, having the following form: 
formula
(8)

For a given confidence level , let and be the 5% and 95% quantiles of the random variable h, respectively. Then the 90% confidence interval of h is used to reflect the forecast uncertainty.

Updating models

AR updating model

The most widely used and the simplest updating model is the AR model, whose procedure consists of three steps. First, the forecast errors series and the corresponding distance average series can be calculated by 
formula
(9)
 
formula
(10)
where is the forecast error at the current time t, ht and st denote the observed and forecasted flow discharge at the current time t, respectively, is the mean value of forecast error series, and is the distance average series of et, where the mean value of equals zero.
Then, the estimates of for one-step ahead can be expressed as follows: 
formula
(11)
where are the coefficients of the AR model, which can be estimated by the ordinary least square method; is the estimate of , and p is the order of the AR model, which can be determined by autocorrelation and partial autocorrelation and Bayesian information criterion (BIC).
Finally, the updated forecasts can be obtained by adding up the forecast error estimates and the original forecast value, having the following form: 
formula
(12)
where is the updated forecast and is the output of the deterministic hydrological model.
Equation (11) only describes the one-step-ahead updating procedure, while the lead-2 and lead-3 updating procedures can be expressed as follows: 
formula
(13)
 
formula
(14)
where and are calculated by Equations (11) and (13), respectively.

ARXM model

In the ARXM updating procedure, the recent values of the observed discharge together with the non-updated outputs of the hydrological model are used as inputs of the ARXM model to produce the updated discharge forecast. The ARXM updating procedure for one-step lead-time can be expressed as follows: 
formula
(15)
where denotes the one-step-ahead updated flow discharge, denotes the antecedent observed flow discharge, denotes the current and recent output of the hydrological model, and p and q are the orders of auto-regressive and the exogenous input parts of the ARXM, respectively. are the coefficients of the ARXM model, which are estimated by the ordinary least square method.
For a lead-time , the updated flow discharge forecasts can be given by 
formula
(16)

The ARXM updating model differs from the AR updating model. In the ARXM updating procedure, the model acquires the updated forecasts directly, whereas during the AR updating process, the forecast errors at the current time are firstly estimated, and then added to the original model output to obtain the updated forecasts.

Adaptive neural fuzzy inference system

According to the types of fuzzy reasoning and fuzzy if-then rules, the ANFIS can be categorized into three types (Jang 1993). In this study, Takagi and Sugeno's fuzzy if-then rules are used, according to which the output of the whole fuzzy inference system is a weighted average of each rule's output, the weight values of each rule output are determined by the membership function (a mathematical tool used to represent the ‘degree of truth’ of the element belonging to a fuzzy set).

As shown in Figure 1, the simplest example is given to explain the structure of ANFIS. We assume that the ANFIS has two inputs x1 and x2 and the fuzzy system has only two if-then rules of Takagi and Sugeno's type. Then, its typical fuzzy inference rule can be described as follows: 
formula
(17)
 
formula
(18)
where A1 and A2 are the linguistic label (small or large) of the inputs and B1 and B2 are the linguistic label (small or large) of the inputs .
Figure 1

The structure of a two-rule two input ANFIS (adapted from Jang 1993).

Figure 1

The structure of a two-rule two input ANFIS (adapted from Jang 1993).

As shown in Figure 1, the ANFIS network structure consists of five layers, the functions of each layer node are explained as follows:

Layer 1: Each node in this layer is an adaptive node, which is denoted by a square node. The node function describes the membership degree to which the given input x1 belongs to the corresponding fuzzy set Ai, having the following form: 
formula
(19)
where is the membership function of the fuzzy set Ai, the Gaussian function is chosen in this study.
Layer 2: Each node in this layer is a fixed node, denoted by a circle node. The output of the node is the multiple of all the incoming signals. 
formula
(20)
Layer 3: The output of the ith node is the ratio of the ith rule confidence wi to the sum of all the fuzzy if-then rules confidence, having the following form: 
formula
(21)
Layer 4: Each node in this layer is also an adaptive node. The output of each node in this layer can be expressed as follows: 
formula
(22)
Layer 5: The node in this layer is a fixed node. It calculates the sum of all incoming signals as the total final output: 
formula
(23)
The input and output of the ANFIS updating model in each lead time are the same as the ARXM updating model. For the lead times of 1d, 2d, and 3d, the relationship between the input and output of the ANFIS model can be expressed as follows, respectively: 
formula
(24)
 
formula
(25)
 
formula
(26)
where FANFIS denotes a non-linear function.

Evaluation criteria

Deterministic forecast

  • 1.
    Nash-Sutcliffe efficiency (NSE) 
    formula
    (27)
    where and denote the observed and forecasted flow discharges at time t, respectively; and is the average value of observed flows. The closer the value of NSE is to 1, the better the model performs. NSE indicators are sensitive to the evaluation of high flow, while being insensitive to the additive and proportional differences between forecasts and observations.
  • 2.
    Relative error (RE) 
    formula
    (28)
    where RE represents the water balance error, which is an evaluation of the total water volume. The value closer to zero indicates a better forecast. The forecasting errors of flood volume are also very important for large reservoirs.
  • 3.
    RMSE 
    formula
    (29)
    where RMSE represents the standard deviation of the differences between the observed values and simulated values, which is sensitive to outliers. The value closer to zero indicates a better forecast. In the evaluation of medium and high flows, the RMSE indicator can better reflect the pros and cons of the model.
  • 4.
    Mean absolute error (MAE) 
    formula
    (30)
    where MAE is the average absolute error between the observed and forecasted value, it can reflect the forecast error of each point more intuitively. Small MAE values are preferred. In the evaluation of low and medium flows, the MAE indicator can better reflect the pros and cons of the model. The units of RMSE and MAE are m3/s.

Probabilistic forecast

In this section, continuous ranked probability score (CRPS), quantile-quantile (Q-Q) plot, and other indicators are used to evaluate the quality of probabilistic forecast.

  • 1.

    CRPS

CRPS is the most widely used criterion to give a comprehensive evaluation of the quality of the probabilistic forecast, which can address the reliability and sharpness of predictive distribution simultaneously (Gneiting & Raftery 2007). The CRPS is defined as follows: 
formula
(31)
where is the cumulative probability function of the forecast at time t; is the observed flow at time t; and denotes the Heaviside step function, which is calculated by 
formula
(32)

For the deterministic forecasts, the CRPS reduces to generalize the absolute error, which is equivalent to the MAE (Zhao et al. 2015; Liu et al. 2018). Thus, it can also be used to directly compare the accuracy of the deterministic forecast and probabilistic forecast. Smaller CRPS values indicate better performance.

  • 2.

    Q-Q plot

The reliability of the probabilistic forecast can be evaluated by the Q-Q plot (Laio & Tamea 2007). The closer the Q-Q plot to the 1:1 line indicates a more reliable probabilistic forecast. In order to evaluate the distance between the Q-Q plot and 1:1 line, the indicator presented by Renard et al. (2010) is used and calculated as follows: 
formula
(33)
where is the value that the predictive CDF derives at the observation and is the non-exceedance uniform probability of the observation. The value of the ranges from 0 to 1, large values mean it has greater reliability.
  • 3.

    Assess the 90% confidence interval

The 90% confidence interval is used to reflect the uncertainties of forecast results. To assess the quality of the 90% confidence interval, four indices were evaluated: containing ratio (CR), average band-width (B), average relative band-width (RB), and percentage of observations bracketed by the unit confidence interval (PUCI). These indices are used and calculated as follows: 
formula
(34)
where nc is the number of measured data contained in the prediction interval. It describes the ratio of the number of observed points within the confidence interval to the total number of the observed points. The larger the CR value indicates the higher percentage of the observed flows within the confidence interval 
formula
(35)
where denotes the upper bound of the 90% confidence interval corresponding to the 95% quantile and denotes the lower bound of the 90% confidence interval corresponding to the 5% quantile 
formula
(36)
where is the observed flow at time t. The average band-width and the average relative band-width are also commonly used to assess the quality of the confidence interval. For a given confidence level , it is best that the values of B and RB should be as small as possible.
According to the above, we can know that the smaller RB value indicates that the average relative band-width of the confidence interval is narrower; the larger the CR value means higher percentage of the observations bracketed by the confidence interval. However, Xiong et al. (2009) indicated that as the CR increased, the RB value of the confidence interval also increased. The two indices are often contradictory when assessing the quality of the 90% confidence interval. Therefore, a comprehensive index of the percentage of observations bracketed by the unit confidence interval (PUCI) proposed by Li et al. (2011) was used to evaluate the quality of the uncertainty interval estimate. It can be expressed as follows: 
formula
(37)
The larger the PUCI value, the better the performance of the 90% confidence interval is.

CASE STUDIES

Study area and data

The Yangtze River is the largest river in China and the third longest river in the world. The Three Gorges Reservoir (TGR), located in the upper reaches of the Yangtze River, was selected as a case study. It is the largest and the most important water conservancy project in China. The river channel above the dam site of the TGR is 4,500 km long with a contributing drainage area of . The TGR plays an important role in flood control in the middle and lower sections of the Yangtze River, with the combined benefits of power generation, navigation improvement, etc.

The TGR intervening basin, shown in Figure 2, has a catchment area of 55,907 km2, accounting for 5.6% of the TGR basin area. The mean annual rainfall and evaporation in the basin are 1,100 mm and 900 mm, respectively. The mean annual inflow is 14,300 m3/s. The observed and forecasted daily inflows of Cuntan and Wulong hydrological stations, and the daily rainfall data in the intervening basin from 2012 to 2015 are collected for the inflow forecasting of the TGR. The data series from 2012 to 2014 are used for calibration, and the data of 2015 are used for validation. Numerical precipitation forecasting is necessary to prolong the forecast horizon. The numeric weather prediction daily precipitation predictions from the European Center for Medium-range Weather Forecast (ECMWF) in lead times of 1d, 2d, and 3d are used as a hydrological model input in this study. The statistical characteristics of calibration and validation data are also shown in Figure 3. The maximum and minimum flow during the validation period are included in the calibration period, and the samples in the validation period have similar statistical characteristics during the calibration period, indicating that the calibration and validation period are reasonably divided and can be used for model calibration.

Figure 2

Sketch map of the TGR intervening basin.

Figure 2

Sketch map of the TGR intervening basin.

Figure 3

Daily inflow series of the TGR and its statistical box plot.

Figure 3

Daily inflow series of the TGR and its statistical box plot.

Hydrological forecasting model

The inflow of TGR consists of three components, the main streamflow measured by the Cuntan gauge station, the tributary inflow from the Wu River measured by the Wulong gauge station, and the rainfall runoff of the TGR intervening basin, as shown in Figure 2. The flowchart of deterministic and probabilistic flow forecasting for the TGR is shown in Figure 4, in which the deterministic hydrological forecasting of TGR is implemented by a mixed model. The rainfall–runoff relationships of the intervening basin are modeled by the Xinanjiang model (Zhao 1980). The Muskingum method is applied for flood routing (Cunge 1969). Then, the inflow forecast of the TGR can be calculated by 
formula
(38)
where Qsx(t), Qct(t), and Qwl(t) denote the inflow of TGR, the flow discharges of Cuntan station and Wulong station at time t, respectively.
Figure 4

Flowchart of deterministic and probabilistic flow forecasting for the TGR.

Figure 4

Flowchart of deterministic and probabilistic flow forecasting for the TGR.

The Xinanjiang model proposed by Zhao (1980) is suitable for humid and semi-humid regions and has been widely used for flood simulation and forecasting in China. The model is mainly characterized by the concept of runoff formation as natural storage in humid regions, which implies that runoff is produced after the soil water content reaches field capacity. The storage capacity curve is used to solve the problem of the uneven space distribution of the runoff formation throughout the basin. There are 15 parameters, and each parameter has a defined physical meaning. A more detailed description of the model is available in Zhao (1980).

In the model simulation calibration and verification periods, the observed flood series of the Cuntan and Wulong hydrological stations as well as the observed rainfall of the interval basin area are used as model inputs. The output of the model is the inflow of the TGR. The parameters of the Xinanjiang model and Muskingum method are optimized by the genetic algorithm (Goldberg 1989). The simulation results for the NSE and RE are 99.05% and −0.57% in the calibration period, 98.30% and −1.57% in the validation period, respectively. The simulation results indicate that the established model can be used as a forecasting scheme in practice.

RESULTS AND DISCUSSION

The order of the three updating models

The order of the three updating models can be determined by using the autocorrelation function (ACF), partial ACF (PACF) and the BIC. In the ARXM and ANFIS updating procedure, the recent values of the observed discharge together with the non-updated outputs of the hydrological model are used as inputs, whereas during the AR updating process, the distance average series of the forecast error series is used as its input, the order of models is determined separately. The BIC results of the AR model under different model orders are plotted in Figure 5 (red line). The ACF and PACF of the error time series are shown in Figure 6. In Figure 5, the BIC has the minimum value at the order p = 5. In Figure 6, the ACF of the error time series tends to be zero with the increase in the lags. However, the partial autocorrelation coefficient is truncated, which occurs roughly at the order p = 5. Therefore, the order of the AR model is determined as 5.

Figure 5

BIC results. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2019.094.

Figure 5

BIC results. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2019.094.

Figure 6

The ACF and PACF of errors time series.

Figure 6

The ACF and PACF of errors time series.

For the ARXM and ANFIS updating models, the order p and q need to be determined. The ACF and PACF of the observed TGR inflow series are calculated and shown in Figure 7. In Figure 7, the PACF is truncated at the order p = 3. The BIC results of different orders q at order p = 3 are shown in Figure 5 (blue line), on which the BIC has the minimum value at the order q = 3. Thus, the order p = 3 and q = 3 are used to update the forecasts by the ARXM and ANFIS models.

Figure 7

The ACF and PACF of the observed TGR inflow series.

Figure 7

The ACF and PACF of the observed TGR inflow series.

The marginal distributions

First, the marginal distributions of the observed flow, raw forecasts, and updated forecasts need to be determined. Seven commonly used distributions, including Gamma, Lognormal, Generalized Extreme Value (GEV), Pearson type III (P-III), Gumbel, Weibull, and Log-logistic, are chosen as candidate distributions to fit the marginal distributions of all variables. The parameters of each distribution are estimated by the maximum likelihood estimation method. The modified Shapiro-Wilk test developed by Ashkar & Aucoin (2012) is applied to evaluate the goodness of fit of the probability distributions. The test results of each candidate distribution are listed in Table 2. The best marginal distribution is the one that has the maximum test values. As shown in Table 2, the P-III distribution performs the best fitting with the maximum test values for all variables. The upper tail of the observed flow series is relatively large, which is consistent with the characteristics of the P-III distribution and has been widely used by many authors (Zhong et al. 2017; Yin et al. 2018a, 2018b). The estimated parameters of P-III distributions for all variables are listed in Table 3.

Table 2

Test results for seven selected marginal distributions

Variable Gamma Lognormal GEV P-III Gumbel Weibull Log-logistic 
H 0.9150 0.9278 0.9272 0.9412 0.9007 0.9049 0.9164 
S1 0.9071 0.9246 0.9340 0.9451 0.8918 0.8961 0.9142 
S2 0.9029 0.9183 0.9282 0.9355 0.8882 0.8928 0.9079 
S3 0.9060 0.9183 0.9260 0.9320 0.8927 0.8968 0.9070 
SAR1 0.9125 0.9278 0.9304 0.9442 0.8979 0.9016 0.9169 
SAR2 0.9082 0.9217 0.9247 0.9364 0.8959 0.8972 0.9111 
SAR3 0.9106 0.9210 0.9222 0.9336 0.9016 0.8998 0.9098 
SARXM1 0.9129 0.9274 0.9295 0.9428 0.8982 0.9023 0.9164 
SARXM2 0.9101 0.9223 0.9229 0.9348 0.8962 0.9005 0.9111 
SARXM3 0.9162 0.9249 0.9224 0.9345 0.9039 0.9076 0.9127 
SANFIS1 0.9151 0.9273 0.9259 0.9395 0.9011 0.9053 0.9158 
SANFIS2 0.9166 0.9278 0.9257 0.9383 0.9028 0.9071 0.9159 
SANFIS3 0.9178 0.9265 0.9237 0.9356 0.9042 0.9100 0.9142 
Variable Gamma Lognormal GEV P-III Gumbel Weibull Log-logistic 
H 0.9150 0.9278 0.9272 0.9412 0.9007 0.9049 0.9164 
S1 0.9071 0.9246 0.9340 0.9451 0.8918 0.8961 0.9142 
S2 0.9029 0.9183 0.9282 0.9355 0.8882 0.8928 0.9079 
S3 0.9060 0.9183 0.9260 0.9320 0.8927 0.8968 0.9070 
SAR1 0.9125 0.9278 0.9304 0.9442 0.8979 0.9016 0.9169 
SAR2 0.9082 0.9217 0.9247 0.9364 0.8959 0.8972 0.9111 
SAR3 0.9106 0.9210 0.9222 0.9336 0.9016 0.8998 0.9098 
SARXM1 0.9129 0.9274 0.9295 0.9428 0.8982 0.9023 0.9164 
SARXM2 0.9101 0.9223 0.9229 0.9348 0.8962 0.9005 0.9111 
SARXM3 0.9162 0.9249 0.9224 0.9345 0.9039 0.9076 0.9127 
SANFIS1 0.9151 0.9273 0.9259 0.9395 0.9011 0.9053 0.9158 
SANFIS2 0.9166 0.9278 0.9257 0.9383 0.9028 0.9071 0.9159 
SANFIS3 0.9178 0.9265 0.9237 0.9356 0.9042 0.9100 0.9142 
Table 3

Estimated parameters of P-III distributions for all variables

Variables    Variables    
H 1.09 0.000096 2,919.72 SARXM1 1.03 0.000093 3,059.47 
S1 0.94 0.000087 3,340.38 SARXM2 1.07 0.000096 2,885.52 
S2 0.97 0.000091 3,150.95 SARXM3 1.25 0.000108 2,679.98 
S3 1.09 0.000104 2,915.85 SANFIS1 1.11 0.000098 2,857.52 
SAR1 1.01 0.000091 3,151.49 SANFIS2 1.17 0.000102 2,702.51 
SAR2 1.05 0.000096 3,247.06 SANFIS3 1.25 0.000107 2,495.05 
SAR3 1.21 0.000111 3,351.23     
Variables    Variables    
H 1.09 0.000096 2,919.72 SARXM1 1.03 0.000093 3,059.47 
S1 0.94 0.000087 3,340.38 SARXM2 1.07 0.000096 2,885.52 
S2 0.97 0.000091 3,150.95 SARXM3 1.25 0.000108 2,679.98 
S3 1.09 0.000104 2,915.85 SANFIS1 1.11 0.000098 2,857.52 
SAR1 1.01 0.000091 3,151.49 SANFIS2 1.17 0.000102 2,702.51 
SAR2 1.05 0.000096 3,247.06 SANFIS3 1.25 0.000107 2,495.05 
SAR3 1.21 0.000111 3,351.23     

Construction of the joint distribution based on copula

Three commonly used Archimedean copula functions, including Gumbel, Clayton, and Frank, are selected to construct the joint distribution of the forecasts and observations. The Kendall correlation coefficient between observed and forecasted series is calculated, and then the parameters of three copula functions are computed, according to the relationship between and . The RMSE values for evaluating the goodness of fit between theoretical and empirical joint distributions are calculated for each copula function. The results listed in Table 4 show that the Gumbel copula has the best fit with the smallest RMSE values for lead times of 1d, 2d, and 3d. Thus, the Gumbel copula is applied for constructing the joint distribution between observations and forecasts. It can also be seen in Table 4 that the Kendall correlation coefficient between observations and updated forecasts is greater than the original ones, which indicates that the Kendall correlation coefficient between observed series and forecasted series has been improved after error correction. Furthermore, it can be seen in Table 4 that the Kendall correlation coefficient decreases as the lead time increases as expectedly, since forecast quality decreases as the lead time increases.

Table 4

Estimated parameters of copula function and corresponding RMSE values

Variables Copula   RMSE 
H,S1 Gumbel-Hougaard 0.940 16.67 0.0034 
Clayton 31.34 0.0120 
Frank 64.99 0.0119 
H,S2 Gumbel-Hougaard 0.904 10.37 0.0029 
Clayton 18.75 0.0065 
Frank 39.78 0.0076 
H,S3 Gumbel-Hougaard 0.877 8.10 0.0055 
Clayton 14.20 0.0115 
Frank 30.66 0.0126 
H,SAR1 Gumbel-Hougaard 0.954 21.90 0.0030 
Clayton 41.80 0.0098 
Frank 85.92 0.0089 
H,SAR2 Gumbel-Hougaard 0.908 10.91 0.0035 
Clayton 19.83 0.0122 
Frank 41.94 0.0115 
H,SAR3 Gumbel-Hougaard 0.883 8.52 0.0069 
Clayton 15.05 0.0111 
Frank 32.36 0.0135 
H,SARXM1 Gumbel-Hougaard 0.955 22.36 0.0042 
Clayton 42.71 0.0103 
Frank 87.75 0.0091 
H,SARXM2 Gumbel-Hougaard 0.905 10.52 0.0039 
Clayton 19.04 0.0105 
Frank 40.36 0.0075 
H,SARXM3 Gumbel-Hougaard 0.879 8.28 0.0035 
Clayton 14.56 0.0061 
Frank 31.38 0.0061 
H,SANFIS1 Gumbel-Hougaard 0.957 23.33 0.0081 
Clayton 44.65 0.0164 
Frank 91.62 0.0124 
H,SANFIS2 Gumbel-Hougaard 0.919 12.38 0.0026 
Clayton 22.76 0.0079 
Frank 47.82 0.0097 
H,SANFIS3 Gumbel-Hougaard 0.887 8.87 0.0035 
Clayton 15.75 0.0060 
Frank 33.77 0.0085 
Variables Copula   RMSE 
H,S1 Gumbel-Hougaard 0.940 16.67 0.0034 
Clayton 31.34 0.0120 
Frank 64.99 0.0119 
H,S2 Gumbel-Hougaard 0.904 10.37 0.0029 
Clayton 18.75 0.0065 
Frank 39.78 0.0076 
H,S3 Gumbel-Hougaard 0.877 8.10 0.0055 
Clayton 14.20 0.0115 
Frank 30.66 0.0126 
H,SAR1 Gumbel-Hougaard 0.954 21.90 0.0030 
Clayton 41.80 0.0098 
Frank 85.92 0.0089 
H,SAR2 Gumbel-Hougaard 0.908 10.91 0.0035 
Clayton 19.83 0.0122 
Frank 41.94 0.0115 
H,SAR3 Gumbel-Hougaard 0.883 8.52 0.0069 
Clayton 15.05 0.0111 
Frank 32.36 0.0135 
H,SARXM1 Gumbel-Hougaard 0.955 22.36 0.0042 
Clayton 42.71 0.0103 
Frank 87.75 0.0091 
H,SARXM2 Gumbel-Hougaard 0.905 10.52 0.0039 
Clayton 19.04 0.0105 
Frank 40.36 0.0075 
H,SARXM3 Gumbel-Hougaard 0.879 8.28 0.0035 
Clayton 14.56 0.0061 
Frank 31.38 0.0061 
H,SANFIS1 Gumbel-Hougaard 0.957 23.33 0.0081 
Clayton 44.65 0.0164 
Frank 91.62 0.0124 
H,SANFIS2 Gumbel-Hougaard 0.919 12.38 0.0026 
Clayton 22.76 0.0079 
Frank 47.82 0.0097 
H,SANFIS3 Gumbel-Hougaard 0.887 8.87 0.0035 
Clayton 15.75 0.0060 
Frank 33.77 0.0085 

Comparison of deterministic and probabilistic forecasting results

Deterministic forecasting results

According to the distribution range of the observation, let hmin = 3,000 m3/s and hmax = 90,000 m3/s (the cumulative probability of which is above 0.99999). It is theoretically better to select a small discrete spacing. However, smaller discrete spacing introduces more computational complexity into the calculation. Therefore, the discrete spacing is selected in this study. Then the discrete H of each time step can be calculated by Equation (6), where the number of discretization is n = 43,501. Based on the principle of mathematical statistics, the probability median forecast is calculated by formula (8) and used as the point estimate results of probabilistic forecast. Four evaluation indexes, namely NSE, RE, RMSE, and MAE, are used to evaluate the accuracy of the median forecast. The comparison with the raw forecast is presented in Figure 8, which shows that in terms of all evaluation indicators, the accuracy of probability median forecast is significantly improved compared with the raw forecast. Compared with the raw forecast results, the RE of the probability median forecast are reduced by 0.12%, 3.38%, and 5.01%, and the NSE are improved by 0.1%, 0.6%, and 1.6% for lead times of 1d, 2d, and 3d in the validation period, respectively.

Figure 8

Comparison of probability median forecast and raw forecast for lead times of 1d, 2d, and 3d.

Figure 8

Comparison of probability median forecast and raw forecast for lead times of 1d, 2d, and 3d.

The deterministic evaluation indexes of raw forecast and three updating models for lead times of 1d, 2d, and 3d in calibration and verification periods are listed in Table 5. The comparison between the forecast results of the three updating models and the original model is also shown in Figure 9. It is obvious that the performances of three updating models are superior to that of the raw forecasts in terms of almost all evaluation indicators. The ANFIS updating model performs the best for almost all lead times, while the ARXM updating model performs slightly better than the AR model for 1d, 2d, and 3d lead times. The ANIFIS model performs the best because it not only can effectively identify the non-linear relationship in hydrological process but also effectively classify the inputs into high and low levels, which can consider the relationship between errors and flood magnitude. In addition, according to the results of NSE and RE indexes, the improvement degree becomes more obvious with the increase in lead time. However, for the RMSE and MAE indexes, the improvement degree decreases with lead time increases. For the ANFIS updating model, the MAE indicator in 1d, 2d and 3d lead times has improved compared with the original forecast, while for the AR and ARXM models, the improvement in the MAE index is no longer achieved for long lead times (3d) compared to the original forecast. This indicates that the ability of the updating model for correcting medium and low flow is gradually lost as the lead time increases. Comparing the results of indicators in Figures 8 and 9, it can be found that the traditional updating models, including AR, ARXM, and ANFIS, are better than the copula-based probability model.

Table 5

Deterministic evaluation indexes of raw forecast and three updating models for lead times of 1d, 2d, and 3d in calibration and verification periods

Lead time Model NSE (%)
 
RE (%)
 
MAE
 
RMSE
 
Cal. Val. Cal. Val. Cal. Val. Cal. Val. 
1d RAW 98.61 98.35 −0.41 −0.79 742.52 540.05 1,241 868 
AR 98.99 98.92 0.06 0.05 546.24 411.58 1,058 702 
ARXM 99.00 99.03 −0.57 −0.71 522.11 380.46 1,053 665 
ANFIS 99.10 99.06 0.00 0.23 508.28 373.40 1,000 654 
2d RAW 95.59 95.57 −2.72 −4.23 1,209.10 818.53 2,209 1,422 
AR 96.16 96.33 0.07 0.73 1,117.40 765.04 2,062 1,295 
ARXM 96.19 96.60 −0.85 −2.36 1,096.22 756.68 2,055 1,246 
ANFIS 96.43 96.79 0.03 −0.96 1,044.33 715.91 1,988 1,211 
3d RAW 90.43 91.94 −5.81 −5.97 1,720.63 1,042.29 3,254 1,918 
AR 91.55 93.30 0.11 1.27 1,741.54 1,103.87 3,059 1,749 
ARXM 92.09 93.75 0.29 0.99 1,644.52 1,059.29 2,959 1,689 
ANFIS 92.60 93.60 0.02 0.31 1,542.81 997.71 2,863 1,709 
Lead time Model NSE (%)
 
RE (%)
 
MAE
 
RMSE
 
Cal. Val. Cal. Val. Cal. Val. Cal. Val. 
1d RAW 98.61 98.35 −0.41 −0.79 742.52 540.05 1,241 868 
AR 98.99 98.92 0.06 0.05 546.24 411.58 1,058 702 
ARXM 99.00 99.03 −0.57 −0.71 522.11 380.46 1,053 665 
ANFIS 99.10 99.06 0.00 0.23 508.28 373.40 1,000 654 
2d RAW 95.59 95.57 −2.72 −4.23 1,209.10 818.53 2,209 1,422 
AR 96.16 96.33 0.07 0.73 1,117.40 765.04 2,062 1,295 
ARXM 96.19 96.60 −0.85 −2.36 1,096.22 756.68 2,055 1,246 
ANFIS 96.43 96.79 0.03 −0.96 1,044.33 715.91 1,988 1,211 
3d RAW 90.43 91.94 −5.81 −5.97 1,720.63 1,042.29 3,254 1,918 
AR 91.55 93.30 0.11 1.27 1,741.54 1,103.87 3,059 1,749 
ARXM 92.09 93.75 0.29 0.99 1,644.52 1,059.29 2,959 1,689 
ANFIS 92.60 93.60 0.02 0.31 1,542.81 997.71 2,863 1,709 
Figure 9

Deterministic evaluation indexes of three updating models for lead times of 1d, 2d, and 3d. Note: AR1, ARXM1, ANFIS1 models with 1d lead time, AR2, ARXM2, ANFIS2 models with 2d lead time, AR3, ARXM3, ANFIS3 models with 3d lead time.

Figure 9

Deterministic evaluation indexes of three updating models for lead times of 1d, 2d, and 3d. Note: AR1, ARXM1, ANFIS1 models with 1d lead time, AR2, ARXM2, ANFIS2 models with 2d lead time, AR3, ARXM3, ANFIS3 models with 3d lead time.

Probabilistic forecasting results

The conditional probability density and distribution curve of the observation H, given Sk = sk, can be obtained based on the conditional probability for given forecasts. Take 1 July 2015 (S1 = 38,462 m3/s, S2 = 33,469 m3/s, and S3 = 26,386 m3/s) as an example, the PDF and CDF of observation H are presented in Figure 10. As can be seen in Figure 10, the conditional density function becomes fatter with the increase in lead time, indicating that the flood forecasting uncertainty increases as the lead time increases.

Figure 10

The PDF and CDF of observation H on 1 July 2015.

Figure 10

The PDF and CDF of observation H on 1 July 2015.

In order to evaluate the reliability of the probabilistic forecast, the predictive Q-Q plots for three updating models with different lead times are shown in Figure 11. As the Kolmogorov–Smirnov test is highly sensitive when applied to the large data sets (Serinaldi 2009), none of the Q-Q plots lie within the significance bands. However, overall the Q-Q plots for the updating models are closer to the 1:1 line than the original model for all lead times. It is noted that we cannot obtain enough information from the visual inspection of the Q-Q plot to evaluate their reliability. Therefore, the CPRS and are used to quantitatively measure the reliability of the probabilistic forecasts. The values of CRPS and for the probability raw forecast and updating models for lead times of 1d, 2d, and 3d are listed in Table 6. The values of for three updating models are closer to 1 than those of the original model, which implies that the probabilistic forecasts are more reliable after forecast correction. It is also shown that the CRPS values increase as the lead time increases, which indicates that the probabilistic forecast is less sharp and reliable at longer lead times. Compared with the corresponding raw deterministic forecast, the CRPS values of the probabilistic forecast in the validation period are reduced by 25.76%, 26.86%, and 26.37% for lead times of 1d, 2d, and 3d, respectively. It indicates that the conditional probability model based on copula is a useful tool for quantifying and assessing forecast uncertainty. Table 6 also shows that as the accuracy of deterministic forecast increases after error correction, the corresponding probabilistic forecast becomes more reliable. This is because the conditional probability model achieves probabilistic forecast on the basis of given deterministic forecast and the reliability of probabilistic forecast depends largely on the accuracy of deterministic forecast. In terms of the comprehensive index CRPS, the ANFIS model performs the best in almost all lead times.

Table 6

The values of CRPS and of probability raw forecast and updating models for lead times of 1d, 2d, and 3d in calibration and verification periods

Lead time Model Probability forecast/CRPS
 
Deterministic forecast/MAE
 
CRPS decline (%)
 

 
Cal. Val. Cal. Val. Cal. Val. Cal. Val. 
1d Raw 513.48 400.91 742.52 540.05 30.85 25.76 0.916 0.826 
AR 406.12 312.18 546.24 411.58 25.65 24.15 0.928 0.875 
ARXM 402.07 294.15 522.11 380.46 22.99 22.69 0.917 0.868 
ANFIS 391.48 287.24 508.28 373.40 22.98 23.07 0.928 0.887 
2d Raw 919.38 598.66 1,209.10 818.53 23.96 26.86 0.900 0.861 
AR 853.29 584.17 1,117.40 765.04 23.64 23.64 0.913 0.892 
ARXM 847.58 572.37 1,096.22 756.68 22.68 24.36 0.912 0.886 
ANFIS 814.47 549.51 1,044.33 715.91 22.01 23.24 0.930 0.891 
3d Raw 1,365.86 767.42 1,720.63 1,042.29 20.62 26.37 0.911 0.855 
AR 1,293.59 749.53 1,741.54 1,103.87 25.72 32.10 0.930 0.884 
ARXM 1,265.91 766.36 1,644.52 1,059.29 23.02 27.65 0.932 0.891 
ANFIS 1,212.78 763.28 1,542.81 997.71 21.39 23.50 0.930 0.888 
Lead time Model Probability forecast/CRPS
 
Deterministic forecast/MAE
 
CRPS decline (%)
 

 
Cal. Val. Cal. Val. Cal. Val. Cal. Val. 
1d Raw 513.48 400.91 742.52 540.05 30.85 25.76 0.916 0.826 
AR 406.12 312.18 546.24 411.58 25.65 24.15 0.928 0.875 
ARXM 402.07 294.15 522.11 380.46 22.99 22.69 0.917 0.868 
ANFIS 391.48 287.24 508.28 373.40 22.98 23.07 0.928 0.887 
2d Raw 919.38 598.66 1,209.10 818.53 23.96 26.86 0.900 0.861 
AR 853.29 584.17 1,117.40 765.04 23.64 23.64 0.913 0.892 
ARXM 847.58 572.37 1,096.22 756.68 22.68 24.36 0.912 0.886 
ANFIS 814.47 549.51 1,044.33 715.91 22.01 23.24 0.930 0.891 
3d Raw 1,365.86 767.42 1,720.63 1,042.29 20.62 26.37 0.911 0.855 
AR 1,293.59 749.53 1,741.54 1,103.87 25.72 32.10 0.930 0.884 
ARXM 1,265.91 766.36 1,644.52 1,059.29 23.02 27.65 0.932 0.891 
ANFIS 1,212.78 763.28 1,542.81 997.71 21.39 23.50 0.930 0.888 
Figure 11

The predictive Q-Q plots and Kolmogorov 5% significance band for lead times of (a) 1d, (b) 2d, and (c) 3d.

Figure 11

The predictive Q-Q plots and Kolmogorov 5% significance band for lead times of (a) 1d, (b) 2d, and (c) 3d.

The evaluation index results of 90% confidence interval at different lead times are illustrated in Figure 12. The average and relative band-widths gradually increase as lead times increase, indicating that the inflow forecast uncertainty increases. The B and RB values of three updating models are smaller than that of the original forecast. The ANFIS updating model has the lowest B and RB values among the three updating models. In terms of CR values, the coverage rate of different updating models in all lead times is close to or exceeds the given confidence level 90%, indicating that the 90% forecast interval is reasonable and reliable. The PUCI indicator is used to compare the performance of different updating models. As shown in Figure 12, the PUCI values of probabilistic forecast decrease as lead times increase, which means that the forecast uncertainty of reservoir inflow increases with lead times. On the other hand, the PUCI value of updated forecast increases significantly compared with that of the raw deterministic forecast, which implies that the probabilistic forecasts are more reliable after error correction. In terms of PUCI value, the ANFIS model performs the best with the largest PUCI value among the three updating models. There is no significant difference between the ARXM and AR models.

Figure 12

Evaluation index of 90% confidence interval for lead times of 1d, 2d, and 3d.

Figure 12

Evaluation index of 90% confidence interval for lead times of 1d, 2d, and 3d.

The 90% confidence intervals of the raw forecast model and ANFIS updating model in the flood season in the validation period are plotted in Figure 13. It is observed that the 90% confidence interval covers most of the observed flows, indicating that the attained probabilistic forecast interval is reliable. The band-width of 90% confidence intervals attained by the updating models is narrower than the original. Furthermore, the 90% confidence intervals become wider at a longer lead time, indicating that the forecast uncertainty becomes larger as lead time increases.

Figure 13

The observed and median forecasted flows and 90% confidence interval for lead times of 1d, 2d, and 3d in the validation period.

Figure 13

The observed and median forecasted flows and 90% confidence interval for lead times of 1d, 2d, and 3d in the validation period.

CONCLUSIONS

Quantifying the inherent uncertainty in flood forecasting is of great importance for flood control decision-making and water resources management. Properly communicating uncertainties in flood forecasting provides the basis for reservoir operation, and thus prevents casualties and economic losses. However, traditional flood forecasting does not consider predictive uncertainty. This study develops a conditional probability model based on copula to quantify and estimate predictive uncertainty. Three updating models, namely the AR model, ARXM model, and ANFIS model, are applied to correct raw forecasts. The TGR in China is selected as a case study. The main conclusions of this study can be summarized as follows:

  • 1.

    Three updating models can significantly improve the raw forecast results in terms of the deterministic evaluation criteria. The degree of improvement in forecast accuracy by the ANFIS updating model outperforms the ARXM and AR models.

  • 2.

    The conditional probability model based on copula is a useful tool for post-processing deterministic forecast to quantify forecast uncertainty, providing both a deterministic forecast value along with a reasonable and reliable uncertainty interval. The updated probability forecasting results show higher accuracy and reliability than that of raw forecasts in terms of CRPS, , and PUCI, indicating that the updating models are a useful tool to improve the reliability of probabilistic forecast. The ANFIS model is superior to the ARXM and AR models in terms of comprehensive indicator CRPS and PUCI in almost all lead times.

  • 3.

    The conditional probability forecast model based on copula can quantitatively assess forecast uncertainty, and the three updating models can improve the accuracy and reliability of the probability forecast. In terms of deterministic flood forecasting, the conditional probabilistic forecasting model based on copula can improve the forecast accuracy; however, the degree of improvement by the three updating models is more significant.

In summary, the developed conditional probability model combined with the updating models quantify the hydrological forecasting uncertainty accurately and reliably. Like most existing studies, this paper focuses on evaluating the pros and cons of probabilistic forecasting schemes; however, how to effectively take advantage of the risk information provided by probabilistic forecasts in the reservoir operation decision-making process needs further study.

ACKNOWLEDGEMENTS

This study is financially supported by the National Key R&D Project (2016YFC0402206) and National Natural Science Foundation (51539009, 51879192) of China. We are very grateful to the editor and anonymous reviewers for their valuable comments and constructive suggestions, and Dr Bryan O'Reilly for proofreading that helped us to greatly improve the manuscript.

REFERENCES

REFERENCES
Box
G. E. P.
&
Cox
D. R.
1964
An analysis of transformations
.
J. R. Stat. Soc.
26
(
2
),
211
252
.
Duan
Q.
,
Ajami
N. K.
,
Gao
X.
&
Sorooshian
S.
2007
Multi-model ensemble hydrologic prediction using Bayesian model averaging
.
Adv. Water Resour.
30
(
5
),
1371
1386
.
Gneiting
T.
&
Raftery
A. E.
2007
Probabilistic forecasts, calibration and sharpness
.
J. Roy. Stat. Soc. Ser. B Stat. Met.
69
(
2
),
243
268
.
Gneiting
T.
,
Raftery
A. E.
,
Westveld
A. H.
&
Goldman
T.
2005
Calibrated probabilistic forecasting using ensemble model output statistics and minimum CRPS estimation
.
Mon. Weather Rev.
133
(
5
),
1098
1118
.
Goldberg
D. E.
1989
Genetic Algorithms in Search, Osptimization and Machine Learning
,
Vol. 13
(
7
).
Addison Wesley
,
Reading
, pp.
2104
2116
.
He
S.
,
Guo
S.
,
Liu
Z.
,
Yin
J.
,
Chen
K.
&
Wu
X.
2018
Uncertainty analysis of hydrological multi-model ensembles based on CBP-BMA method
.
Hydrol. Res.
49
(
5
),
1636
1651
.
Jang
J. S. R.
1993
ANFIS: adaptive-network-based fuzzy inference system
.
IEEE Trans. Syst. Man Cybern.
23
(
3
),
665
685
.
Kelly
K. S.
&
Krzysztofowicz
R.
1997
A bivariate meta-Gaussian density for use in hydrology
.
Stoch. Hydrol. Hydraulics
11
(
1
),
17
31
.
Krzysztofowicz
R.
&
Kelly
K. S.
2000
Hydrologic uncertainty processor for probabilistic river stage forecasting
.
Water Resour. Res.
36
(
11
),
3265
3277
.
Liu
Y.
&
Gupta
H. V.
2007
Uncertainty in hydrologic modeling: toward an integrated data assimilation framework
.
Water Resour. Res.
43
(
7
),
W07401
.
Liu
Z.
,
Guo
S.
,
Xiong
L.
&
Xu
C. Y.
2018
Hydrologic uncertainty processor based on copula function
.
Hydrol. Sci. J.
49
(
3
),
332
342
.
Madadgar
S.
&
Moradkhani
H.
2014
Improved Bayesian multimodeling: integration of copulas and Bayesian model averaging
.
Water Resour. Res.
50
(
12
),
9586
9603
.
Montanari
A.
&
Grossi
G.
2008
Estimating the uncertainty of hydrological forecasts: a statistical approach
.
Water Resour. Res.
44
(
12
),
681
687
.
Nelsen
R. B.
2006
An Introduction to Copulas
,
2nd edn
.
Springer
,
New York
.
Palmer
T. N.
,
Shutts
G. J.
,
Hagedorn
R.
,
Doblas-Reyes
F. J.
,
Jung
T.
&
Leutbecher
M.
2005
Representing model uncertainty in weather and climate prediction
.
Ann. Rev. Earth Planet. Sci.
33
,
163
193
.
Raftery
A. E.
,
Gneiting
T.
,
Balabdaoui
F.
&
Polakowski
M.
2005
Using Bayesian model averaging to calibrate forecast ensembles
.
Mon. Weather Rev.
133
(
5
),
1155
1174
.
Renard
B.
,
Kavetski
D.
,
Kuczera
G.
,
Thyer
M.
&
Franks
S. W.
2010
Understanding predictive uncertainty in hydrologic modeling: the challenge of identifying input and structural errors
.
Water Resour. Res.
46
(
5
),
1187
1191
.
Shen
J. C.
,
Chang
C. H.
,
Wu
S. J.
,
Hsu
C. T.
&
Lien
H. C.
2015
Real-time correction of water stage forecast using combination of forecasted errors by time series models and Kalman filter method
.
Stoch. Environ. Res. Risk A
29
(
7
),
1903
1920
.
Sklar
M.
1959
Fonctions de repartition an dimensions et leurs marges
.
Publ. Inst. Statist. Univ. Paris
8
,
229
231
.
Weerts
A.
,
Winsemius
H.
&
Verkade
J.
2010
Estimation of predictive hydrological uncertainty using quantile regression: examples from the National Flood Forecasting System (England and Wales)
.
Hydrol. Earth Syst. Sci.
19
(
7
),
3181
3201
.
Yin
J.
,
Guo
S.
,
Liu
Z.
,
Yang
G.
,
Zhong
Y.
&
Liu
D.
2018a
Uncertainty analysis of bivariate design flood estimation and its impacts on reservoir routing
.
Water Resour. Manag.
32
(
5
),
1795
1809
.
Yin
J. B.
,
Guo
S. L.
,
He
S. K.
,
Guo
J. L.
,
Hong
X. J.
&
Liu
Z. J.
2018b
A copula-based analysis of projected climate changes to bivariate flood quantiles
.
J. Hydrol.
566
,
23
42
.
Zappa
M.
,
Rotach
M. W.
,
Arpagaus
M.
,
Dorninger
M.
,
Hegg
C.
,
Montani
A.
,
Ranzi
R.
,
Ament
F.
,
Germann
U.
,
Grossi
G.
,
Jaun
S.
,
Rossa
A.
,
Vogt
S.
,
Walser
A.
,
Wehrhan
J.
&
Wunram
C.
2008
MAP DPHASE: real-time demonstration of hydrological ensemble prediction systems
.
Atmos. Sci. Lett.
9
(
2
),
80
87
.
Zhang
L.
&
Singh
V. P.
2006
Bivariate flood frequency analysis using the copula method
.
J. Hydrol. Eng.
11
(
2
),
150
164
.
Zhao
R. J.
1980
The Xinanjiang model.Proceedings Oxford Symposium
, Vol. 129.
IAHS Publ.
pp.
351
356
.
Zhao
T.
,
Wang
Q. J.
,
Bennett
J. C.
,
Robertson
D. E.
,
Shao
Q.
&
Zhao
J.
2015
Quantifying predictive uncertainty of streamflow forecasts based on a Bayesian joint probability model
.
J. Hydrol.
528
,
329
340
.
Zhong
Y.
,
Guo
S.
,
Liu
Z.
,
Wang
Y.
&
Yin
J.
2017
Quantifying differences between reservoir inflows and dam site floods using frequency and risk analysis methods
.
Stoch. Environ. Res. Risk Assess.
32
(
2
),
419
433
.