Uncertainty Analysis of Multi-Model Flood Forecasts

This paper demonstrates, by means of a systematic uncertainty analysis, that the use of outputs from more than one model can significantly improve conditional forecasts of discharges or water stages, provided the models are structurally different. Discharge forecasts from two models and the actual forecasted discharge are assumed to form a three-dimensional joint probability density distribution (jpdf), calibrated on long time series of data. The jpdf is decomposed into conditional probability density distributions (cpdf) by means of Bayes formula, as suggested and explored by Krzysztofowicz in a series of papers. In this paper his approach is simplified to optimize conditional forecasts for any set of two forecast models. Its application is demonstrated by means of models developed in a study of flood forecasting for station Stung Treng on the middle reach of the Mekong River in South-East Asia. Four different forecast models were used and pairwise combined: forecast with no model, with persistence model, with a regression model, and with a rainfall-runoff model. Working with cpdfs requires determination of dependency among variables, for which linear regressions are required, as was done by Krzysztofowicz. His Bayesian approach based on transforming observed probability distributions of discharges and forecasts into normal distributions is also explored. Results obtained with his method for normal prior and likelihood distributions are identical to results from direct multiple regressions. Furthermore, it is shown that in the present case forecast accuracy is only marginally improved, if Weibull distributed basic data were converted into normally distributed variables.


Purpose of Study
Among the most important tasks of applied hydrology is the forecasting of river stages or discharges, for which hydrologists and meteorologists have developed increasingly more complex models.The problem considered here is how to obtain the best possible forecast by using outputs from a given set of two calibrated forecast models.The best forecast is defined as that linear combination of model outputs from the set which minimizes the forecast error.More specific, our objective is to provide a systematic error analysis based on Bayes formula in order to obtain the optimum forecast from up to two models, and to determine the probability density function (pdf) of forecast errors eF(i+m).The result is a simplified version of the Bayesian framework developed by Krzysztofowicz [1], and Todini [2].This approach is generally applicable to any two forecast models.For illustration purposes, it is applied to data and models from a study of forecasting discharges for station Stung Treng on Mekong River in South-East Asia, which is described in Shahzad and Plate [3].

Background
Setting up a multi-model forecast system requires a sequence of steps.First step is the selection of the set of models to be used for forecasting.The second step is model calibration because hydrological models generally depend on empirical parameters which are specific for the basin considered.For each of the flood forecast models that set of model parameters is determined which minimizes the error between forecast and observation of long historical time series.The third step is application, i.e., the use of each model for real time forecasting.Finally, in a fourth and final step, the models are suitably combined to yield the optimum real time forecast.Due to model-, parameter-, and other uncertainties only estimates for the true record are obtained, and the purpose of the analysis in this paper is to minimize total uncertainty, which shall be expressed through the error of observed record minus forecast.We shall distinguish between design and operational (conditional) uncertainties, which are of a very different nature [4], as will be discussed briefly in Section 1.2.3.

Model Selection
Hydrological models for river discharges can only approximately reproduce the true hydrological situation of real hydrological basins, because they generally require to combine runoff increments from many poorly defined, aggregated and connected hydrological area elements-local contributions which are generated by climate and precipitation inputs of uncertain magnitude and poorly known spatial distribution acting on complex local hydrological processes.Therefore, forecast models for river discharges (or stages) can only yield estimates, with an uncertainty band which reflects the total uncertainty resulting from data input, model structure, and model parameters.
For a flood forecaster, the physical correctness of a model is of secondary importance as compared to the goodness of fit of its outputs.In principle, any data-based hydrological model for flood calculations could be used."Data based" implies that for calibration sufficiently long records of observations are available to obtain statistically significant outputs from the model.There exist a large number of potentially useful models for this purpose (as summarized by Singh and Woolhiser [5], or see the extensive list in Gouweleeuw et al. [6]).They are basically either time series models, where the output is determined from known input data by means of transition probabilities (time series models), sometimes as artificial neural networks (ANN) [7] or conceptual rainfall-runoff models, usually based on applying linear systems theory (one or more Nash Cascades, see [8] for an extensive discussion).A strong case can be made to not rely on a given model but to develop models based on basic hydrological principles, which optimally account for physical and climatic conditions of the region [9][10][11], perhaps starting with an elementary model, which may be upgraded with increased observational evidence (model development along the "axis of complexity" [12]).

Model Calibration
At the outset, the modeller has to decide which model to use, a process which is guided, among other criteria, by the available data basis, and by the purpose for which the forecast is to be made-agriculture (ecological flooding), navigation, or flood protection are typical topics.Important is selection of time intervals, ∆ t , for example days, and of maximum forecast time, m∆t .It is assumed in this paper that forecast specifications call for models to forecast the complete time series of discharges for m = 5 (days) for one season.
The selected forecast model is calibrated by means of a given set of time series of discharges and, for rainfall-runoff models, precipitation from measuring stations located in or near the river basin.Generally, the whole time series of available data of discharge and rainfall fields should be used.However, in our case the analysis is restricted to the monsoon season of South-East Asia, from June to November of each year.The model is calibrated to determine the set of parameters that finds Water 2015, 7, 6788-6809 the optimum of some objective function, which should depend on the planned use of the forecast.For example, weighted errors are useful for selecting the best among different models, obtained by introducing appropriate objective functions [13,14], which may be based on economic criteria [14], or on risk for people living along the river.However, standard procedure, adapted also in this paper, is to find that set of parameters for the hydrological model which minimizes linear least squares of system errors eM k (i+m) = Q(i+m) ´QF k (i+m).Here, Q(i+m) is the historical discharge time series, observed at time i for m time intervals ahead, and QF k (i+m) is the discharge calculated from forecast model k.Parameters are optimized by means of split sampling techniques, yielding a validated set of optimum parameters for the selected model.
Methods of parameter estimation differ mainly in how parameter spaces are defined and which objective functions are used for determining goodness of fit [15].Because there is generally no physical reason for assigning particular numerical values to most model parameters, different parameter sets may fit the model equally well.This problem of "equi-finality" [16] is inherent in hydrological models.One way to proceed is to define a range in which parameters may vary, assign suitable empirical probability distributions to each of the parameters in the parameter space and determine possible parameters through a Monte Carlo procedure, with different sets used to generate an output pdf.This approach is used for their GLUE method by Beven and Binley [17], Freer et al. [18], or for the method of Gupta et al. [19].Other authors transfer parameters from other catchments, for example Hundecha and Bardossy [20].If the data series is long enough, some final best fit set of parameters for the model can be estimated, whereas for shorter series, a preliminary set obtained from the available data may be upgraded if additional information becomes available, for example by means of Bayesian methods [21].
System error eM k (i+m) sets the limit of forecast quality obtainable with a given model k, as it corresponds to the best forecast obtainable from the model.Further improvement can be obtained only if model structure and/or calibration data base, or both, are improved.However, the forecast error may be reduced if more than one model is used, which can be arranged in series or in parallel.For example, models are used in series, where a decision is made to switch from one model to another depending on season or catchment characteristics [22].Or the models are arranged in parallel, applied to the same period of forecast, and separately calibrated on the same basic data, as in [1][2][3] and used in this paper.For parallel models, optimum forecasts are obtained by combining outputs from all models into a final forecast.The simplest way of using more than one model is to average the forecasts from all models.A direct improvement on this approach is obtained by multiple linear regression of model outputs, giving weights to the individual outputs by means of least squares analysis [1,2,23].In the remainder of this paper the Bayesian version of this approach will be discussed.

