Hydrological uncertainty processor based on a copula function

ABSTRACT

Quantifying the uncertainty in hydrological forecasting is valuable for water resources management and decision-making processes. The hydrological uncertainty processor (HUP) can quantify hydrological uncertainty and produce probabilistic forecasts under the hypothesis that there is no input uncertainty. This study proposes a HUP based on a copula function, in which the prior density and likelihood function are explicitly expressed, and the posterior density and distribution obtained using Monte Carlo sampling. The copula-based HUP was applied to the Three Gorges Reservoir, and compared with the meta-Gaussian HUP. The Nash-Sutcliffe efficiency and relative error were used as evaluation criteria for deterministic forecasts, while predictive QQ plot, reliability, resolution and continuous rank probability score (CRPS) were used for probabilistic forecasts. The results show that the proposed copula-based HUP is comparable to the meta-Gaussian HUP in terms of the posterior median forecasts, and that its probabilistic forecasts have slightly higher reliability and lower resolution compared to the meta-Gaussian HUP. Based on the CRPS, both HUPs were found superior to deterministic forecasts, highlighting the effectiveness of probabilistic forecasts, with the copula-based HUP marginally better than the meta-Gaussian HUP.

KEYWORDS:

• hydrological forecasting
• uncertainty analysis
• probabilistic forecasts
• copula function
• Three Gorges Reservoir

Editor D. Koutsoyiannis Associate editor E. Gargouri

1 Introduction

Hydrological forecasting is a crucial non-structural flood mitigation measure and provides an essential basis for flood warning, flood control and reservoir operation (Guo et al. 2004, Calvo and Savi 2009, Chen et al. 2014c, Zhang et al. 2015). The forecasting models that are widely used at present are typically deterministic and model outputs are provided to users in the form of deterministic values (Chen and Yu 2007, Coccia and Todini 2011, Ma et al. 2013, Bergstrand et al. 2014, Li et al. 2014). However, a hydrological forecasting model is only a simulation of the real hydrological processes and is therefore imperfect and not precise (Ravines et al. 2008, Wetterhall et al. 2013). These models accept hydrological input, meteorological input, etc., and utilize conceptualized model parameters; these complex factors inevitably cause uncertainties in the hydrological forecasts (Freer et al. 1996, Montanari 2007, Montanari and Grossi 2008, Renard et al. 2010, Chen et al. 2014a). The principle of rational decision making under uncertainty indicates that when a deterministic forecast turns out to be wrong, the consequences will probably be worse than a situation where no forecast is available (Krzysztofowicz 1999, Verkade and Werner 2011, Ramos et al. 2013, Wetterhall et al. 2013). A rational decision maker who wants to make optimal decisions should therefore take forecast uncertainty explicitly into account (Berger 1985, Verkade and Werner 2011, Ramos et al. 2013). Therefore, quantitative assessment of inherent uncertainty is a critical issue. Hydrological forecasting services are trending toward providing users with probabilistic forecasts, in place of traditional deterministic forecasts.

The transition from a deterministic forecast to a probabilistic forecast is based on quantification of the uncertainty inherent in the deterministic forecast. The Bayesian Forecasting System (BFS) proposed by Krzysztofowicz (1999) provides a general framework to produce probabilistic forecasts via any deterministic hydrological model. Various probabilistic forecasting systems suited to different purposes have been developed within this framework (Reggiani and Weerts 2008, Calvo and Savi 2009, Biondi et al. 2010, Weerts et al. 2011, Sikorska et al. 2012, Pokhrel et al. 2013).

In the BFS, the total uncertainty is decomposed into input uncertainty and hydrological uncertainty. The hydrological uncertainty processor (HUP) is a component of the BFS that quantifies the hydrological uncertainty and produces a probabilistic forecast under the hypothesis that there is no input uncertainty (Krzysztofowicz and Kelly 2000). Through Bayes’ theorem, the HUP combines a prior distribution, which describes the natural uncertainty about the realization of a hydrological process, with a likelihood function that quantifies the uncertainty in model forecasts, and outputs a posterior distribution, conditional upon the deterministic forecasts. This posterior distribution provides a complete characterization of uncertainty, including quantiles, prediction intervals and probabilities of exceedence for specified thresholds, which are needed by rational decision makers and information providers who want to extract forecast products for their customers.

