Improving Operational Short- to Medium-Range (SR2MR) Streamflow Forecasts in the Upper Zambezi Basin and Its Sub-Basins Using Variational Ensemble Forecasting

Abstract

The combination of Hydrological Models and high-resolution Satellite Precipitation Products (SPPs) or regional Climatological Models (RCMs), has provided the means to establish baselines for the quantification, propagation, and reduction in hydrological uncertainty when generating streamflow forecasts. This study aimed to improve operational real-time streamflow forecasts for the Upper Zambezi River Basin (UZRB), in Africa, utilizing the novel Variational Ensemble Forecasting (VEF) approach. In this regard, we describe and discuss the main steps required to implement, calibrate, and validate an operational hydrologic forecasting system (HFS) using VEF and Hydrologic Processing Strategies (HPS). The operational HFS was constructed to monitor daily streamflow and forecast them up to eight days in the future. The forecasting process called short- to medium-range (SR2MR) streamflow forecasting was implemented using real-time rainfall data from three Satellite Precipitation Products or SPPs (The real-time TRMM Multisatellite Precipitation Analysis TMPA-RT, the NOAA CPC Morphing Technique CMORPH, and the Precipitation Estimation from Remotely Sensed data using Artificial Neural Networks, PERSIANN) and rainfall forecasts from the Global Forecasting System (GFS). The hydrologic preprocessing (HPR) strategy considered using all raw and bias corrected rainfall estimates to calibrate three distributed hydrological models (HYMOD_DS, HBV_DS, and VIC 4.2.b). The hydrologic processing (HP) strategy considered using all optimal parameter sets estimated during the calibration process to increase the number of ensembles available for operational forecasting. Finally, inference-based approaches were evaluated during the application of a hydrological postprocessing (HPP) strategy. The final evaluation and reduction in uncertainty from multiple sources, i.e., multiple precipitation products, hydrologic models, and optimal parameter sets, was significantly achieved through a fully operational implementation of VEF combined with several HPS. Finally, the main challenges and opportunities associated with operational SR2MR streamflow forecasting using VEF are evaluated and discussed.

1. Introduction

1.1. Decisions and Limitations of Hydrologic Forecasting

At any spatiotemporal scale, critical decisions about the design, functionality, and operability of a Hydrologic Forecasting System (HFS) need to be made to reduce the total hydrologic uncertainty (THU) propagated from different components of a hydrologic modelling paradigm (HMP). In fact, reducing the total hydrological uncertainty is key to developing reliable Integrated Water Resources Management (IWRM) strategies across spatial and temporal scales. For river basins across the world, the allocation of water resources largely relies on accurate streamflow forecasts. In Africa, for instance, the waters of the Upper Zambezi River Basin (UZRB) are shared by eight countries: Angola, Namibia, Zambia, Botswana, Malawi, Tanzania, Zimbabwe, and Mozambique. However, the administrative complexities created by the transnational nature of the Zambezi Basin (Figure 1) result in inconsistencies in the operation and maintenance of the rain gauges and stream gauges (see Table 1), and consequently a lack of reliable hydrologic data for the implementation of an HFS. For example, rainfall or streamflow time series with missing records can undermine the effectiveness of calibration and validation schemes, consequently increasing the propagation of total hydrological (meteorological) uncertainty for the establishment of hydrologic processing strategies. Therefore, an appropriate identification and quantification of uncertainty (at any level) can help reducing the THU for the final development of streamflow forecast products.

1.2. Knowledge Gaps and Justification of the Study

The first HFS for Zambezi was applied in 2011, where the use of the Kalman filter was combined with a simple two-layer conceptual hydrological model to forecast streamflow in three sub-basins [3]. Satellite-based soil moisture estimates were used to calibrate the aggregated hydrologic model, which was able to generate daily streamflow forecasts up to 40 days into the future. Meier [4] argued that the spatial resolution of the satellite data needed to be improved together with the implementation of more sophisticated hydrologic models. In 2014, the daily floodplain behavior of the Zambezi was simulated by [5] applying a modified reservoir approach for the SWAT model [6]. Their results showed that the modified version of SWAT improved the simulation of daily streamflow and floodplain development in the Zambezi basin. Several other hydrological models have been satisfactorily calibrated and validated in other sub-basins of the Zambezi with available records (i.e., [7,8,9,10,11,12]). Despite all these modelling efforts in the Zambezi Basin, an operational HFS for the undisturbed flows of the poorly gauged UZRB has yet to be established. The primary objective of an operational HFS is to generate an accurate short- to medium-range (SR2MR) streamflow forecast that can inform water distribution schemes at the relevant spatial and temporal scales. The forecasts can be obtained from multiple ensembles constructed from Variational Ensembles Forecasting (VEF) approaches.

1.3. Variational Ensemble Forecasting (VEF) to Improve Operational Streamflow Forecasts

An operational HFS can only be implemented if the realtime and short-term forecasts of the input data (e.g., rainfall, temperature, etc.) are readily available. Therefore, the use of multiple satellite precipitation products (or regional climate models), combined with multiple conceptual and physic-based hydrologic models, can provide important insights about the practical and scientific aspects of implementing operational HFS in poorly gauged basins [13,14]). In this regard, many Variational Ensemble Forecasting (VEF) algorithms have been proposed to improve the representation of the components included in a Hydrologic Forecasting System (HFS). Previous studies have applied VEF based on multiproduct, multimodel, or multi-initialization schemes (see for example [13,14,15,16,17,18,19,20,21]). However, the evaluation of VEF including an additional dimension with multiple optimal parameter sets has not been explored in simulation schemes or operational forecasting yet’ neither has the role of Hydrologic Processing Strategies to improve the assimilation of VEF applications in an operational HFS context.