Model Application for Discharge Forecasting
For the operational stage, models are assumed given and calibrated.They become a part of the forecasting system, usually embedded in a flood forecasting package, such as by Werner et al. [24].When making a forecast new errors occur in addition to system errors yielding estimates QF k (i+m) for Q(i+m).In the forecasting case, Q(i+m) is the unknown predictand, and QF k (i+m), the output of forecast model k, is the calculated predictor.The difference eF k (i+m) = Q(i+m) ´QF k (i+m) is the forecast error.In contrast to system error eM k (i+m), which is based on average performance with known input and output, forecast error eF k (i+m) is conditional on known initial conditions but depends on unknown forecasts of future inputs [4].For time series models without rainfall input, eF k (i+m) and eM k (i+m) are basically the same, because dependence on basic data is the same in model building and forecasting.For rainfall-runoff models neither rainfall nor lateral inflows from tributaries of the future are known and therefore must be forecast by means of a rainfall forecast model.This yields additional uncertainties.Consequently, forecast errors eF k (i+m) of conditional forecasts from model k are generally larger than system errors eM k (i+m) of the same model.
Water 2015, 7, 6788-6809 Forecast quality of conceptual forecast rainfall-runoff models depends to a large extent on adequacy of the sub-models used for generating rainfall inputs.Many sub-models of different degrees of complexity exist.Simple models use present day rainfall, or rainfall averaged over a number of days [25], as input into the model.Other methods of rainfall forecasting based on more refined time series models have been suggested, [25,26].These methods are empirical, estimating future rainfall records on the basis of past measurements, some by means of different types of artificial neural networks (ANN) [26][27][28][29][30]. Physically based rainfall forecasts are the alternative, which extend existing weather pattern into the future [31], or forecasts using downscaling from large scale weather models, such as the global gridded weather model GME [32].Examples are the sub-grid models COSMO [33], or COSMO-LEPS [33,34] developed for Central Europe.

Performance Indicators
Model calibration and application require different performance indices for skill evaluation.Generally, in both cases at least two quantities [35] should be used: a measure of forecast error magnitude, and a measure of forecast quality.The former measure gives an indication of the absolute magnitude of performance errors, whereas the latter is a measure of comparison of the forecast against that from a "benchmark" model [36].Benchmark models should be models which can be used without using forecast information, such as the simple persistence model, or use of averages for each day of the year taken over many years.
During calibration, indicators for error magnitude should be based on the objective function used for model fitting.Typically measures are rms values of error eM k (i+m), or similar absolute quantities, as summarized by Gupta et al. [19] or Dawson et al. [29], such as mean absolute error values.For forecasting, the criterion for error magnitude requires use of conditional forecasting errors eF k (i+m) [4] instead of eM k (i+m), which occur when a model k is applied to generate forecasts.
As measures of calibration quality the error magnitude is compared with the error of a benchmark model.Nash-Sutcliffe Index (NSI) has been used extensively, i.e., [24,37] among many others, which compares the error variance with the variance of the total record.NSI is useful for comparison purposes, for example to facilitate the selection of the best among different models, or as composite index for many different performance measures Gupta et al. [38].However, it is not very appropriate for forecast performance evaluation.Evaluation of forecast application skills requires use of conditional forecasting errors eF k (i+m) [4], which occur when a model k is applied to generate forecasts.A better reference is persistence index (PI) of Kitanides and Bras [39], which has been used in recent studies, for example [4,7].It compares the forecast error variance with that of the error from assuming persistence as benchmark condition, when the discharge of today is used as predictor for the discharge of the future.Other quality parameters may be based on binary information-correct or false forecast of flood level exceedance, or the performance scores used in meteorology [40,41].
In conclusion, forecasters have to face the problem of having to infer useful forecasts from an unreliable information base-a task that can be realized only if it is acknowledged that one has to live with uncertainty.To compensate uncertainty, in addition to using forecast models, forecasters tend to rely on intuition and experience to adjust forecasts [42].On a more rational basis conditional forecasts should be cast into a probabilistic framework, by means of a probability density function, with the mean as a crisp forecast, and different error bounds as indications of degree of uncertainty.These quantities can be improved by means of efficient use of more than one model.It is shown in this paper that Bayes formula, as suggested by Krzysztofowicz [1], allows to combine the outputs from several models in a meaningful way (see also [2]).For this purpose, a system for classifying forecast errors is developed in Section 2, and in Section 3 these error classifications are applied to outputs of simple forecast models described in Shahzad and Plate [3] for the middle reach of Mekong River.

Bayesian Error Classification and Analysis
Consider the task of obtaining forecasted discharges QF k (i+m) at some station j = 0 for a river.Predictand is the discharge Q(i+m) to be expected in the future, at m time intervals (days) from time instant i when the forecast is made.Assume that for this purpose the forecaster has available a number of forecast models, identified by index k (k = 0, 1, 2 . . .).When using these models for flood forecasting the two types of errors discussed in Section 1.2.2 and Section 1.2.3, i.e., system error eM k (i+m) and conditional error eF k (i+m) occur.The purpose of this section is to explore the nature of these errors by means of Bayes formula, as originally used by Krzysztofowicz [1], and further developed by Todini [2].Applications of the results are demonstrated in Section 3, by means of four models from Shahzad and Plate [3], which are briefly described in Section 3.1.

Preliminary Definitions
Models shall be identified by classes and types, which are investigated independently or in combination of any two models.Classes denote the number of models used for making a forecast.Model class 0 actually works without model, class 1 with 1 model, and class 2 with two models.Model type identifies models according to the model generation process.Type 0 is the persistence model with Q(i) as predictor for Q(i+m).Type 1 refers to a model based on time series of the historical record of discharges, and type 2 to rainfall-runoff models.For type 2, subtypes 2.1 and 2.3 are defined, as will be further explained in Section 2.4.3.
For simplifying multiple regression analysis, it is advantageous to convert all variables into dimensionless quantities with mean values of 0, and standard deviations of 1. Adapting the notation of Krzysztofowicz [1], dimensionless variables h, s k and h 0 are introduced: where, the overbars identify arithmetic means, and ϑ x is the standard deviation of h (with x = h), of s (with x = k) and of h 0 (with x = 0).Mean values and variances are determined during calibration, based on the assumption that records are long enough to yield stable estimates.Quantity h 0 is a special case of h, introduced separately to stay within the terminology of Krzysztofowicz, [43].It is the dimensionless forecast from the persistence model, with Q(i) as predictor for Q(i+m).Variables h, s k and h 0 have co-variances.COV !h, s k q " h ¨sk , COV !h, h 0 q " h ¨h0 , and COV !h 0 , s k q " h 0 ¨sk .In recognition of unavoidable output uncertainty of any model, forecasts should be described by a crisp value and an uncertainty band.The crisp value should be the best estimate, which is the expected dimensionless value h opt (i+m), i.e., the predictor obtained with well calibrated models.The uncertainty band is expressed through the pdf of the conditional forecast error, which for one model with forecast s k φ k th |s k u.For two models with outputs s k1 and s k2 it is φ k1,k2 th |s k1 , s k2 u.
It is interesting to note that in a forecast situation, the uncertainty expressed by the likelihood function must be of a regression type.This may be inferred from the following argument.Let h be the predictand and forecast s 1 obtained by model 1 for this quantity be the predictor.During calibration s 1 is found as the best estimate for actual value h, determined by means of the error pdf f{s 1 ´h} by least squares analysis.Thus, the forecast consists of mean value h and pdf f{s 1 ´h}.In a model calibration procedure this function has been used to minimize the variance of error s 1 ´h.In the forecasting case s 1 still is the best estimate for actual value h, so that when we look at the forecast h we know that h is (for given s 1 ) a random variable with expected value a 0,1 . s 1 , where a 0,1 is the regression coefficient.Consequently conditional pdf f{h|s 1 } changes into the error pdf f{h ´a0,1 . s 1 }.
Water 2015, 7, 6788-6809 Thus, the conditional forecast consists of mean value a 0,1 .s 1 and pdf f{h ´a0,1 . s 1 }.In the sense of Bayesian error analysis, this implies that we have prior information, which includes the mean value, which reduces the problem of Bayesian forecasting to linear regression analysis.