The HUP can be implemented in many ways, as different mathematical models for the prior distribution and likelihood function can be developed. Krzysztofowicz and Kelly (2000) introduced a meta-Gaussian HUP, which was developed on the basis of converting both original observations and model forecasts into a Gaussian space by using the normal quantile transform (NQT). This meta-Gaussian HUP has been used widely by many researchers in the fields of hydrology and meteorology (Chen and Yu 2007, Biondi et al. 2010, Biondi and De Luca 2013, Chen et al. 2013a).

Actually, the prior density and likelihood function are conditional probability distributions. It is well known that copula functions have outstanding capability to model joint distributions and give flexibility in choosing an arbitrary marginal distribution (e.g. non-Gaussian form), nonlinear and heteroscedastic dependence structure. The conditional probability distribution can be expressed in explicit form using a copula function (Favre et al. 2004, Nelsen 2006, Zhang and Singh 2006, 2007a, 2007b, 2007c, Genest and Favre 2007, Bárdossy and Li 2008, Chen et al. 2010, Zhang et al. 2011, 2013, 2012). These advantageous characteristics of a copula function motivated us to develop the prior distribution and likelihood function models in the original space directly, without a data transformation procedure into Gaussian space. The main objective of this paper is to propose a post-processor based on a copula function for a deterministic forecast model to produce probabilistic forecasts within the general framework of the HUP.

The remainder of this paper is organized as follows: the methodology is outlined and presented in Section 2; Section 3 introduces a case study of probabilistic inflow forecasts for the Three Gorges Reservoir (TGR) in China; the results and discussion are included in Section 4; and Section 5 draws the main conclusions of the work.

2 Methodologies

2.1 Hydrological uncertainty processor

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 H₀ represent the observed discharge at the time n = 0 when the forecast is prepared; then H_n (n = 1, 2, …, N) is the observed discharge at lead time n, and S_n (n = 1, 2, …, N) is the corresponding deterministic forecast discharge at lead time n. What the rational decision maker then needs is not a single number s_n, but the distribution function of predictand H_n, conditional on H₀ = h₀ and S_n = s_n. The purpose of the HUP is to supply such a conditional distribution function through Bayesian revision (Liu et al. 2016).

The posterior density function of predictand H_n, conditional on S_n = s_n and H₀ = h₀, is derived via Bayes’ theorem (Krzysztofowicz and Kelly 2000):

(1)

In concept, Bayes’ theorem revises the prior density function, which characterizes the prior uncertainty about H_n given H₀ = h₀. The extent of the revision is determined by the likelihood function, which characterizes the degree to which S_n = s_n reduces the uncertainty about H_n. The result of this revision is the posterior density function , which quantifies the uncertainty about H_n that remains after the deterministic forecast model generates forecast S_n = s_n.

2.2 Meta-Gaussian HUP

The structure of Equation (1) indicates that the posterior density strongly relies on both the prior density function and the likelihood function. At present, the most commonly used technique to construct the prior density and likelihood functions is the meta-Gaussian model, which involves three steps. First, both actual flow H_n and predicted flow S_n are converted into Gaussian space by employing the NQT method (Bogner et al. 2012). Second, the transformed and series are assumed to be linear and normally distributed. Finally, the posterior density function of H_n in the transformed Gaussian space is derived by linear regression, from which the posterior density of H_n in the original space can be found based on the Jacobian inverse transformation. For more detailed descriptions and mathematical formulations of the meta-Gaussian HUP, readers are referred to Krzysztofowicz and Kelly (2000) or Liu et al. (2016).

2.3 Copula-based HUP