1.4. Purpose of This Paper

The main motivation of this study is to describe, analyze, and discuss the main steps required to design and implement an HFS, aimed to improve SR2MR streamflow forecasts in the UZRB and its sub-basins during the operational stage. The technical functionality and operability of the HFS is assessed using Hydrologic Processing Strategies (HPS) within the context of a Variational Ensemble Forecasting (VEF) hydrologic modelling paradigm [22], i.e., optimal combination of multiple precipitation products, multiple hydrologic models, and multiple parameters sets. The THU propagated by the HFS is evaluated and compared according to different sources of uncertainty, i.e., satellite-based, or model-based rainfall estimates, hydrological models, and optimal parameters. As detailed later in this paper, a more comprehensive overview of the whole modeling paradigm used to implement an operational HFS can help significantly to generate an improved streamflow forecast products that are closer to the hydrological conditions of the basin under study.

2. Methods

2.1. The Upper Zambezi River Basin (UZRB) Domain

The operational HFS implementation for the UZRB and its sub-basins with records (Table 1), considered the drainage area (~339,521 km²) delineated above the Katima Mulilo streamgauge (Table 1 and Figure 1a). The mean annual streamflow at the UZRB (measured at Victoria Falls) represents about 25% of the mean annual streamflow measured at the Zambezi Delta outlet (~4200 m³/s), the largest contribution of all tributaries within the whole Zambezi Basin [23,24]. The historic mean daily flows (1942–2017) measured at Katima Mulilo stream gauge were about 1174 m³/s (Table 1); however, during extreme episodes the maximum daily streamflows can exceed more than six times the mean daily streamflows. The basin is dominated by broadleaved trees (~66%), herbaceous (16.1%), and shrubs (14.8%), whereas only a little (~0.6%) of the area is managed or represents agriculture (Figure 1c). The spatial distribution of these vegetation types is consistent with the elevational pattern of the basin, which ranges between approximately 731 and 1671 m above sea level (Figure 1d). The UZRB domain for the operational HFS implementation was set up using grid cells at 0.25° of spatial resolution, approximately 25 km at the Equator (Figure 1e). This area above Katima Mulilo was selected, because it is the unique portion of the Zambezi that does not have ongoing hydropower or water storage infrastructure projects (see Figure 1b). This is an advantage for HFS implementation since the presence of dams or any other anthropogenic water regulations can be an important limitation when the primary objective is to simulate natural streams.

Table 1. Description of streamgauges use in this study. The HFS column shows those streamgauges (sub-basins) used for the calibration of hydrologic models.

Country	Streamgauge	South Latitude	East Longitude	Area (km²)	Altitude (m.a.s.l.)	Average Flow (m3/s)	Period	Missing (%)	HFS
Zambia	Kalene Hill Road Bridge	−11.13	24.25	764	1261	12.3	1977–2004	34.81	No
Zambia	Chivata Village	−13.33	23.15	3354	1065	17.4	1962–2004	23.32	Yes
Zambia	Luanginga-Kalabo	−14.96	22.68	34,621	1021	59.0	1958–2004	8.94	Yes
Zambia	Kabompo Pontoon	−13.60	24.21	42,740	1029	252.2	1990–2005	51.04	Yes
Zambia	Lukulu	−14.38	23.233	206,531	1012	772.0	1950–2004	12.44	Yes
Zambia	Senanga	−16.11	23.25	284,538	992	972.6	1947–2004	8.40	Yes
Namibia	Katima Mulilo	−17.48	24.3	339,521	746	1174.5	1942–2017	13.54	Yes

Table 1. Description of streamgauges use in this study. The HFS column shows those streamgauges (sub-basins) used for the calibration of hydrologic models.

Country	Streamgauge	South Latitude	East Longitude	Area (km²)	Altitude (m.a.s.l.)	Average Flow (m3/s)	Period	Missing (%)	HFS
Zambia	Kalene Hill Road Bridge	−11.13	24.25	764	1261	12.3	1977–2004	34.81	No
Zambia	Chivata Village	−13.33	23.15	3354	1065	17.4	1962–2004	23.32	Yes
Zambia	Luanginga-Kalabo	−14.96	22.68	34,621	1021	59.0	1958–2004	8.94	Yes
Zambia	Kabompo Pontoon	−13.60	24.21	42,740	1029	252.2	1990–2005	51.04	Yes
Zambia	Lukulu	−14.38	23.233	206,531	1012	772.0	1950–2004	12.44	Yes
Zambia	Senanga	−16.11	23.25	284,538	992	972.6	1947–2004	8.40	Yes
Namibia	Katima Mulilo	−17.48	24.3	339,521	746	1174.5	1942–2017	13.54	Yes

2.2. Forecasting Timescales and Water Management Activities