Class 0 Error Analysis
For class 0, one has no information that helps to improve a forecast, i.e., s 1 and s 2 are zero.The forecast is made by using only historical records of the predictand.Interpreted in terms of Bayesian estimation theory, the prior is non-informative, with posterior = likelihood function = Φ 00 thu .This is the situation of forecasters who have no information except past records of some time ago, or who must make a forecast for large values m and therefore can only use the pf or pdf of the total record.The best forecast that can be made on that kind of data base is mean value QF 00 (i+m) = Q pmq of the historical record, with error e 00 = Q(i+m) ´Q pmq and varianc ϑ 2 k .If this "model" (= model 00) is used as benchmark model, to assess forecast quality, it leads to the well-known Nash Sutcliffe criterion for model performance evaluation.

Class 1 Error Analysis
For class 1, one forecast model k exists with dimensionless output s k .Bayes formula for this case yields: This is the classical Bayes estimator used for predicting uncertainties due to changes [44].The left side of Equation ( 2) simply describes the conditional error pdf for a forecast with given model output, which is forecast s k .It is the conditional pdf for operational forecasts with model k, as shall be discussed in the remainder of this section.This is a strict result of applying the Bayes formula.In contrast, Equation (2) also may be used as basis for Bayesian estimation, in which case the right side is of interest for estimating effects of basic changes.It separates effects of changing h from effects of changes in model k.f{h} is the prior distribution obtained from historical data, which in an uncertain environment may reflect expectations of changes in the structure of the time series of predictands, i.e., due to climate change.Changes in likelihood function g{s k |h} reflect changes of the model.

Model Class 1 Type k = 0
In an actual forecasting situation, forecasters always have the value of the quantity to be forecast at time of forecast, i.e., they know value h 0 .In the context of forecasting h 0 is the output of model class 1 type 0 (called model 0), yielding s k = h 0 .We define this to be the case when predictor QF 0 (i+m) is discharge Q(i) of today, (persistence assumption) which in some cases, and in particular for short forecasting times or for large rivers, may be a good first guess.Forecasters can make a class 1, type 0 forecast by time shifting the observed record over historical times, and establishing a functional relationship h = q(h 0 ) between h and h 0 .In principle, q(h 0 ) could have any functional form depending on h 0 , but actually, since model outputs are optimum estimators, a linear relationship is plausible and will be postulated.For a linear correlation h|h 0 the left side of Equation ( 2) yields the conditional forecast probability density function (cpdf): Function ψ 0 { } is the error pdf.During calibration, this function can be obtained empirically, but in general errors for model class 1, type 0 would be purely random, so that ψ 0 { } is a normal pdf.Linear regression analysis of h with h 0 yields a 0 = h ¨h0 , so that the optimum dimensionless forecast is h opt = h ¨h0 ¨h0 , with conditional forecast error ε 0 " h ´h ¨h0 ¨h0 , which has variance σ 2 0 " σ 2 h r1 ´pa 0 q 2 s.Since σ 2 h = 1, this leads to dimensional forecast error e 0 = ε 0 ¨σ0 ¨ϑh with variance Water 2015, 7, 6788-6809 As will be shown for Mekong data in Section 3.1, the usefulness of this forecast in practical cases is rapidly reduced for longer forecast times.

Model Class 1 Type k =1
This type of forecast is generated by a time series model which produces runoff forecasts at station 0 using only past records of discharges-in the simplest case as a linear Markov regression model, or a higher level ARMA model based on discharges observed at different upstream gaging stations.A special feature of a model in class 1 type 1 (called model 1) is that no parameters have to be estimated during real time operation, all parameters are optimized during calibration and remain fixed during forecasting.For this case, Equation (3) applies, with h 0 replaced by s 1 .Since h and s 1 are random variables with mean 0 and variances 1, the conditional probability distribution of h is also given by Equation ( 3): where a 1 " s 1 h Error ε 1 = h ´a1 ¨s1 is a random variable with mean zero, and variance For application of forecasts with model 1, the dimensional predictand is QF 1 (i+m), which has conditional forecast error e 1 = ε 1 ¨σ1 ¨ϑh with variance σF 2 1 " ϑ 2 h σ 2 1 .

Model Class 1 Type k = 2
This class defines models based on inputs which have to be estimated in real time before a forecast can be made, for example models depending on inputs from meteorological variables.Typical for this class are rainfall-runoff models, for which rainfall inputs have to be forecast.In such a case two types of uncertainty arise.The first occurs during calibration.It is that of model class 1 type 2.1 (model 2.1), which has dimensionless forecast s 2.1 as output.For this, the future rainfall is known from historical records, but it has to be converted into the input to model 2.1 by means of a rainfall forecast model.s 2.1 is the best value that can be obtained from the model.It is not exact, because input calculated from this rainfall and input from true rainfall may be widely different.It includes errors both from the rainfall forecast model, and from the uncertainty due to estimating the actual rainfall distribution from the precipitation records.Errors due to observed rather than true rainfall become part of the model uncertainty and are reflected in error eF 2.1 (i+m).Consequently, forecast s 2.1 is a random variable with variance σ 2 2.1 , which is found empirically during calibration.Its conditional probability density cpdf φ 2.1 th |s 2.1 u has the functional form of Equation ( 4) for random variable h given s 2.1 , with indices appropriately changed.Function φ 2.1 th |s 2.1 u is determined during calibration by using the model to reconstruct historical data.It differs from the hydrologic uncertainty processor (HUP) of Krzysztofowicz [1], which yields outputs from the hydrological model under perfectly known rainfall inputs, whereas in our case rainfall inputs for Nash-cascades are estimates, although they are generated from actually observed (area averaged) future rainfall inputs from each of the sub-areas.
The second type of uncertainty occurs when the calibrated model is used during actual forecasting (case of model class 2 type 3 = model 2.3), resulting in forecasts s 2.3 .The conditional probability φ 2.3 th |s 2.3 u is that of operational forecasts s 2.3 , when a rainfall model has to be used and future rainfall inputs for the rainfall model have to be forecast.φ 2.3 th |s 2.3 u also is of the form of Equation ( 4), and determined by means of historical data in the calibration phase.It uses a forecast for the rainfall, and a rainfall model to convert observed rainfall into model input.We assume that random variable s 2.3 is bias free, and has variance σ 2 2.3 .