Copula functions are an effective tool used to develop prior distribution and likelihood function models, in which the predictand and the deterministic forecast are allowed to have distribution functions of any form, along with nonlinear and heteroscedastic dependence structure. Therefore, they can be implemented in the original space directly without a data transformation procedure into Gaussian space.

2.3.1 Copula theory

A copula function connects multivariate probability distributions to their one-dimensional marginal distributions (Nelsen 2006). Let () be the cumulative distribution function (cdf) of X_i. The multivariate distribution can be expressed in terms of its marginal and the associated dependence function using Sklar’s theorem:

(2)

where C, the copula function, captures the essential features of the dependence among the random variables. A more detailed description of copula function properties can be found in Favre et al. (2004), Nelsen (2006) and Zhang et al. (2013).

2.3.2 Prior density

The prior cdf of H_n given H₀ = h₀ can be expressed as:

(3)

where is the conditional cdf, and P is the non-exceedence probability.

The prior density function is the corresponding probability density function (pdf) of and can be defined as:

(4)

Let H₀ and H_n be random variables with marginal cdfs, and . Then, U₁ and U₂ are uniformly distributed random variables; u₁ denotes a specific value of U₁, and u₂ denotes a specific value of U₂. Using the copula function, the joint cdf is expressed by .

Using the copula function, the conditional cdf and pdf can be rewritten as (Zhang and Singh 2006):

(5)

(6)

where is the density function of , and ; is the pdf of H_n. Equation (6) is the expression of the prior pdf.

2.3.3 Likelihood function

It is considered that S_n is a random variable with marginal cdf and pdf . The conditional cdf of S_n given H₀ = h₀ and H_n = h_n can be expressed as:

(7)

where is the conditional cdf.

The corresponding pdf of is defined as:

(8)

Using the copula function, the joint cdf of H₀, H_n and S_n, denoted as, can be expressed as . So the conditional cdf and pdf are rewritten as (Zhang and Singh 2007c):

(9)

(10)

where is the density function of . From another point of view, given H₀ = h₀ and S_n = s_n, the likelihood function of H_n can be calculated using Equation (10).

2.3.4 Posterior density

Substituting the Equations (6) and (10) into Equation (1), the pdf of H_n can be rewritten as:

(11)

For fixed realizations H₀ = h₀ and S_n = s_n, u₁ and u₃ are constants, while u₂ varies from 0 to 1. Since the denominator cannot be obtained directly by an analytic method, the Monte Carlo sampling technique (Xiong et al. 2014, Yu et al. 2014) is applied by the following steps:

1. Generate M random numbers u₂ from uniform distribution U(0, 1);

2. Compute the value of ;

3. Calculate the mean of the M calculated values, which is approximately equal to the definite integral (Kroese et al. 2013, Robert and Casella 2013).

Subsequently, the pdf can also be estimated.

2.4 Candidate marginal distributions and trivariate copulas

Specifying and determining marginal distributions of the actual flow H₀, H_n and predicted flow S_n (n = 1, 2, …, N) is an essential step. The actual flows {H_n: n = 0, 1, …, N} are considered as random variables. Given only such a record, there is usually no basis for assigning a probability distribution to flow H_n that differs from the distribution assigned to flow H₀, for any n = 1, 2, …, N within a few days. In other words, there is no a statistical difference between these 1 + N flow series (Koutsoyiannis and Montanari 2015). Therefore, we hold the opinion that the variables H_n follow the same marginal cumulative distribution functions (cdf) with H₀, and thus only the cdf of H₀ needed to be fitted. The predicted flows {S_n: n = 1, …, N} are considered as different random variables and different cdfs needed to be fitted for variable S_n.

The main purpose of this study is not hydrological frequency analysis, which aims to extrapolate the extreme events far beyond the observations (very low probabilities, e.g. return period exceeds 1000 years). Instead, the probability distribution of daily flows used in this paper refers to the flow–duration curve, which gives a summary of flow variability at a site and is interpreted as a relationship between any discharge value and the percentage of time that this discharge is equalled or exceeded during a given period (Vogel and Fennessey 1994, Castellarin et al. 2004, Shao et al. 2009). The flow–duration curve has been widely used by engineers and hydrologists around the world in numerous applications, such as hydropower generation, inflow forecasting, and designing of irrigation systems (Vogel and Fennessey 1995, Yokoo and Sivapalan 2011, Gottschalk et al. 2013).