An HFS can be implemented for either of the three (or a combination thereof) principal forecasting timescales: (1) Realtime Monitoring or Short- to Medium-Range Forecasting (SR2MR), (2) Subseasonal to Seasonal or Short- to Long-Range Forecasting (SR2LR), or (3) Climate Change Predictions. The choice of forecasting timescales is partly determined by the relevant water management goals, and the needs of end users of the forecasts. In this study, the need for better water management schemes for flood warning and water allocation in the lower Zambezi Basin prompted the implementation of an HFS for realtime streamflow monitoring and short- to medium-range (SR2MR) forecasting in the UZRB and its sub-basins. As mentioned in the previous section, the HFS was implemented to forecast streamflow in the UZRB because it is the only part of the whole Zambezi River Basin without water regulation infrastructure, i.e., natural streamflow patterns can be observed in the UZRB. The SR2MR scheme allows streamflow forecasting from realtime up to eight days into the future (Figure 2).

2.3. The Operational Context of a Hydrologic Modeling Paradigm

A functional HFS requires the design of a hydrologic modelling paradigm under an operational scheme. The design can range from a simple propagation of meteorological forecasts using a single hydrological model to more advanced techniques which can include multiparameter, multimodel, or multi-initialization schemes (Figure 3). The use of a simple propagation scheme does not allow for the quantification of the hydrological uncertainty propagated from the model structure or from the model parameters. However, by adding more parameter sets for a single hydrological model (i.e., [15]), by adding more hydrological models [13,14,16,17,18], or simply by performing more initializations of the initial starting conditions of a single model [19,20], the hydrological uncertainty either from the parameters, model assumptions, or initial conditions can be quantified.

In the UZRB, the operational HFS is designed as a VEF modelling paradigm [22] in which the best SR2MR forecasts are derived from the combination of multiple precipitation products, multiple hydrologic models, and multiple parameters sets (Figure 3). This technique makes it possible to increase the range of possible streamflow forecast ensembles that can be used to evaluate and select the best representations of the hydrological states and fluxes for the UZRB and its sub-basins.

2.4. Selection of Hydro Climatological Forcings

The input data can have significant impacts on the propagation of meteorological and total hydrological uncertainty, and consequently, on the final streamflow forecasts. In this regard, it is well known that rainfall is the most important variable for streamflow simulation. At this stage, either point-based instrumental records or gridded-based climatology products can be used to establish baseline conditions and to correct SPPs or RCMs. The UZRB is a poorly gauged basin that lacks a consistent hydrometric network with continuous records (only nine rain gauges with discontinuous records were available for the HFS implementation). Therefore, the hydrological preprocessing hypothesis (also known as meteorological postprocessing hypothesis) [22] was approximated using rainfall climatology provided by CHIRPS (Climate Hazards Group InfraRed Precipitation with Station data) [25,26], and temperature climatology provided by the Global Meteorological Forcing Dataset (GMFD) [27].

To minimize the probability of errors propagated during the implementation of the operational HFS, the quality of the climatological forecast products used in an HFS needs to be evaluated and validated before they can be reliably used for hydrological applications. In the UZRB, three SPPs were evaluated, corrected, and then used to provide short-range (realtime) rainfall estimates (1) TMPA-RT [28]; (2) CMORPH [29]; and (3) PERSIANN [30,31]. To complete the SR2MR scheme, medium-range rainfall forecasts from the Global Forecasting System (GFS) [32] were used. The original GFS product provides rainfall forecasts in six-hour intervals (00, 06, 12, and 18 UTC) up to 16 days into the future. However, since the quantity of missing records between 9 and 16 days in the archives is larger (at least in this region), only rainfall forecasts provided at 00 UTC and up to eight days into the future were selected for operational implementation. (

Table 2. Climatology data used to correct (preprocess) rainfall forecasts from SPPs and RCMs.

Product	Institution	Spatial Resolution	Temporal Resolution	Global Coverage	Period
CHIRPS 1 [26]	UCSB	0.25° × 0.25°	Daily	50° N–50° S 180° E–180° W	1981 to present
GMFD 2 [27]	Princeton	0.25° × 0.25°	Daily	50° N–50° S 180° E–180° W	1981–2012
PERSIANN-CCS 3 [31]	UCI	0.25° × 0.25°	3-hourly	37.8° N–40.6° S 28.0° W–56.2° E	1998 to present
CMORPH 4 [29]	NOAA-CPC	0.25° × 0.25°	3-hourly	60° N–60° S 180° E–180° W	1998 to present
TMPA-RT 5 [28]	NASA GES DISC	0.25° × 0.25°	3-hourly	50° N–50° S 180° E–180° W	1998 to present
GFS 6 [32]	NOAA-NCEI	0.25° × 0.25°	3-hourly	90° N–90° S 180° E–180° W	2014 to present
Katima Mulilo	GRDC	Streamgauge	Daily	17.48° S–24.3° W	1942 to present

¹ Climate Hazards Infrared Precipitation with Station data (CHIRPS). ² Global Meteorological Forcing Dataset (GMFD) ³ Precipitation Estimation from Remotely Sensed Information using Artificial Neural Networks-Cloud Classification System (PERSIANN) ⁴ Climate Prediction Center morphing method (CMORPH) ⁵ Tropical Rainfall Measuring Mission (TRMM) Multi-satellite Precipitation Analysis-Real Time ⁶ Global Forecast System (GFS).

).

2.5. Hydrologic Models for Operational HFS