Relative Importance of Rainfall and Model Uncertainties
In the operational mode, model 2.  6) are assumed to be (approximately) Gaussian, and because a linear correlation h " a 2.3 ¨s2.3 exists, then one obtains, by integration of Equation ( 6): which also is a normal density distribution with variance a 2  9) for a rainfall runoff model in a given situation is small, then an effort to obtain better rainfall forecast data offers no advantage, and primary efforts should be directed to improve the models and/or their parameters.If the difference is large, improving forecasts of precipitation inputs or rainfall forecast should have priority.As an example, in Section 3.3, Equation (9) will be used to assess the relative importance of the rainfall forecast model on the total variance of the rainfall-runoff model by means of data from Shahzad and Plate [3].

Class 2 Error Analysis
For class 2 models, dimensionless predictand h is calculated by means of predictors s 1 and s 2 from two models, type k 1 and k 2 (called model k 1 ,k 2 ).Uncertainty bands of forecasts are derived from the joint probability density function (jpdf) f{h,s 1 ,s 2 } for the triplet h, s 1 , s 2 .Following where for any given pair s 1 and s 2 the integral on the right is a constant, which is known when making forecasts with both models k 1 and k 2 .With this expression one obtains from Equations ( 10) and ( 11): where φ th|s 1 , s 2 u is the conditional forecast cpdf to be determined.
Using the terminology of the Bayes estimation theory, function f p {h|s 2 } for any calculated pair s 1 s 2 , is the conditional prior, and g{s 1 |h,s 2 } the conditional likelihood function.However, in contrast to its use in Bayesian estimation, Equation ( 12) is the classical Bayes formula [44], connecting probabilities from multiple sources.

Model Class 2: Type k,0
Well studied is the case of linear combination of one (conceptual) forecast model with output s 1 and model 0 with output h 0. This is model class 2 type 0,k (called model 0,k).It was introduced by Krzysztofowicz [1] as part of his HUP (hydrological uncertainty processor).It is based on the fact that whenever one has a mathematical forecast model k yielding some output s k , one also has available the observed value of h 0 as additional predictor.Inclusion of initial values h 0 into forecasting models yields a three-dimensional pdf f{h, s k , h 0 } with variables h, s k and h 0 , so that Equation (12) becomes: This is the conditional probability density of h for given values of s k and h 0 .It is possible to determine φ k,0 th |s 1 , h 0 u either directly by means of multiple regression analysis, as used by Todini, [2] with his model conditional processor MCP, or by means of separate analysis of functions g ts k |h, h 0 u and f p th |h 0 u , as done for the HUP of Krzysztofowicz [43], and Krzysztofowicz and Herr [45].All functions have to be found from parallel calibration data for quantities h, s 1 , and h 0 .
If linear correlations exist among the three variables, with h dependent and h 0 and s 1 as independent variables, it is straightforward to use φ k,0 th |s 1 , h 0 u and to apply linear multiple regression analysis to h as function of s 1 and h 0 .Let the distribution of this linear combination be: Determination of coefficients of linear regression by means of least square fitting yields: Water 2015, 7, 6788-6809 Error ε k,0 " h ´ak,0 ¨sk ´bk,0 ¨h0 has variance: It should be noted that for this analysis of error ε k,0 no assumption was made on types of probability distributions of the three variables apart from assuming linearity of the relationship among them.Note: although during calibration the error variances can be found directly from the data, it is useful, as a control for accuracy, to also use Equation ( 16).

Model Class 2 Type w,0
A different analysis of class 2 type k,0 error is obtained if the problem is formulated as a Bayesian problem by means of the right side of Equation ( 13), as suggested by Krzysztofowicz and Herr [45] for their HUP.In this way they were able to use the properties of conjugate normal distributions [46].Krzysztofowicz named g{.}likelihood function, and f p {.} prior distribution, in analogy to Bayes estimation.Basic data are s k , h and h 0 , which for many rivers are found to be Weibull distributed [43].The data for the present study also were well represented by a Weibull distribution, so that the approach of [45] could be used.A method of fitting Weibull distributions to observed data is briefly described in an appendix.As an alternative, Krzysztofowicz and Herr [45] recommend to convert variables s, h, and h 0 into normal variables z, w, w 0 , respectively, through equality of normal probability function and Weibull probability function.

Prior Density Distribution
The Bayes formulation does not depend on specific forms of prior and likelihood functions, thus f p (w|w 0 ) and g tz|w, w 0 u may be any empirically found distribution to be fitted by a suitable distribution.Relation w|w 0 of prior f p {w|w 0 } is an unknown function, which in general also could have any functional form.The simplest type of dependency is a linear dependency w = a p w 0 with w 0 as independent variable, in which case function f p (w|w 0 ) also is a normal distribution.Krzysztofowicz [43] derived coefficient a p by connecting w and w 0 as a Marcov chain leading from w 0 to forecast w 1 for m = 1, where transition probability a p1 is based on past observation records.The same transition probability is then used to connect any two time steps up to time m, so that w|w 0 remains linear.As was shown by Todini [2] better results are obtained if instead of this Marcov chain model a direct linear regression of w (= w(m)) on w 0 is assumed, as will be used here.It leads to a normally distributed prior distribution: which is of the same form as the pdf of error class 1 Type 0, so that we obtain a p " ww 0 and σ 2 p " 1 ´pa p q 2 " 1 ´ww 0 2 , which yield the dimensionless error ε W 1,0 " w ´ww 0 ¨w0 .

Likelihood Density Distribution
Likelihood function g{.} is based on calibration data, for which w and w 0 are known, so that for purpose of forecasting z can be expressed as a function of w and w 0 .The simplest relationship again is a linear dependency of z on w and w 0 to yield a distribution: g tz|w, w o u " ψ g ˆz ´al ¨w ´bl ¨w0 σ l