In our study, we concentrated on the daily flow–duration curve. If the daily streamflow is assumed to be a random variable, the flow–duration curve may also be viewed as the complement of the cumulative distribution function used in hydrological frequency analysis when identifying the percentage of time with probability (Castellarin et al. 2004). As a consequence, the flow–duration curve is also a very practical tool used to describe hydrological regimes and represents the relationship between magnitude and frequency of flow (Vogel and Fennessey 1995, Liucci et al. 2014, Xiong et al. 2015).

Six commonly used distributions in hydrology, namely normal, gamma, Gumbel, Pearson type III, log-normal and log-Weibull, were selected as candidate models for H₀ and S_n (n = 1, 2, …, N). These univariate probability distributions are summarized in Table 1. In this study, the distribution parameters were estimated by L-moments method (Hosking 1990, Liu et al. 2017). The Kolmogorov-Smirnov (K-S) test statistic, D, was adopted to measure the goodness of fit between the candidate theoretical distribution and the empirical distribution (Tsai et al. 2001, Rahman et al. 2010). The 5% significance level was selected to reject or accept a fitted distribution, and then the probability distribution that provides the minimum D value was chosen as the best fitting distribution.

Table 1. Candidate univariate distributions adopted in this study.

Distribution	Probability distribution function (pdf)	Range
Normal
Gamma
Gumbel
Pearson type III
Log-normal	,
Log-Weibull

Different families of copulas have been described by Nelsen (2006). The Archimedean copula family is the most desirable for hydrological analyses, because it can easily be constructed and it has been applied frequently by many authors (Favre et al. 2004, Zhang and Singh 2006, 2007a, 2007b, 2007c, Salvadori and De Michele 2007, Chebana and Quarda 2011, Li et al. 2013, Chen et al. 2013b, 2013c, 2014b). To estimate the pdfs expressed in Equation (11), three-dimensional (3D) joint distributions of H₀, H_n and S_n need to be constructed. In this study, the symmetric copulas were not considered because the dependence among the three variable pairs (H₀, H_n), (H₀, S_n) and (H_n, S_n) is not the same, which will be tested against data for the case study. Hence, we used three widely used asymmetric trivariate Archimedean copulas, namely Gumbel-Hougaard, Frank and Clayton as candidates. These three trivariate Archimedean copulas are described in Table 2. Dependence parameters of the trivariate copula functions were estimated using the maximum pseudo-likelihood method (Zhang and Singh 2007b, 2007c, Chen et al. 2010). The root mean square error (RMSE) was used to measure the goodness of fit of the copula distribution (Zhang and Singh 2007a). The copula that has the smallest RMSE value is preferred.

Table 2. Summary of the three candidate trivariate Archimedean copulas.

Copula
Gumbel-Hougaard
Frank
Clayton

2.5 Evaluation criteria used

2.5.1 Performance of deterministic forecasts

Two widely applied criteria, namely the Nash-Sutcliffe efficiency (NSE) and relative error (RE), were adopted to evaluate the performance of the deterministic forecast model (Nash and Sutcliffe 1970, Xiong and Guo 1999, Akhtar et al. 2008). The mathematical expressions of NSE and RE can be found in these references. An NSE value close to 1 and RE value close to zero indicate a good simulation (Liu et al. 2010, Ma et al. 2010).

2.5.2 Performances of probabilistic forecasts