The selection of a model structure to be included in the hydrologic modelling paradigm (HMP) and ultimately in the operational HFS is as important as the selection of the model parameters to be used in the calibration process [16]. For an operational HFS, many HMP options can be implemented, including a single selection or an optimal combination of multiple precipitation or climate products, multiple hydrological models, multiple model state initializations, and/or multiple parameters sets. For example, the operational VEF approach implemented in the UZRB, and its sub-basins utilized three distributed hydrologic models (at 0.25° of spatial resolution): (1) HBV_DS, (2) HYMOD_DS, and (3) VIC 4.2.b (Figure 4). The first two models are distributed versions of two traditional well-known hydrological models: HBV [33,34] and HYMOD [35,36]. The VIC model used in this study is a modified version of the well-known Variable Infiltration Capacity (VIC) land surface model [37] that can resolve both the water and the energy balances. The modification allows postprocessing of VIC model outputs with the Lohmann’s model for routing [38]; a Gamma distribution to represent the catchment’s unit hydrograph; and the linearized version of the Saint-Venant Equations for final river routing. Additional details of the model states, fluxes, and parameters used in this study are provided in the Appendix A.

2.6. Calibration of Models Included in the HFS

The models described in the previous section were calibrated and validated for the whole UZRB and its sub-basins at 0.25 degrees of spatial resolution (approximately 25 km at the Equator). Historical available daily streamflow records for the UZRB and its sub-basins were used as reference data (Figure 5 and Figure 6). The genetic algorithm [39] was used to optimize the parameter sets of the three models. However, at this stage, any suitable optimization scheme can be implemented (e.g., [39,40,41,42,43]; and many others) based on the availability of time, resources, and expertise. In this study, a daily pooled calibration considering all observed daily records in the UZRB, and its sub-catchments (Figure 5) was applied at the Massachusetts Green High-Performance Computing Center (MGHPCC). The algorithm was run in parallel processing mode using 100,000 iterations, with population sizes ranging between 100 and 1000 generations. With this approach, the evaluation of the most appropriate population of parameter sets within each generation could be conducted in a more efficient manner. The Kling-Gupta efficiency [44] and the Nash-Sutcliffe efficiency [45], among other measures, were used to evaluate the degree of agreement between observed and simulated streamflow (Figure 6).

2.7. Operational Variational Ensemble Forecasting (VEF)

A group of i hydroclimatological inputs, j hydrologic models, and k calibrated model parameters can be used to establish simple or variational ensemble forecasts. A VEF approach can provide a larger number of forecasts than a simple assembling approach, because it evaluates all possible modelling chain sequences that can be arranged to construct an HFS (see Figure 7). At this stage, data assimilation techniques to support the “perfect model assumption” can be also applied to generate deterministic streamflow forecasts. Ref. [46] reviewed and evaluated many available methods to perform deterministic forecasts. However, in an operational context, the experience of streamflow forecasters suggests that uncertainty bands can better support decision making for water management schemes [47,48,49,50,51,52]. The intention of a VEF approach is to use all available weighted or non-weighted components of a modelling chain (HFS) to generate streamflow forecasts. The final hypothesis about the probability distribution of streamflow forecasts can be approximated using several procedures applied to the Multiproduct, Multimodel, and Multiparameter sets, through the implementation of a VEF approach that can be trained during T_W and used during T_F as:

$$ℳ(z_{y_{T_F}}^{i,j,k}|{\hslash }_{u_{T_F}}^i,{\hslash }_{p_{T_W}}^j,{\hslash }_{{\theta }_{T_W}}^k)=\mathscr{H}\left(\hslash (z_{y_{T_F}}^{i,j,k}|{\hslash }_{u_{T_F}}^i,{\hslash }_{p_{T_W}}^j,{\hslash }_{{\theta }_{T_W}}^k)\right)\forall T_W\in \left\{t_0,\dots ,t\right\} \mathrm{and} \forall T_F\in \left\{t,\dots ,t_f\right\}$$

(1)

$${\hslash }_{u_{T_F}}^i=\hslash (u_{T_F}^i|z_{u_{T_W}}^i,z_{u_{T_F}}^i)$$

(2)

$${\hslash }_{p_{T_W}}^j=\hslash ({\mathrm{S}}_{T_W}^j|{\hslash }_{{\theta }_{T_W}}^k)$$

(3)

$${\hslash }_{{\theta }_{T_W}}^k=\hslash ({\theta }_{T_W}^k|z_{{\theta }_{T_W}}^k)$$

(4)

$${\hslash }_{y_{T_F}}^{i,j,k}=\hslash (z_{y_{T_F}}^{i,j,k}|{\hslash }_{u_{T_F}}^i,{\hslash }_{p_{T_W}}^j,{\hslash }_{{\theta }_{T_W}}^k)$$

(5)

• where, $ℳ$ is a Multiproduct, Multimodel, and Multiparameter Variational Ensemble Streamflow Forecast for the forecast period T_F.
• $\mathscr{H}$ is a family of hypothetical ensemble components for the warmup period T_W and used to forecast the period T_F.
• ${\hslash }_{y_{T_F}}^{i,j,k}$ is a hypothesis of the hydrologic process from a family of input data i, hydrologic model j, and parameter set k, about the HFS for the forecast period T_F.
• $z_{y_{T_F}}^{i,j,k}$ is the streamflow prediction i, j, k about the HFS response for the forecast period T_F.
• ${\hslash }_{u_{T_F}}^i$ is a hypothesis of the input data from a family of input data i about the HFS for the forecast period T_F.
• ${\hslash }_{p_{T_W}}^j$ is a hypothesis of the hydrologic process from a family of hydrologic models j and for the warmup period T_W.
• ${\hslash }_{{\theta }_{T_W}}^k$ is a hypothesis of the parameter sets from a family of parameter sets k for the warmup period T_W.
• $u_{T_F}^i$ is the control variable for the forecast period T_F.
• $z_{u_{T_F}}^i$ is a family of input data i about the HFS for the forecast period T_F.
• $z_{u_{T_W}}^i$ is a family of input data i about the HFS for the warmup period T_W.
• ${\mathrm{S}}_{T_W}^j$ is the model structure j for the warmup period T_W.
• ${\theta }_{T_W}^k$ is the parameter set k for the warmup period T_W.
• $z_{{\theta }_{T_W}}^k$ is a family of parameter sets k for the warmup period T_W.