˙(18)
Water 2015, 7, 6788-6809 which has the same form as Equation ( 14), so that we obtain best fit estimates: a l " zw 0 ´w0 w ¨w0 z 1 ´w0 w ¨ww 0 b l " w 0 z ´wz ¨w0 w 1 ´ww 0 ¨w0 w Posterior Distribution Posterior distribution φ th|s k , h 0 u is obtained as combination of the two distributions into one conditional forecast pdf: g tz|w, w 0 u ¨fp tw|w 0 u .Assuming this to be the product of two normal distributions, the solution is also a normal distribution, where exponents of both distributions are added and the result is expressed in terms of w as function of z and w o : g tz |w, w 0 u ¨fp tw |w 0 u " ψ w ˆw ´aw,0 ¨z ´bw,0 ¨w0 i.e., difference w ´aw,o ¨z ´bw,o w 0 is a random variable with variance σ 2 w,0 .As for Equation ( 14), linear least squares fitting yields coefficients a w,0 , b w0 and variance σ 2 w,0 : Finally the result is re-transformed into original variables.For this last step, original observed data Q(i+m) were graphically correlated with their normal counterparts and a nonlinear regression curve was fitted, which was used to re-transform normally distributed variables into the original Weibull distributed variables.
These relationships have been extensively used by Krzysztofowicz and his students.However, given the assumption of linearity and normal distributions, separation into prior and likelihood functions actually adds no additional information to the error analysis of class 2 type k,0 error.For if one inserts coefficients a l , b l and σ 2 l from Equations ( 18) and ( 19) into expressions a w,0 .b w,0 and σ 2 w,0 from Equation ( 21), equalities a w,0 = a 1,0 .b w,0 = b 1,0 and σ 2 w,0 = σ 2 1,0 are readily established [4].Thus, the Bayesian approach offers no advantage, if normality and linearity can be assumed.Differences may come only from conversion of original distributions of variables into normal distributions, because linear regressions in original space may not equal linear regressions in normal space.
In their application of these results, Krzysztofowicz and Herr [45] have added a model for forecasting rainfall based on conditions of no rain or with rain on day i.This modification of the basic theory was not used in the application of this model to the Mekong River, because analysis was restricted to rainy seasons with only very few dry days, and the simple persistence assumption for rainfall of today = rainfall of the future was made.

Model Class 2 Type 1,2
The influence of h 0 fades very quickly with increase of forecast time.But other models which have a better forecast quality can also be used in class 2 mode, as was suggested, for example, in [2] and used in Shahzad and Plate [3].A combination of output s 1 from a regression model of type 1 and output s 2 obtained from a rainfall-runoff model of type 2.3 is such a case of applying two forecast Water 2015, 7, 6788-6809 models for the same forecast situation.If a linear combination among h, s 1 and s 2 is assumed, then the left side of Equation ( 14) in the form φ 1,2 th |s 1 , s 2 u with s k = s 1 and h 0 replaced by s 2 leads to: The best forecast estimate then is: h opt " a 1,2 s 1 `b1,2 s 2 (23) with forecast error: ε 1,2 " h ´hopt " h ´a1,2 s 1 ´b1,2 s 2 and variance σ 2 1,2 .

Results of the Error Analysis
In the final form the forecast results are expressed through two quantities: a crisp forecast, i.e., expected value of the forecast and the standard deviation as a measure of the error band.In a dimensional form, the expected value of the predictand QF opt (i+m) is found by replacing s 1 and s 2 (or h 0 in combination with either k = 1 or k = 2) into Equations ( 22) and ( 23), and then converting h opt into dimensional form by means of Equation ( 1) to yield: with variance: h of the forecast error e k1,k2 " ε k1,k2 ¨ϑh where indices k identify model types.If the error pdf is Gaussian (as should be checked), then probabilities of any range can easily be obtained from mean and standard deviation.If the pdf is non-Gaussian, its empirical distribution may be found by ranking data from the lowest value with rank n = 1 to the highest value with rank n = n max , and assigning an empirical probability function n/(n max + 1) to the data of rank n.Alternatively, a transformation can be introduced that converts error data into Gaussian distributions, as was outlined in Section 2.4.2 for a Weibull distribution.For Gaussian distributions the distribution is fully specified by mean and variance of the data.From Gaussian exceedance probabilities these pdfs are re-transformed into original data space to obtain upper bounds, for example upper bounds for 80% of all errors.

Criterion for Forecast Quality
As measure of forecast error bands we use standard deviations σF k pmq of forecast pdfs.Additionally needed is a measure of forecast quality.A meaningful index is obtained from the persistence assumption, Q(i+m) = Q(i), which uses the fact that during a forecast, the present day value is always known and can be included as a special model case (as "benchmark" model, [36]).It leads to variance σF 2 0 , and to persistence index PI of Kitanides and Bras [39], defined as: For a perfect forecast error variance σF 2 k is zero, and PI is 1.PI is 0 if assuming persistence is just as good as, and is negative, if not better than the forecasting model.By using h 0 as reference, a forecaster has improved his forecasts as compared to that of the no data case (i.e., using the Nash-Sutcliffe index) by reducing the forecast error variance by a factor of: Water 2015, 7, 6788-6809

Application to Forecast Data for Mekong River
In order to illustrate the concepts of Section 2, data obtained from three models for forecasting Mekong River flows for the middle reach of the Mekong River (see Figure 1) were used.Two data based models, using very different model structures, are taken from Shahzad and Plate, [3].Here, we add a third model, the persistence model, which is called model 0. Models 1 and 2 were developed by Shahzad [47] using records of daily discharges Q j (i+m) from seven gaging stations j (shown in Figure 1) and rainfall sums for each day i from 37 rainfall stations, which were available for 15 years from 1991 to 2005.Approximate flow times between gages are one or two days.First a basic flow network was set up, which connects all gages j.An algorithm based on the continuity equation was developed [3], which allows the routing of area averaged contributions from the subareas between adjacent stations to reference station 0, where forecasts for up to m = 5 days were made.It is used for both models, which were calibrated using the same observed discharge differences between adjacent upstream gages.The difference between models 1 and 2 lies in their different forecast sub-models for the sub-areas between gages.Model 1 of Shahzad and Plate [3] is a typical time series model.Discharge increments for each sub-area are forecast by regression analysis with known increments of the past.Model 2 is a unit hydrograph rainfall-runoff model for the discharge increments due to runoff from the sub-areas between gages.A special feature of model 2 is the seasonally varying adjustment factor KN, which converts rainfall into effective precipitation.For details, paper [3] should be consulted.
Water 2015, 7 16 hydrograph rainfall-runoff model for the discharge increments due to runoff from the sub-areas between gages.A special feature of model 2 is the seasonally varying adjustment factor KN, which converts rainfall into effective precipitation.For details, paper [3] should be consulted.During calibration empirical parameters of both models were fitted to minimize variances of calibration errors e = Q(i+m) − QF(i+m), assuming Q(i+m) given (from historical data) and QF(i+m) the result of model calculation.In forecast applications, the situation is reversed, forecast QF(i+m) is given, and Q(i+m) is estimated by means of its conditional pdf.The purpose of the present study is to investigate this reverse process and to study effects of various error correction routines.The models originally were applied to all seven stations along the middle Mekong River [3].In the following, data During calibration empirical parameters of both models were fitted to minimize variances of calibration errors e = Q(i+m) ´QF(i+m), assuming Q(i+m) given (from historical data) and QF(i+m) the result of model calculation.In forecast applications, the situation is reversed, forecast QF(i+m) is Water 2015, 7, 6788-6809 given, and Q(i+m) is estimated by means of its conditional pdf.The purpose of the present study is to investigate this reverse process and to study effects of various error correction routines.The models originally were applied to all seven stations along the middle Mekong River [3].In the following, data and model outputs for station Stung Treng, located at the end of the middle reach of the Mekong shown in Figure 1, are used as a reference station, (identified as station j = 0) for illustrating the error analysis of Section 2. For this station, error statistics for all years, from 1991 to 2000, are derived for each forecast time of m days.For illustrations, discharge forecasts for three days ahead for year 1997 at station Stung Treng are used.