Several methods, e.g. predictive quantile–quantile (QQ) plot, α-index and π-index have been proposed in the literature to evaluate probabilistic forecasts (see e.g. Gneiting et al. 2007, Laio and Tamea 2007, Thyer et al. 2009, Engeland et al. 2010, Renard et al. 2010, Evin et al. 2014, Madadgar et al. 2014, Smith et al. 2015) and were used in this study. The interpretation of the predictive QQ plot can be found in Thyer et al. (2009). For the mathematical expressions of α-index and π-index, please refer to the previous references. The value of α-index varies between 0 (worst reliability) and 1 (perfect reliability). A greater value of π-index indicates greater resolution (lower uncertainty) of forecasts.

The trade-off between reliability and sharpness has been discussed in previous research (Xiong et al. 2009, Li et al. 2010a, Kasiviswanathan et al. 2013), which showed that these two desirable objectives could not be achieved simultaneously. The continuous rank probability score (CRPS) is a standard measure that combines reliability and sharpness (Hersbach 2000, Gneiting et al. 2005) and is used for selecting the preferred model. The CRPS is defined as (Hersbach 2000, Pappenberger et al. 2015):

(12)

where denotes the Heaviside step function and takes the value 0 when and 1 otherwise. For a deterministic forecast system, the CRPS reduces to the mean absolute error (MAE). This is an advantage of CRPS and consequently allows the comparison of deterministic and probabilistic forecasts (Pappenberger et al. 2015, Zhao et al. 2015). The smaller the CRPS value, the better the prediction performance.

3 Case studies

3.1 Study area

The Yangtze River is one of the largest rivers in the world, and the Three Gorges Reservoir (TGR) is a vitally important, backbone project in the development and harnessing of the Yangtze River in China, as shown in Figure 1. The basin area of the TGR is 1 000 000 km², and the Yichang hydrological station is the control station for the TGR. The annual average discharge and runoff volumes at the dam site are 14 300 m³/s and 4510 × 10⁸ m³, respectively. The total storage capacity of the TGR is 393 × 10⁸ m³, of which 221.5 × 10⁸ m³ is flood control storage. Since the primary role of the TGR is flood prevention for the downstream area of the Yangtze River basin, flood forecasting is of great importance for reservoir operation, flood prevention and disaster relief (Li et al. 2013).

Figure 1. Sketch map of the Three Gorges Reservoir and Upper Yangtze River basin.

3.2 Deterministic inflow forecasts of the TGR

As shown in Figure 1, the inflow of the TGR has three components: the main upstream inflow (Cuntan station), the tributary inflow from the Wu River (Wulong station), and the lateral flow from the TGR intervening basin calculated by the Xinanjiang model (Zhao 1992, Cheng et al. 2006, Lin et al. 2014, Si et al. 2015). The multiple-input single-output linear systematic model was adopted to simulate the inflow of the TGR (Liang et al. 1992); for more details please refer to Li et al. (2010b) and Chen et al. (2015).

The whole data set ranged from 2003 to 2009, of which the period 2003–2007 was used for model calibration and 2008–2009 was used for validation. The simulation results for the NSE and RE in the calibration period are 97.72% and −1.04%, respectively, and in the verification period 95.84% and −0.21%, respectively. These results show that the deterministic forecast model is proven to be quite efficient in simulating the inflow series for the TGR.

4 Results and discussion

The forecasting in this study can be viewed as hindcasting essentially, in which the future rainfall is treated as perfect, and thus focuses on quantifying the hydrological uncertainty. The lead times are 24 h (n = 1), 48 h (n = 2) and 72 h (n = 3). The observed rainfalls are treated as the “perfect rainfall forecasts”, and are used to produce deterministic model inflows (s₁, s₂, s₃). They are attached to actual inflows (h₀, h₁, h₂, h₃) to obtain one joint realization of the model–actual inflow process. The dataset from 2003 to 2009 was used to calibrate and compare the meta-Gaussian HUP and copula-based HUP.

4.1 Determination of marginal distributions