The VEF approach (Figure 7) makes it possible to explicitly represent how different components of each model ensemble (modelling chain) vary according to all possible combinations of biogeophysical representations of the climate system and the hydrologic system (e.g., [53]). In addition, one of the novel aspects of this research is introduction of the assumption that the forecast skill of hydrologic models to represent underlying hydrologic processes can be better captured by assessing the generalization capabilities of multiple optimal parameters sets.

2.8. Strategies to Reduce Uncertainty and Improve VEF in an Operational Environment

In the development of an operational streamflow forecasting paradigm, it must be determined how total hydrological uncertainty will be reduced to improve streamflow forecasts. This topic is still a matter of discussion among many hydrologists and researchers around the world [47,48,49,50,51,52]. For instance, standardized processes to decompose, quantify, or evaluate the meteorological or the total hydrological uncertainty propagated from a modelling paradigm are required (see, i.e., Figure 8). The implementation of new techniques for the decomposition of uncertainty can take advantages from VEF approaches to explore all available sources of climate data and physical representations (i.e., model structures and parameters) for hydrologic modelling. Taking this into account, a combination of VEF and hydrological processing hypotheses can help us better understand how the propagation of errors occurs. For example, from multiple climate products the amount of meteorological uncertainty propagated through the modelling chain can be identified; then comparisons can be made to estimate the amount of total hydrological uncertainty (THU) propagated from the same system. This more systematic method to identify uncertainties can be useful to inform additional pre- or postprocessing of THU.

Assuming that an operational HFS can be evaluated as a VEF approach, three main strategies (Figure 8) to establish hydrological processing hypotheses can be applied to evaluate and improve the forecast skill of any hydrological model: (1) Hydrological Pre-processing (HPR) (or Meteorological Postprocessing), which can include the application of bias correction techniques (i.e., physical or statistical improvements) to reduce the propagation of errors or the application of Bayesian or non-Bayesian approaches to estimate an optimal weighted combination of precipitation products; (2) Hydrological Processing (HP), which can include improvements in the models’ structure, states, or parameters; and (3) Hydrological Postprocessing (HPP), which can include bias correction techniques applied over the streamflow forecasts or the use of Bayesian approaches to estimate an optimal weighted combination of streamflow forecasts.

3. Results

The three main strategies for the evaluation of hydrological hypotheses (HPR, HP, and HPR) were applied to evaluate and improve the forecast skill of operational streamflow forecasts in the UZRB. The implementation of one or another strategy can have significant impacts on the final forecasting products. For example, precipitation bias correction methods (HPR strategies) can dramatically perturb the volume of water entering to the system, in the HP strategies the most sensitive model parameters can significantly modify water routing though the model structure, and HPP strategies can have a large impact by directly perturbing the forecasts to adjust scaling and fitting issues derived from a poor model representation. Details on how these three processing strategies can impact operational streamflow forecasts are discussed in the following sections of this paper.

3.1. Strategy 1: Hydrological Pre-Processing (HPR)

Rainfall forecasts or any other climate forecasts derived from satellites or climatological models are prone to errors that must be corrected. The propagation of these errors is more significant when the forecasts are biased, especially those for rainfall, since this is the most important variable for hydrological modeling. Corrections applied over rainfall records allow for the identification of sensitivities or gaps associated to the improvement of satellite data, for the calibration and validation of variational hydrological models, and for the identification or selection of the best ensembles for any operational implementation of an HFS.

The corrections correspond to the application of any selected HPR hypothesis, and they contribute to the reductions of the propagated meteorological uncertainties through the VEF implementation. One of the most popular HPR strategies is Quantile Mapping [54,55,56,57], a technique that has been previously evaluated in the UZRB and compared to Principal Components Analysis (PCA) [13]. In the operational context, the Quantile Mapping (QM) technique is applied at a daily time scale with the assumption that the probability density functions (PDFs) of the rainfall observations and forecasts follow Gamma PDFs. The key idea behind this technique is to swap the quantiles of the simulated data with the quantiles of the observed data. The application of this technique has shown that daily estimates from SPPs and RCMs can be satisfactorily corrected at the catchment scale or at more regional scales (i.e., [54,55,56,57,58,59,60]). The previous studies agree with our findings for the UZRB where precipitation estimates from three satellites were significantly improved after the application of the Quantile Mapping (QM) method, used with a Gamma Probability Distribution Function (PDF). The results showed that all raw SPPs (Figure 9a–c) could be satisfactorily corrected at daily time scales (Figure 9d–f). Here, it is important to notice that many missing and false detections of rainfall can be corrected (see all raw to corrected scatters in Figure 9); however, the selection and fitting of a fixed PDF to the whole rainfall dataset can also reduce the performance of rainfall forecasts in some areas if (1) the selected PDF is not a good representation of the local climate conditions; if (2) the application of QM is exclusively tied to regional parameters instead of cell-by-cell parameters or vice versa, and if (3) only fixed temporal parameters are used instead of temporally varying parameters. All these factors can have result on successful or inadequate rainfall corrections that can have a significant impact on the next steps related to model calibration, final structural design, and generation of final operational streamflow forecasts products.