Application of Models 0, Model 1 and Model 0,1
Results from the analysis of errors for all models involving regressions (model 0, model 1 and model 0,1) are summarized in Table 1.Parameters for models listed in Table 1 were calculated from dimensionless discharge data for station Stung Treng for 1991 to 2000.Forecasts QF k (i+m) according to model k were obtained, which yield errors e k .Combinations of models are denoted by index k 1 ,k 2 .Typical results for forecasts by means of model 1 are shown in Figure 2 for m = 3. Top curves are observed and forecast seasonal discharges, and bottom curves conditional forecast errors: e 1 for single model dependency on model 1, e 0 shows single model dependency on the persistence model, and e 1,0 two model dependency on optimum combination of model 1 and persistence model.Included in Table 1, column 11 is the skew coefficient of the pdf of QF 0,1 (i+m), and the last column, 12.The skew coefficient is needed for Weibull fitting (see Appendix).Water 2015, 7 17

Application of Models 0, Model 1 and Model 0,1
Results from the analysis of errors for all models involving regressions (model 0, model 1 and model 0,1) are summarized in Table 1.Parameters for models listed in Table 1 were calculated from dimensionless discharge data for station Stung Treng for 1991 to 2000.Forecasts QFk(i+m) according to model k were obtained, which yield errors ek.Combinations of models are denoted by index k1,k2.Typical results for forecasts by means of model 1 are shown in Figure 2 for m = 3. Top curves are observed and forecast seasonal discharges, and bottom curves conditional forecast errors: e1 for single model dependency on model 1, e0 shows single model dependency on the persistence model, and e1,0 two model dependency on optimum combination of model 1 and persistence model.Included in Table 1, column 11 is the skew coefficient of the pdf of QF0,1(i+m), and the last column, 12.The skew coefficient is needed for Weibull fitting (see Appendix).

Table 1.
Results from combinations of model 1 and persistence model 0. Column 7 lists standard deviations of calibration error e = Q(i+m) − QF1(i+m), columns 2, 3 and 4 show correlation coefficients for all error equations, column 5 and 6 are parameters for Equation (15).All standard deviations (std) are in m³/s.The resulting forecast hydrograph QF1(i+m) has a shape very similar to that of model 0, as illustrated by means of the two error hydrographs in Figure 2. In both figures, substantial forecast errors are observed, in particular for longer forecast times, and for times with large discharge changes The resulting forecast hydrograph QF 1 (i+m) has a shape very similar to that of model 0, as illustrated by means of the two error hydrographs in Figure 2. In both figures, substantial forecast Water 2015, 7, 6788-6809 errors are observed, in particular for longer forecast times, and for times with large discharge changes over short times, where model outputs always lag behind observed hydrographs.Outputs of model 1 and of the combination of model 1 and model 0 (model 0,1, not shown) are practically identical, as seen from comparing columns 9 and 10.The reason for this behavior is that both model 1 and model 0 are regression models, and reduce errors in the same way.However, model 1 shaves off some part of the error peaks of the persistence model 0. This, perhaps, is not so evident in Figure 2, but it is documented in the reduction of the standard deviation listed in Table 1, where model averages over all 10 years of records are shown.The standard deviation of error e 1 relative to standard deviation of model 0 error e 0 is about 68% for m = 2, 76% for m = 3, and 82% for m = 5, leading to persistence indices PI of about 0.3 for m = 5 and of about 0.4 for m = 3.

Application of Model 2
Two different types of model 2 of Shahzad and Plate [3] are important: model 2.1, which is used to calibrate model 2 with historical input and output data, and model 2.3, which uses parameters of model 2.1, but with rainfall forecasts (applying the rainfall persistence assumption P(i+r) = P(i), (r = 1,2 . . .m) as rainfall forecast model).
Consider first, application of model 2.3 and the combination of model 0 and model 2.3, yielding forecasts QF 2.3,0 (i+m).Typical results are shown in Figure 3 for m = 3 for 1997 at gage Stung Treng.By comparing columns 7 and 9 in Table 2 we notice first that standard deviations of calibration errors e are almost the same as conditional forecast errors e2.3, as was also observed in other cases by Todini [48].In the forecast-hydrograph strong oscillations of QF 2.3 (i+m) are noted, which are caused by the strong variation of rainfall from day to day.Due to the persistence assumption for rainfall estimation the daily variation of rainfall is fully preserved.The oscillations could be reduced if the rainfall input was smoothed over some time period, i.e., by taking the average of the last three or four days.Such a smoothing filter leads to a reduction of the variance, but also to larger time lags for all stations and all times.
Water 2015, 7 18 over short times, where model outputs always lag behind observed hydrographs.Outputs of model 1 and of the combination of model 1 and model 0 (model 0,1, not shown) are practically identical, as seen from comparing columns 9 and 10.The reason for this behavior is that both model 1 and model 0 are regression models, and reduce errors in the same way.However, model 1 shaves off some part of the error peaks of the persistence model 0. This, perhaps, is not so evident in Figure 2, but it is documented in the reduction of the standard deviation listed in Table 1, where model averages over all 10 years of records are shown.The standard deviation of error e1 relative to standard deviation of model 0 error e0 is about 68% for m = 2, 76% for m = 3, and 82% for m = 5, leading to persistence indices PI of about 0.3 for m = 5 and of about 0.4 for m = 3.

Application of Model 2
Two different types of model 2 of Shahzad and Plate [3] are important: model 2.1, which is used to calibrate model 2 with historical input and output data, and model 2.3, which uses the parameters of model 2.1, but with rainfall forecasts (applying the rainfall persistence assumption P(i+r) = P(i), (r = 1,2…m) as rainfall forecast model).
Consider first, application of model 2.  2 we notice first that standard deviations of calibration errors e are almost the same as conditional forecast errors e2.3, as was also observed in other cases by Todini [47].In the forecast-hydrograph strong oscillations of QF2.3(i+m) are noted, which are caused by the strong variation of rainfall from day to day.Due to the persistence assumption for rainfall estimation the daily variation of rainfall is fully preserved.The oscillations could be reduced if the rainfall input was smoothed over some time period, i.e., by taking the average of the last three or four days.Such a smoothing filter leads to a reduction of the variance, but also to larger time lags for all stations and all times.It is evident that inclusion of model 0 in the combined forecasts of model 2.3 and 0 also has such a smoothing effect, as seen in Figure 3. Significantly smaller standard deviations of e0,2.3 as compared with e2.3 are documented in columns 9 and 10 of Table 2.This is reflected in column 12 of Table 2, which shows the improved persistence index for best results, corresponding to curve QF 0,2.3 (i+3) of   3. Errors due to rainfall forecast uncertainty are rapidly increasing with increased forecast time, reaching more than 50% for m = 5, indicating the need for improved rainfall forecasts for larger m, while the model error dominates short times m.

Application of Model 1 and Model 2 in Combination
The success of combining model 2 type 3 and model 0 suggests using model 1 instead of model 0 as part of a class 2 model, because model 1 and model 0 are both regression models, and model 1 actually is an improved version of model 0. Let this combination be called model 3, with index 1,2.3.Table 3 is a summary of model 3 parameters and standard deviations.It had already been suggested in [3] to use such a combination, and the generalization used here is a logical extension of class 2 error analysis, as it results in weighted averages of both models.The weights are listed in column 5 and 6 of Table 4, with a 1,2 the weight given to model 1, and b 1,2 the weight for model 2. Weights are seen to shift with increase of m from 0.66 for model 2 type 3 for m = 1 to 0.39 for m = 5.