The sample series of H₀ was taken from 1 June to 27 September every year, S₁ from 2 June to 28 September, S₂ from 3 June to 29 September, and S₃ from 4 June to 30 September; thus, all these four variables have a data length of 833. The parameters of the six candidate distributions were estimated by the L-moment method and the results are listed in Table 3. The K-S test was used to verify the null hypothesis and the corresponding K-S statistic (D) values are reported in Table 4, which shows that the null hypothesis could not be rejected at the 5% significance level (critical value: 0.0471) for all six candidate distributions except the normal distribution. For the four hydrological variables, the Gumbel distribution provided the minimum D value and was chosen as the best fitting distribution. Figure 2 shows the empirical cdf values obtained from the Gringorten plotting-position formula (Zhang and Singh 2006) and theoretical cdf values calculated by the Gumbel distributions. It can be seen that the theoretical values fit the empirical values very well. For comparison purposes, the copula-based HUP used the same marginal distributions as the meta-Gaussian HUP in this paper.

Table 3. Parameter estimation results for four marginal distributions.

Variable	Normal		Gamma		Gumbel		Pearson type III			Log-normal		Log-Weibull
H0	22 354.30	8454.22	6.74	3318.06	6881.34	18 382.28	4.08	0.00023	4741.82	9.94	0.39	9.79	19.48	6109
S1	22 314.90	8393.89	6.81	3275.31	6832.24	18 371.22	3.92	0.00022	4639.73	9.94	0.39	9.79	19.52	6132
S2	22 240.16	8140.83	7.21	3084.93	6626.25	18 415.39	4.43	0.00024	3979.82	9.94	0.38	9.80	20.36	5991
S3	21 897.76	7763.05	7.70	2842.82	6318.76	18 250.47	5.95	0.00031	2670.46	9.93	0.37	9.78	20.88	5985

Table 4. Results of Kolmogorov-Smirnov test statistic D for four marginal distributions.

Variable	Normal	Gamma	Gumbel	Pearson type III	Log-normal	Log-Weibull
H0	0.0780	0.0306	0.0172	0.0219	0.0219	0.0244
S1	0.0780	0.0266	0.0183	0.0246	0.0246	0.0207
S2	0.0717	0.0467	0.0255	0.0264	0.0264	0.0263
S3	0.0720	0.0467	0.0178	0.0181	0.0181	0.0270

Figure 2. Empirical and theoretical values fitted by Gumbel distributions.

4.2 Calibration of copula-based HUP

The rank-based correlation (Kendall’s coefficient) matrix of variables H₀, H_n and S_n is shown in Table 5. It is demonstrated that the dependence among the three variable pairs (H₀, H_n), (H₀, S_n) and (H_n, S_n) is not the same. Furthermore, the highest correlation coefficient is exhibited in the variable pair (H_n, S_n). This result indicates that the asymmetric trivariate copula functions may be more appropriate than symmetric ones for use in the 3D joint distributions of H₀, H_n and S_n. When constructing the 3D joint distributions using the asymmetric copula functions, the structures (H_n, S_n)H₀ were applied. Specifically, the copula was first built for (H_n, S_n) and then for H₀ – .

Table 5. Ranked-based correlation matrix of the variables.

Lead times	Variables
24 h	H0, H1	0.823
	H0, S1	0.824
	H1, S1	0.929
48 h	H0, H2	0.694
	H0, S2	0.711
	H2, S2	0.883
72 h	H0, H3	0.600
	H0, S3	0.658
	H3, S3	0.828

The 3D joint distributions of H₀, H_n and S_n (n = 1, 2, 3) were constructed using the three candidate trivariate copula functions. Dependence parameters of the trivariate copula functions were estimated using the maximum pseudo-likelihood method and the results are listed in Table 6. It was found that the Frank copula performed best, with the smallest RMSE values for the three joint distributions. Empirical cdfs obtained from the Gringorten plotting-position formula and theoretical cdfs calculated from the Frank copula for the three joint distributions are plotted in Figure 3. An overall satisfactory agreement between the empirical and theoretical cdfs is shown. Hence, the asymmetric trivariate Frank copula functions have good performance in modelling the joint distributions of H₀, H_n and S_n.

Table 6. Estimated parameters of the three candidate copulas.