3.2. Strategy 2: Hydrological Processing (HP)

The VEF hydrologic modeling paradigm requires bridging science and engineering for the design of functional and operational HFS. Therefore, all available precipitation datasets utilized during HPR, must be used to evaluate the propagation of meteorological uncertainty into final SR2MR streamflow forecasts (see Figure 10a₁–a₇). Generally, the objective of HP is to constrain the range of valid model outcomes for the application of HPP strategies and for the generation of forecast products. Multimodel Ensemble approaches are a popular alternative to propagate uncertainty into future forecasts; however, they can yield boundless errors in inference, which can produce unbounded uncertainty bands [61]. Although this is an acceptable argument, the science and engineering of streamflow forecasting can play an important role in defining the best alternatives for the management of total hydrologic uncertainty (THU). In doing this, Variational Ensemble Forecasting (VEF) approaches have the advantage that they can be implemented using all available sources of input data, hydrologic models, and optimal parameter sets, to improve the assimilation of forecasts. VEF can also be coupled with regularization techniques to constrain the forecasting range based on the classification and evaluation of historic events to define the best ensembles for HPP.

Constraining and selecting the best ensembles for operational forecasting should be understood as a procedure that can vary at any model run. The variation of ensembles depends on the historic forecast skill performance estimated for all available hydrological events and their classified characteristic responses. In the UZRB and its sub-basins, the implementation of an operational VEF approach was used to generate SR2MR streamflow forecasts derived from all possible combinations of SPPs, hydrologic models, and optimal parameter sets (Figure 10a₁–a₇). The propagation of meteorological uncertainty can be quite large in the resulting streamflow traces when the whole input-model-parameter space is mapped and used for operational forecasts. To avoid unbounded uncertainty bands, the VEF approach allows improving the accuracy of streamflow forecasts through a ranking evaluation and posterior identification of the best hydrologic ensembles for the UZRB (Figure 10b₁–b₇ for forecasts ranked using skill analysis of R² and Figure 10d₁–d₇ for forecasts ranked using total uncertainty defined as RMSE). Then the best ten raw streamflow forecasts are ranked using skill and/or uncertainty measures (Figure 10c₁–c₇,e₇–e₇). The application of specific or combined verification techniques on streamflow forecasts usually leads to an improvement of the HPP hypotheses. This multiple evaluation is helpful to identify how skill and uncertainty perform over space and time but also to evaluate how the spatial resolution of precipitation products can have a large effect on the operational forecasts. For instance, the smallest sub-basins inside the UZRB (see Kabombo, Kalabo, Chivata, and Kalene in Figure 10) resulted in higher total uncertainty compared to large basins (see Katima, Senanga, and Lukulu in Figure 10). This reduction in performance in the smallest sub-basins is probably associated with the smallest number of available precipitation grid cells, which produces larger averaged forecasting errors if missing or false detections are present in the forecasts.

3.3. Evaluating Pre-Operational SR2MR Streamflow Forecasts

All 72 possible SR2MR forecasts that were generated in the UZRB by the combination of multiple rainfall products, hydrological models, and optimal parameter sets, enabled the identification of differences between the satellite precipitation products (SPPs) and their HPR strategies (or bias corrections) but also between the hydrological models and their optimal parameter sets (or HP strategy). From the VEF approach, the best raw model ensembles (Figure 10) were retained (20 out of 72 ensembles were selected for this study) and used for the generation of final streamflow forecasts products and reports for end users.

In general, an expected outcome is to have climatology products with better forecast skill (see for example Figure 11a,b) because of the instrumental corrections applied during the Hydrological Preprocessing (HPR) stage. However, one disadvantage of the climatology products is that they are not available in the SR2MR domain. For this reason, SPPs and RCMs are still needed for SR2MR streamflow forecasts in the UZRB or any other basins around the World. Therefore, the selection of final streamflow forecasts is also operationally based on SR2MR records. Furthermore, the quantification and propagation of retrospective meteorological uncertainties might be required for the release of final forecasts depending on the needs of end users. In the UZRB, multiples ensembles for deterministic streamflow forecasts and the spread of uncertainty by means of a probability density function were established for each daily SR2MR forecast. These forecasts were provided for Namibian Hydrological Services, one of the relevant African institutions that manages water resources in the Upper Zambezi River.

3.4. Strategy 3: Hydrologic Post-Processing (HPP) for Raw Streamflow Forecasts