Probability Densities of Final Errors
Probability densities were found by first normalizing all errors by division through their standard deviation.Their mean was zero in all cases.Then the data were ordered in classes of width 0.5 times standard deviation, relative numbers of values found per class were used as estimators for class probability, and empirical pdfs determined by drawing smooth curves through normalized relative magnitudes.
Results for all three optimum types considered in Sections 3.1 to 3.3 are shown in Figure 5.All data yield shapes distinctly different from normal distributions, also shown in the figures.They are much more concentrated around the center-showing a large proportion of small errors-and are somewhat skewed.Their skew coefficients are listed in columns 11 of Tables 1, 2, and 4. All hydrographs from model 2 show (see Figures 3 and 4) that the largest errors occur due to rapid changes in rainfall input, which is not forecasted by the rainfall model.

Probability Densities of Final Errors
Probability densities were found by first normalizing all errors by division through their standard deviation.Their mean was zero in all cases.Then the data were ordered in classes of width 0.5 times standard deviation, relative numbers of values found per class were used as estimators for class probability, and empirical pdfs determined by drawing smooth curves through normalized relative magnitudes.
Results for all three optimum types considered in Section 3.  3 and 4) that the largest errors occur due to rapid changes in rainfall input, which is not forecasted by the rainfall model.
Probability densities were found by first normalizing all errors by division through their standard deviation.Their mean was zero in all cases.Then the data were ordered in classes of width 0.5 times standard deviation, relative numbers of values found per class were used as estimators for class probability, and empirical pdfs determined by drawing smooth curves through normalized relative magnitudes.
Results for all three optimum types considered in Sections 3.1 to 3.3 are shown in Figure 5.All data yield shapes distinctly different from normal distributions, also shown in the figures.They are much more concentrated around the center-showing a large proportion of small errors-and are somewhat skewed.Their skew coefficients are listed in columns 11 of Tables 1, 2, and 4. All hydrographs from model 2 show (see Figures 3 and 4) that the largest errors occur due to rapid changes in rainfall input, which is not forecasted by the rainfall model.It is a special characteristic of the hydrological regime of the middle Mekong that these changes are difficult to anticipate.Rainfall maxima typically are generated by typhoons moving in from the South China Sea and occur almost simultaneously over large stretches along the Annamite (Truong Son) mountains, which form the boundary between Laos and Vietnam.The catchment characteristic of the middle Mekong are such that, instead of spreading the runoff from these rainfall extrema over time, they are superimposed and lead to large discharge peaks at downstream stations on the middle Mekong.The analysis of the final errors does not consider the heteroscedasticity of the error.Since critical forecasts are associated with strong jumps in the hydrograph, where unfortunately errors are largest, tail ends of distributions are of special significance.It is evident that probabilities for errors larger than two standard deviations are larger than normal.This must be taken into account if weighted decisions based on these tails are to be made.