Variables	Gumbel-Hougaard		Frank		Clayton
	[θ1, θ2]	RMSE	[θ1, θ2]	RMSE	[θ1, θ2]	RMSE
H0, H1, S1	[9.08, 10.65]	0.0116	[25.44, 35.23]	0.0103	[15.16, 20.78]	0.0150
H0, H2, S2	[4.25, 7.58]	0.0141	[13.45, 20.38]	0.0112	[6.49, 14.33]	0.0186
H0, H3, S3	[2.98, 6.57]	0.0144	[9.58, 16.21]	0.0117	[4.58, 9.34]	0.0188

Figure 3. Plots of empirical and theoretical values estimated by Frank copulas for three joint cdfs. (Note: rank represents number of ordered pair, ranked in ascending order in terms of theoretical joint cdf, respectively.].

4.3 Comparison of the meta-Gaussian HUP and copula-based HUP

4.3.1 Posterior median forecasts

For 24-, 48- and 72-h lead times, the NSE and RE calculated by both the deterministic forecast model and posterior median forecasting associated with the meta-Gaussian HUP and copula-based HUP are listed in Table 7. It is shown that the results of both the meta-Gaussian HUP and the copula-based HUP are slightly better than those of the deterministic forecast model, and the copula-based HUP is comparable to the meta-Gaussian HUP. Compared with deterministic forecasts, the NSE and the RE of the copula-based HUP for 24-, 48- and 72-h lead time forecasts are improved by 1.24, 1.26 and 1.26% and reduced by 0.17, 0.57 and 1.72%, respectively. It is also noted that the accuracy of posterior median forecasts of the both HUPs decreases as the lead time increases.

Table 7. Comparison of performance evaluation criteria for deterministic forecasts.

Lead time	Deterministic model		Meta-Gaussian HUP		Copula-based HUP
	NSE (%)	RE (%)	NSE (%)	RE (%)	NSE (%)	RE (%)
24 h	97.55	–0.34	97.65	–0.31	98.79	–0.17
48 h	94.10	–0.85	94.54	–0.68	95.36	–0.28
72 h	88.52	–2.51	89.14	–1.14	89.78	–0.79

4.3.2 Probabilistic forecasts

The predictive QQ plot, α-index, π-index and CRPS were adopted to evaluate the probabilistic forecasts. Figure 4 presents the predictive QQ plots with respect to the meta-Gaussian HUP and copula-based HUP for 24-, 48- and 72-h lead times. Using the guide presented in Thyer et al. (2009) to assess the results, it is clear that the overall performances of all predictive QQ plots are acceptable. Both the meta-Gaussian HUP and the copula-based HUP systematically under-predicted the inflows, since the observed p values at the theoretical median are a bit higher than the theoretical quantiles. In addition, the plots also show that the observed p values cluster around the tails (i.e. a high slope around theoretical quantile 0.4–0.6). This finding means that the predictive uncertainty is somewhat underestimated for both HUPs. The overall behaviours of the meta-Gaussian HUP and the copula-based HUP are found to be similar. The QQ plot for the copula-based HUP is slightly closer to the 1:1 line than for the meta-Gaussian HUP. That is to say, the copula-based HUP performs marginally better in terms of reliability. Nonetheless, these underestimations for both meta-Gaussian and copula-based HUPs are in zones where p values are relatively high, indicating that such differences may not be statistically significant.

Figure 4. Predictive QQ plots of meta-Gaussian HUP and copula-based HUP.

The results of α-index, π-index and CRPS are summarized in Table 8. For both the meta-Gaussian HUP and the copula-based HUP, it is clearly shown that the α-index value increases (higher reliability) when the lead time increases. However, it should be noted that this is at the expense of decreasing π-index values (lower resolution). Besides, the copula-based HUP has slightly larger α-index values, but smaller π-index values compared with the meta-Gaussian HUP. In terms of CRPS value, both HUPs outperform the deterministic forecasts, which demonstrates the effectiveness of probabilistic forecasts. Comparison results also indicate that the copula-based HUP is marginally better than the meta-Gaussian HUP. The CRPS value of the copula-based HUP for 24-, 48- and 72-h lead times is improved (decreased) by 16.6, 21.2 and 23.3%, respectively.