Our inability to generate exact physical representations of natural hydrologic systems creates the need for streamflow forecasts that can quantify and reduce the total hydrologic uncertainty (THU) propagated from a hydrologic modelling paradigm (HMP). Hydrological Postprocessing (HPP) hypotheses focus on establishing standardized methods to quantify and propagate total hydrological uncertainty (which is the sum of all uncertainties i.e., input, parameter, or structural uncertainties propagated into the final streamflow forecasts). Similar to HPR and HP, different methods can be applied for HPP hypotheses, i.e., stochastic, Bayesian, or machine learning methods can be used as post-processors for the final ensemble of SR2MR streamflow forecasts. The objective at this stage is to minimize the error of the deterministic forecast but also the spread of total hydrological uncertainty around the raw streamflow forecasts. The deterministic forecast can be obtained as an optimal weighted SR2MR streamflow forecast, which is used to propagate the total hydrological uncertainty (Figure 12). In the UZRB, a Multivariate Combinatorial Linear Regression (MCLR) approach was applied as a regularization technique to inform the best selection of hydrologic ensembles that minimized the spread of total hydrologic uncertainty. The MCLR was also used as a first level hydrologic postprocessor to resolve scaling issues of the raw streamflow forecasts (Figure 12). The final ensemble of raw and corrected SR2MR streamflow forecasts was then evaluated using two inference-based hydrologic postprocessors: (1) Inverse-Variance Weighting (IVW), and (2) Inverse-Probability Weighting (Figure 12). Both approaches were compared to the best streamflow prediction (BSP) and the Average Streamflow Prediction (AVSP) using the retrospective total hydrologic uncertainty quantified using the root mean square error in millimeters per day (Figure 12). The application of hydrologic postprocessors (MCLR, IVW, and IWP) revealed an improved efficiency of both the forecasts and the propagation of THU in the UZRB and its sub-basins (Figure 12).

The combination of HPP with short and long-term memory windows (for periods ranging between 5 and 180 days before the event) also proved to be more efficient in improving the performance of SR2MR streamflow forecasts (see Appendix A Figure A3 for details). The effect of HP for operational streamflow forecasting was mainly observed in two aspects: (1) a better performance of the deterministic streamflow forecasts, i.e., correction of scaling issues, and (2) narrower total hydrologic uncertainty bands around the deterministic forecast, i.e., more practical and realistic uncertainty bands for decision makers (see details in Figure 12).

4. Discussion

To establish operational hydrological forecasting systems (HFS) it is important to first define who will be the main and potential end users of the streamflow forecasts products within the basin under consideration. In the UZRB (Figure 1), the need for better water management schemes for flood warning and water allocation schemes has required the generation of daily SR2MR streamflow forecasts (Figure 2). Several water users and water authorities of the countries sharing waters from the UZRB can take advantage of the benefits of this operational HFS implementation, which can largely help developing sustainable water management and allocation activities for this transnational basin. For example, the water users of streamflow forecasts in the UZRB can be informed about the main water management decisions that need to be taken care of for water supply, reservoir management, hydraulic design, etc. These decisions are generally complex; therefore, under such a scenario, the best way to inform the local authorities and end users in the UZRB can be obtained if we had ‘perfect’ weather forecasts, or climatic predictions, that could be combined with ‘perfect’ hydrologic models to generate “almost” perfect streamflow forecasts. If this hypothesis were still true only one weather or climatic prediction (model), and one hydrologic model would be needed for streamflow forecasting. For example, the ‘perfect’ climate estimates of the climate model could be combined with one ‘perfect’ hydrological model, and the streamflow prediction should be close enough to inform both managers and users, so they can apply the right decisions for water management. The problem of this hypothesis is the fact that so far, hydrologists and meteorologists has not been able to create or establish perfect model representations of climate and hydrology, i.e., models are simplified representations of the hydroclimatological processes occurring across spatiotemporal scales. Therefore, we know and assume that during the operational implementation of the HFS in the UZRB there were countless sources of uncertainty, and different tools and schemes were established in this study to quantify how meteorological and hydrological uncertainty propagate through an operational VEF approach. One of the first adopted techniques that emerged was the utilization of single model realizations as those provided in many previous studies (see for example [3,4,5,6,7,8,9,10,11,12]). Then, multi-ensembles of climate predictions were used to quantify and propagate the meteorological uncertainty through the streamflow forecasts. The main issue of this technique is that it did not allow quantifying the uncertainty propagated from the hydrological model either from the model structure or from its parameters. To resolve this problem [15], proposed a way to quantify the hydrological uncertainty from a single model with multiple parameters that can be obtained from the calibration of multiple climate products. The final set of outputs is obtained from a multiproduct and multiparameter scheme. This technique allows quantifying the hydrological uncertainty propagated from multiple parameters sets of a single hydrological model; however, it does not allow quantifying how the uncertainty propagates and varies as a function of the structure of the selected hydrological model. This latter conceptualization of hydrological uncertainty is tied to the modelling assumptions required to establish the physical representation of the hydrologic system. Obviously, these assumptions also vary as a function of the selected hydrological model and its structure. In this context, hydrologists have opted to promote more comprehensive modelling schemes that combine multiproducts and multiple hydrological models [14,16,17,18]. With this technique, it has been possible to quantify the meteorological uncertainty propagated from the climate products and the hydrological uncertainty propagated from the hydrologic models and their corresponding structures. On the other hand, recent studies [62,63,64] have argued that the performance of a single hydrological model can be improved together with the streamflow forecasts if multiple climate products are combined with a model that is perturbed by generating multiple initializations with different initial conditions, i.e., changes in surface storage, groundwater storage, or snow storage, among others. The main objective behind this approach is to capture the current hydrologic condition of a catchment to assimilate streamflow. For example, in this regard it has been traditionally assumed that the warmup period is required for hydrologic modelling; however, this approach can be bypassed given that the initial states are adjusted beforehand, taking in consideration the streamflow assimilation. This scheme has also stablished a new source for the quantification of hydrologic uncertainty just by changing the initial conditions of the hydrological model. With this technique, a new research niche has been recognized by hydrologists and now is also applied to quantify total hydrological uncertainty (THU). Having said that, it is also important to add that all the schemes mentioned above have allowed quantifying the hydrological uncertainty in separated procedures, either from multiple climate products, from multiple hydrological models, from multiple parameters sets, or from multiple initializations. In fact, all these hydrologic modeling paradigms (HMP) are still applied in a systematic manner to identify and quantify different sources of hydrological uncertainty.