Summary and Conclusions
In this paper, we looked at possible forecast improvements by means of combining two different models.We start with given forecast models, which are supposed to have been developed in a design phase.As an example, we use results from the Mekong River, i.e., the models of Shahzad and Plate [3], as optimum models.Parameters for the models had been optimized by means of historical data as part of design.The parameters are kept constant during forecast.
Due to scarcity of rainfall stations and/or model simplicity, models used for the Mekong River suffer from large inherent uncertainty-the error analysis showed that there is much room for model improvement, in particular if the rainfall forecast could be improved.Under any circumstances one must live with the system error eM k (i+m).It sets the limit of forecast quality.It can be reduced only if model or data base, or both, are improved.For a regression model without forecasted components, such as our model 1, error eM k (i+m) statistics are the same as forecast error eFk(i+m) statistics.For rainfall-runoff models, eM k (i+m) and eFk(i+m) are different: the former obtained with historical forecasts, the latter with forecasts based on rainfall forecast modeling, in our case assuming persistence of rainfall.The variance of eM k (i+m) must be used as reference against which to measure forecast error eF k (i+m), which is obtained in the operational or forecast mode as combination of system error and error due to imperfect forecast of input data.If the difference in the two errors for a given situation is small, then an effort in obtaining better forecast input data offers no big Water 2015, 7, 6788-6809 advantage, and the primary effort should be directed to improve the model and/or its parameters.If the difference is large, a major effort is needed for improving the forecast of the input data.The large differences of the outputs of models 2.1 (yielding eM k (i+m)) and model 2.3 (yielding eFk(i+m)) indicate that much could be gained for forecasts with the simple unit hydrograph model 2 if better rainfall forecasts were available.
The error analysis proceeds along concepts developed by Krzysztofowicz.Operational forecasts are identified as depending on the set of conditions existing at the time of forecasting, such as discharge or rainfall data at time i.This concept permits to combine contributions of different models in a logical and systematic fashion.We used three models: the simple persistence model 0 (assuming QF(i+m) =Q(i)), model 1, which is a regression model based on regressing discharge QF(i+m) on inputs from upstream reaches, and model 2, which is a rainfall-runoff model.Results of the error analysis showed that the best results were obtained by means of model 2.3 combined with model 1.The original data are well fitted by Weibull distributions, but no advantage was found from converting these data into normally distributed data, which after conversion are linearly correlated, as had been suggested by Krzysztorowicz [1].
An interesting observation is that the contribution of h 0 to the forecast pdf is much smaller for model 1 than for model 2, although from a direct comparison it would seem that the larger uncertainty of model 1 would make it more likely that a contribution of h 0 would matter.The reason for this behavior is that model 1 is a regression model, and the forecast for model 1 is more likely than model 2 to be comparable to the forecast from the persistence model, which also is a (rather elementary) regression model.Therefore one may conclude that the regression property of model 0 is already included in the general regression model, so that h 0 can make only a small contribution to the forecast by means of model 1.On the other hand, model 2 does not have a regression component.Therefore, model 0 complements the forecast from the rainfall-runoff model.This may imply that the use of a rainfall-runoff model profits more from the inclusion of h 0 -dependency than a regression model.
Water 2015, 7, 6788-6809 A semi-graphic method was developed in Plate [44] which is based on the fact that parameter α of Equation (A-1) depends on C sx only.Pairs of corresponding values of α and C sx are tabulated, and a universal diagram showing pµ th {βq α and σ{µ th as functions of 1/α is given in [44].With α known, the two parameters β and µ th " x ´x0 are determined from the universal diagram.With these parameter F x px, α, β, x 0 q is calculated, and by setting F x px, α, β, x 0 q " N `.x, µ x , σ x ˘normal variables .
x with mean µ x and variance σ 2 x of the untransformed variable x are obtained.

Figure 1 .
Figure 1.Middle reach of the Mekong, showing gaging stations (with permission of Mekong River Commission).

Figure 1 .
Figure 1.Middle reach of the Mekong, showing gaging stations (with permission of Mekong River Commission).

32 Figure 2 .
Figure 2. Forecast results for m = 3 days with model 0 and model 1 for Stung Treng in 1997.

Figure 2 .
Figure 2. Forecast results for m = 3 days with model 0 and model 1 for Stung Treng in 1997.
3 and the combination of model 0 and model 2.3, yielding forecasts QF2.3,0(i+m).Typical results are shown in Figure 3 for m = 3 for 1997 at gage Stung Treng.By comparing columns 7 and 9 in Table

Figure 3 .
Figure 3. Forecast results for 3 days with model 2.3 and linear combination of model 2.3 with model 0 for Stung Treng in 1997.Additionally shown is the forecast of model 2.1.

Figure 3 .
Figure 3. Forecast results for 3 days with model 2.3 and linear combination of model 2.3 with model 0 for Stung Treng in 1997.Additionally shown is the forecast of model 2.1.

Figure 3 .
As comparison of column 9 and 10 shows, the combination of model 0 with model 2.3, in the multiple regression mode, improves the standard deviation of e0,2.3 by 21 % for m = 1 and 11% Water 2015, 7, 6788-6809 for m = 5, yielding improvement of persistence index PI from 0.40 for m = 5 to 0.58 for m = 1 (see column 12).

Figure 4 .
Figure 4. Forecast results for m = 3 days showing model 1, model 2 type 3, and linear combination of the two models to yield model 1,2.3 for Stung Treng in 1997.

Figure 4 .
Figure 4. Forecast results for m = 3 days showing model 1, model 2 type 3, and linear combination of the two models to yield model 1,2.3 for Stung Treng in 1997.
1 to Section 3.3 are shown in Figure 5.All data yield shapes distinctly different from normal distributions, also shown in the figures.They Water 2015, 7, 6788-6809 are much more concentrated around the center-showing a large proportion of small errors-and are somewhat skewed.Their skew coefficients are listed in columns 11 of Tables 1, 2 and 4. All hydrographs from model 2 show (see Figures

Figure 5 .
Figure 5. Probability densities of errors for cases of model combinations.(a) pdf of error e 1,0 ; (b) pdf of error e 2.3,0 ; (c) pdf of error e 1,2.3 .
1 cannot be used, because neither the true rainfall model nor the actual rainfall input is known.However, from a comparison of φ 2.3 th |s 2.3 u and φ 2.1 th |s 2.1 u it is Water 2015, 7, 6788-6809 possible to infer that part of the variance which is due to uncertainty of rainfall information.For this, we recognize that forecasts s 2.3 are related to forecast s 2.1 through a conditional pdf f ts 2.1 |s 2.3 u .In our definition, f ts 2.1 |s 2.3 u is the precipitation uncertainty processor, which is unknown, because rainfall inputs are estimated quantities, which may have large errors.Variance σ 2 02 of f ts 2.1 |s 2.3 u can be determined by means of the following argument.Calibration using the assumed rainfall model and forecast rainfield at time m ∆t yields φ s th |s 2.3 , s 2.1 u .Applying the Bayes formula, this can be written as: φ s th |s 2.3 , s 2.1 u " φ 2.1 th |s 2.1 u ¨f ts 2.1 |s 2.3 u (5) where f ts 2.1 |s 2.3 u is the unknown cpdf for the output s 2.1 due to the unknown true rainfall, conditioned on the known calculated output s 2.3 .Conditional pdf φ 2.1 th |s 2.1 u is obtained from model 2.1 with known rainfall fields.It has mean a 2.1 ¨s2.3 and variance σ 2 2.1 obtained during calibration with historical data.From Equation (5), φ 2.3 ph ˇˇs 2.3 q is obtained by marginalization over s 2.1 to yield: φ 2.1 th |s 2.3 u " th |s 2.1 u ¨f ts 2.1 |s 2.3 u ¨ds 2.1 (6)If both φ 2.1 th |s 2.1 u and f ts 2.1 |s 2.3 u in Equation ( th |s 2.3 u can be obtained directly from the data by means of model 2.3, with the result: Note that during model calibration a comparison of variances for e 2.1 (i+m) and e 2.3 (i+m) indicates where efforts of improvements of model or model input might be most effective.If difference Equation ( Water 2015, 7, 6788-6809Krzysztofowicz [1], the desired error pdf f th |s 1 , s 2 u is obtained by decomposing f{h,s 1 ,s 2 } by means of the Bayes formula: f th, s 1 , s 2 u " φ th |s 1 , s 2 u ¨f ts 1 , s 2 u " g ts 1 |h, s 2 u ¨fp th|s 2 u ¨f2 ts 2 u ´8 g ts 1 |h, s 2 u ¨fp th |s 2 u ¨f2 ts 2 u ¨dh " f 2 ts 2 u ¨ż 8 ´8 g ts 1 |h, s 2 u ¨fp th |s 2 u ¨dh

Table 1 .
(15)lts from combinations of model 1 and persistence model 0. Column 7 lists standard deviations of calibration error e = Q(i+m) ´QF 1 (i+m), columns 2, 3 and 4 show correlation coefficients for all error equations, column 5 and 6 are parameters for Equation(15).All standard deviations (std) are in m³/s.

Table 2 .
(15)lts from combinations of model 2.3 and persistence model.Columns 2-4 show parameters (correlation coefficients) for all error equations, column 5 and 6 are parameters for Equation(15).Column 7 lists standard deviations (std) of calibration error e = Q(i+m) ´QF 1 (i+m).All std are given in m³/s.An advantage of Bayesian analysis for model class 1 type 2 is the separation of model uncertainty of the hydrological model from the uncertainty of the rainfall forecast model, as outlined in Section 2.3.4.Evidently, the difference in discharge forecasts from the two versions of model 2 is attributable to rainfall forecast inaccuracy.The best forecast with model class 2 is obtained with historical rainfall (model 2.1) in combination with forecast from model 1.An analysis of this case shows that the rainfall-runoff model 2.1 with inclusion of model 1 (column 4 of Table3) is practically the same, with and without model 1 (column 3 of Table3).Therefore, we use the standard deviations of column 4 as the best model with rainfall known.According to Section 2.3.4,thedimensionalvariance for model 2.1 is σF22.1 " σ 2

Table 3 .
Effect of rainfall forecast on forecast accuracy for combinations of model class 2 type 2.1 and 2.3. m

Table 4 .
(8)ults from combinations of model 2 and model l.Columns 2, 3, and 4 show parameters (correlation coefficients) for all error equations, column 5 and 6 are parameters for Equations(7)and(8).All std in m³/s.Table4are results from using linear multiple regression of variables w obtained from transforming originally Weibull distributed data into normal.The method described in Section 2.4.3 is applied to the combination of model 1 and model 2.3.A very small improvement of about 1% in variance reduction was obtained, so small that it does not seem to justify the effort of conversion into normal variables, in particular since no theoretical reason can be given why linear regressions among transformed variables is better than linear regressions among the original variables.Typical results illustrating model 3 outputs are shown in Figure 4. Bottom curves are conditional forecast errors: e 1 for single model dependency on model 1, e2.3 shows single model dependency on model 2.3, and e1,2.3 two model dependency on the optimum combination 1,2.3 of model 1 and model 2.3.The results of column 12, with PI values roughly decreasing from 0.6 for small m to 0.47 for m = 5 are the best results obtainable with the present set of models.