Table 8. Comparison of performance evaluation criteria for probabilistic forecasts.

Lead time	Deterministic model	Meta-Gaussian HUP			Copula-based HUP
	CRPS/MAE	α-index	π-index	CRPS	α-index	π-index	CRPS
24 h	688	0.8028	18.55	608	0.8507	16.39	574
48 h	1180	0.8555	12.45	975	0.8916	10.42	930
72 h	1763	0.8879	9.13	1420	0.9184	7.68	1353

Note: α-index and π-index are dimensionless; the unit of CRPS is m³/s.

Although such marginally better performance does not result for each year, for illustrative purposes, the observed and median discharges, and 90% inflow prediction intervals estimated by the meta-Gaussian HUP and copula-based HUP in 2004 are presented in Figures 5 and 6, respectively. It can be seen that most observed inflows are contained within the 90% prediction intervals. This demonstrates that these 90% prediction intervals can effectively capture the forecast uncertainty and provide more information for decision making in flood control and reservoir operation. As lead time increases, the 90% prediction intervals become wider (i.e. greater uncertainty).

Figure 5. Meta-Gaussian HUP: 90% prediction intervals, median and observed discharges in 2004.

Figure 6. Copula-based HUP: 90% prediction intervals, median and observed discharges in 2004.

5 Summary and conclusions

Adequate assessment of uncertainty for flood forecasting is an important issue in hydrological research. Flood forecasting services are trending toward providing users with probabilistic forecasts, in place of traditional deterministic forecasts. This study has developed a copula-based HUP for probabilistic forecasting. The Three Gorges Reservoir (TGR) in China was selected as the case study and the results were compared with those of the widely used meta-Gaussian HUP. The main conclusions are summarized as follows:

1. The output of the HUP is a posterior distribution of the process, conditional upon the deterministic forecast. This posterior distribution provides a complete characterization of uncertainty, including quantiles with specified exceedence probabilities, prediction intervals with specified inclusion probabilities, and probabilities of exceedence for specified thresholds, which are needed by rational decision makers and information providers who want to extract forecast products for their customers.

2. Based on the copula function, the prior density and likelihood function of the HUP are explicitly expressed, and the corresponding posterior density and distribution can be obtained using the Monte Carlo sampling technique. This copula-based HUP can be implemented in the original space directly without a data transformation procedure into Gaussian space and allows for any form of marginal distribution of predictand and the deterministic forecast variable, and a nonlinear and heteroscedastic dependence structure.

3. Comparison results for the case study in China’s TGR demonstrate that the proposed copula-based HUP is comparable to the meta-Gaussian HUP in terms of the posterior median forecasts. It is also shown that probabilistic forecasts produced by the copula-based HUP have slightly higher reliability and lower resolution compared with the meta-Gaussian HUP. On the basis of the CRPS values, it is found that both HUPs are superior to deterministic forecasts, which highlights the effectiveness of probabilistic forecasts, and the copula-based HUP is marginally better than the meta-Gaussian HUP.

It is noted that the HUP discussed in this paper is for a univariate probabilistic forecast. This system independently provides a marginal probability distribution of observed discharge H_n (n = 1, 2, …, N) for each lead time; it does not consider and characterize the stochastic dependence among variables (H₁, H₂, …, H_N). A multivariate probabilistic forecast system (Krzysztofowicz and Maranzano 2004, Engeland and Steinsland 2014) would offer a joint probability distribution of (H₁, H₂, …, H_N) which would require a multivariate HUP and will be considered in further study.

Acknowledgements

The authors would like to thank the Hydrology Bureau of Changjiang Water Resources Commission and the China Three Gorges Corporation for providing the data used in this study.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This study was financially supported by the National Natural Science Foundation of China [NSFC 51539009, 51525902] and the National Key Research and Development Plan of China [2016YFC0402206].