Despite the existence of all these HMP previously mentioned, none of them have proposed a combined technique to quantify both meteorological and hydrological uncertainties propagated from different sources of a VEF implementation (see Figure 7 and Figure 8), i.e., sources as climate forecasts, modelling assumptions, and optimal parameter sets, that can be evaluated for any operational hydrological forecasting system (HFS) implementation with VEF. To accomplish this issue, we have implemented a VEF approach [22] based on multiple satellite precipitation products (and GFS precipitation forecasts), multiple hydrologic models, and multiple optimal parameters sets for SR2MR daily streamflow forecasting. The VEF approach implemented in the UZRB (Figure 10) has allowed increasing the number of possible hydrologic ensembles available for streamflow forecasting, together with an improvement of the streamflow assimilation (observed versus predicted). It also provides a more comprehensive and systematic framework to identify and propagate the spread of total hydrologic uncertainty in an operational hydrologic forecasting system. Now, the natural question is how much room is left to define new modelling paradigms or techniques that can be used to quantify or minimize the propagation of total hydrological uncertainty? To answer this question, we need to differentiate between what we can do to define the best integrated implementation of an HFS, and what adjustments are required at each separated component (i.e., inputs, models, or outputs) of the HFS. The operational HFS implementation in the UZRB, identified and considered three general hydrologic processing strategies (HPS) that can be applied to any VEF approach. If we consider that any hydrologic modelling paradigm (HMP) can either include multiple inputs; hydrologic models, parameters, initializations; and outputs, then, we can hypothesize that at each component of the HMP it is possible to apply additional techniques or methods to improve streamflow forecasting. These strategies were conceptualized as: (1) Hydrological Preprocessing (HPR), (2) Hydrological Processing (HP), and (3) Hydrological Postprocessing (HPP). Taking advantage of the strategies proposed for the establishment of hydrologic processing strategies or hypotheses, standardized methods for HPR (Figure 9) and HPP (Figure 12) were applied in the UZRB, showing that the performance of raw VEF streamflow forecasts can be significantly improved, and the spread of uncertainty can be better constrained by applying regularization processes that combine the strength of Multivariate Combinatorial Linear Regression (MCLR) and Inference-Based approaches (Figure 12).

The science and engineering of future operational streamflow forecasting in the UZRB will continue concentrating efforts in improving the forecasts; however, new needs from end users might also require improving the physical representation at the catchment scale. For this, it will be necessary to establish the role of mathematical, statistical, and/or machine learning methods that can be used to correct and propagate the hydrologic uncertainty from different components of a VEF approach, i.e., Bayesian, stochastic, pattern, or inference-based learning methods, etc. The performance of different rainfall products and methods also needs to be evaluated across different catchment sizes and forecast timescales to determine the space–time variability that propagates total hydrologic uncertainty. This evaluation must also be extended to determine the dependence of the accuracy of the streamflow forecasts and the propagation of hydrological uncertainty from physically based models. All the above will indicate the applicability limits of the VEF approach based on multiple precipitation or climate products, multiple hydrologic models, and multiple optimal parameters sets. All of the above will allow the establishment of new robust theoretical and hypothetical paradigms to quantify, evaluate, reduce, and manage the propagation of total hydrologic uncertainty using VEF approaches.

5. Conclusions

This paper described the main stages and processes required to implement and improve an operational hydrologic forecasting system (HFS) in the UZRB and its sub-basins. The process of implementation is very complex, and important decisions needed to be made about the input data (precipitation from satellites or climate products), the hydrologic models to be included along with their optimal parameter sets, and the timescales required for the generation of streamflow forecasts.

Once the HFS was completely operational in the UZRB, additional improvements to the forecasts were required to improve its performance and reduce the spread of total hydrologic uncertainty into the final streamflow forecast products. In this regard, three general strategies to improve the performance of VEF approach were proposed: from Hydrological Preprocessing to Postprocessing techniques that can improve the input data, the hydrologic models (or their structures), the optimal parameter sets, and the raw streamflow forecasts. The whole range of available techniques for operational HFS will require more detailed and standardized conceptualizations. In this regard, bias corrections or preprocessing (HPR) techniques applied over the input data will still play an important role in operational hydrological forecasts for the UZRB. The operational implementation of the VEF combined with regularization and inference-based methods improved the performance of streamflow forecasts as the primary need from end users in the UZRB; however, new alternatives to improve the physical understanding of the basin are still a pending task. Finally, it is important to add that the science of Hydrological Postprocessing is still under an early stage of development, and it still lacks the standardized methods that can be used for these purposes. Emerging methods will need to be evaluated to establish the real boundary between physical and statistical needs in operational streamflow forecasting. This discipline will also require a natural merging of science and engineering (practical applications) in a real-world context to establish baseline conditions for streamflow forecasting and hydrological uncertainty quantification and propagation.