Verification of short-term runoff forecasts for a small Philippine basin (Marikina)

ABSTRACT

Storm runoff from the Marikina River Basin frequently causes flood events in the Philippine capital region Metro Manila. This paper presents and evaluates a system to predict short-term runoff from the upper part of that basin (380 km²). It was designed as a possible component of an operational warning system yet to be installed. For the purpose of forecast verification, hindcasts of streamflow were generated for a period of 15 months with a time-continuous, conceptual hydrological model. The latter was fed with real-time observations of rainfall. Both ground observations and weather radar data were tested as rainfall forcings. The radar-based precipitation estimates clearly outperformed the raingauge-based estimates in the hydrological verification. Nevertheless, the quality of the deterministic short-term runoff forecasts was found to be limited. For the radar-based predictions, the reduction of variance for lead times of 1, 2 and 3 hours was 0.61, 0.62 and 0.54, respectively, with reference to a “no-forecast” scenario, i.e. persistence. The probability of detection for major increases in streamflow was typically less than 0.5. Given the significance of flood events in the Marikina Basin, more effort needs to be put into the reduction of forecast errors and the quantification of remaining uncertainties.

KEYWORDS:

• short-term flood forecasting
• weather radar
• quantitative precipitation estimation
• Philippines

Editor Z.W. Kundzewicz; Associate editor G. Di Baldassarre

1 Introduction

Runoff forecasting in small basins is a challenge for a number of factors. Compared to large basins, runoff predictions in small catchments are generally more sensitive to external forcing, i.e. rainfall. This is due to short drainage paths and limited storage capacities, resulting in short response times. At the same time, hydrological predictions are more vulnerable to spatially inaccurate estimates/forecasts of rainfall. While such errors may balance out at larger scales, this is not the case in small basins (Arnaud et al. 2011, Vieux and Imgarten 2012). Finally, the probability of extreme areal precipitation generally increases as the considered area becomes smaller (see, e.g. Rodriguez-Iturbe and Mejía 1974, Frederick et al. 1977).

Runoff prediction involves two major aspects: (1) estimation/forecasting of precipitation; and (2) modelling of the basin’s hydrological response. The hydrological models employed may be either empirical or process-based in nature. While the first class of models (see overview in Solomatine and Ostfeld 2008) is easy to calibrate, simple to use, and computationally efficient, empirical models are criticized for their black-box character. Outputs may become unrealistic and even unphysical in extrapolation mode, i.e. when the model is confronted with input signals not contained in its training dataset. This is a particular problem for the prediction of rare extremes, e.g. floods (see Schmitz and Cullmann 2008).

In this paper, the focus is on the second class of models built on the mathematical representation of hydrological processes. These models combine the principle of mass conservation with more or less simplified versions of the energy balance equation. Process-based models are commonly used for the purpose of hydrological forecasting (Kitanidis and Bras 1980, De Roo et al. 2000, Plate 2009). Of special interest are spatially distributed models that discretize river basins into smaller sub-units, e.g. sub-basins. This is a prerequisite for taking advantage of high-resolution precipitation fields as produced, for example, by radar-based quantitative precipitation estimation (QPE). Nonetheless, discretization is necessary to represent spatial heterogeneities in terms of states and parameters and, thereby, it is a key to the proper simulation of the dynamic nonlinearities (Robinson et al. 1995) occurring in rainfall–runoff transformation.

The kind of precipitation data required for runoff forecasting depends on both the intended maximum lead time (LT_max) as well as the basin’s response time (RT). As long as LT_max ≤ RT, real-time QPE data (i.e. recent observations) are sufficient. If LT_max > RT, the hydrological model must also be fed with a quantitative precipitation forecast (QPF). In cases where LT_max $\gg$ RT, the QPF is typically provided by numerical weather prediction models (e.g. Kneis et al. 2012). If LT_max exceeds RT by only minutes to a few hours (at most) the required short-term QPF, also referred to as a nowcast, often originates from extrapolation methods (Toth et al. 2000, Pierce et al. 2004, Ehret et al. 2008). Hybrid approaches to short-term QPF that combine dynamic and statistic rainfall forecasts (e.g. Golding 1998) are also in use (Werner and Cranston 2009).

This paper demonstrates the performance of a short-term hydrological forecast system for a small basin exploiting the limited information that is typically available in a newly industrialized country. The case study was motivated by the remarkable frequency and magnitude of storm events hitting the Philippines in recent years. The basin of interest (380 km², RT ≈ 3 h) is located just upstream of the eastern part of Metro Manila. It has been selected for its high-risk status, i.e. the exceptional loss potential in combination with the frequent occurrence of extreme runoff events.

Hourly deterministic forecasts were generated for two rainy seasons using a process-based hydrological model in conjunction with pragmatic but robust approaches to streamflow data assimilation. The model was forced with QPE data only, i.e. no precipitation forecast was used. Therefore, the analysis is limited to LT_max ≈ RT. Since weather radar stations were recently installed in the Philippines, the possible benefit of radar-based precipitation estimates over a conventional, raingauge-based QPE was of particular interest. The specific value of this study lies in the fact that the alternative QPEs are not just evaluated with respect to the quality of hydrological simulations but with respect to forecasts where deficits in QPE are potentially alleviated through data assimilation.

The major aim of this case study was to evaluate the quality of the hydrological forecasts for the purpose of flood warning and to provide a roadmap for further improvement of the forecast system. Ideally, the study would serve as a precursor to the establishment and verification of runoff forecasting systems for other high-risk basins across the Philippines in the context of the Nationwide Operational Assessment of Hazards (http://noah.dost.gov.ph). It aims especially at fostering the use of data from the country’s recently installed network of weather radar stations. Last but not least, the study demonstrates the productive use of open-source software for radar data processing (Heistermann et al. 2015) and hydrological modelling (Kneis 2015), enhancing transparency and reproducibility of the methodological approach.

The terms simulation, prediction, forecast and hindcast appear several times in the subsequent text. These words are used according to the following convention: (1) Simulation refers to the situation where the hydrological model is fed with observed meteorological forcings to estimate streamflow of the past. (2) The terms forecast and prediction are used synonymously. They refer to the case where the model is used to estimate future streamflow. In an operational context, the model must then be fed with predicted meteorological forcings. (3) The term hindcast denotes a retrospective forecast experiment, carried out with the intention of comparing the predicted outcome of an event to the actual (known) outcome. This is the usual approach to forecast verification.

2 Study area

The study area (Fig. 1) comprises the Marikina River Basin (MRB) up to the gauging station at Montalban (380 km²). The Marikina River and its tributaries originate in the mountain range northeast of the Philippine capital Manila. While the upper parts of the basin are covered by forests, the lower parts are dominated by agricultural use. Urban areas are located near the Montalban gauging station and further downstream.

Figure 1. Upper Marikina Basin with topography (ASTER elevation model), drainage network and sub-basin boundaries. River and raingauges are marked by triangle and squares, respectively. Dashed curves indicate the distance to the radar station at Subic (Table 1). Coordinates are in WGS84

As typical for this part of the Philippines, the climate of MRB is characterized by wet conditions from May to October and a dry season lasting from November to April. During longer periods of drought, streamflow drops to near zero values. Average annual precipitation amounts to 2400 mm. Flood-inducing storm events are often related to typhoons. Strong activity of the southwest monsoon and convective rainfall in the Inter-Tropical Convergence Zone (ITCZ) are relevant as well. Peak discharges of the 10 largest floods for 1959–2000 were reported in Abon et al. (2011), ranging from about 1600 to 2600 m³ s⁻¹ with reference to a catchment area of 540 km². Severe floods were observed in September 2009 (Abon et al. 2011), August 2012 (Heistermann et al. 2013a), August 2013 and more recently in September 2014.

Several streamgauges exist in the MRB but most of them do not provide continuous observations of adequate temporal resolution or the water levels are affected by backwater. Hence, only data from the Montalban gauge (Fig. 1) were used in this study. A typical lag between onset of rainfall in the upstream catchment and rising stage at this gauge is 3 hours.

3 Data and methods

3.1 Radar-based precipitation estimates

The Subic radar station is located about 90 km west of the study area (Fig. 1) at 14.822°N, 120.363°E. The station has been operational since March 2012. Its technical features are summarized in Table 1.

Table 1. Features of the weather radar station near Subic

Bandwidth	S-band
Polarization	Single
Maximum range	120 km
Gate length	500 m
Azimuth resolution	1°
Elevation angles	14 (0.5°–19.5°)
Volume cycle interval	9 min

The Subic radar data were processed using routines from the open-source software library wradlib (Heistermann et al. 2013b). In the first step, a static clutter map was generated by identifying bins with reflectivity persistently greater than 15 dBZ in spite of clear-air conditions. The original reflectivity in those bins was substituted with data from horizontally neighbouring bins not contaminated with clutter.

Second, reflectivity (dBZ) was corrected for a systematic and persistent underestimation (by a factor of 6 dB). The level and persistence of this underestimation were verified by raingauge observations and were presumably caused by hardware miscalibration. This single-parameter correction needs to be distinguished from actual radar-gauge adjustment (e.g. Goudenhoofdt and Delobbe 2009), where parameters typically vary in space and/or time.

In the third step, reflectivity (Z; mm⁶ m⁻³) was transformed into rainfall rate (R; mm/h) using Z = 250R ^1.2. The same relation is applied by the US national weather service NOAA for tropical cyclone conditions.

Finally, a constant altitude plan position indicator (Pseudo-CAPPI) was generated at 2000 m above sea level based on the sweeps at elevation angles of 1.5° and 2.4° (data from the 0.5° angle had to be discarded due to the study basin’s mountainous terrain and the resulting level of ground clutter contamination). Rainfall intensities at the reference altitude were then interpolated to centroids of sub-basins (Fig. 1) as an input for the hydrological model.

The radar-based rainfall estimates were compared with observations at raingauges by Abon et al. (2015). In spite of inherent problems (see, e.g. Kitchen and Blackall 1992) such a comparison can detect major deficits in the radar data and it is a useful complement to hydrological validation (Heistermann and Kneis 2011). According to Table 2, the radar-based data moderately underestimate observed rainfall in the Marikina Basin. In addition to the negative bias, considerable random errors are present. While the radar-based estimates correspond well with gauge data at the daily scale (Nash-Sutcliffe index, E > 0.8), the correlation at the hourly level is rather weak (≤ 0.4). An increased mismatch at shorter time intervals has been found in other studies (Seo and Krajewski 2010, Sebastianelli et al. 2013), and this may be explained partly by the different sampling volumes of raingauge and radar.

Table 2. Error in radar-based precipitation estimates with reference to raingauge observations in the Marikina Basin. E denotes the coefficient of efficiency (Nash-Sutcliffe index)

Interval	Year	Bias (%)	RMSE (mm)	E (–)
Daily	2012	−8.5	15.2	0.81
	2013	−8.4	15.1	0.88
Hourly	2012	−16.2	8.7	0.32
	2013	−7.8	3.8	0.40

3.2 Hydrological model

3.2.1 Model engine

The hydrological model engine used is called HYPSO-RR. It is an object-oriented software built with the ECHSE framework (Kneis 2015). HYPSO-RR is designed for the continuous simulation of streamflow in mesoscale river basins. It belongs to the class of conceptual model engines that combine both physically-based and empirical equations. Since HYPSO-RR has been developed mainly for forecasting purposes, the mathematical representation of hydrological processes was kept simple for the benefit of computational efficiency. This is essential because data assimilation (namely state updating) and uncertainty estimation (e.g. using precipitation ensembles) require a large number of repeated model runs.

A short summary of the classes is contained in Kneis (2015, Table 4). A comprehensive documentation of the mathematical formulations of hydrological processes can be found in Kneis (2012a) Subsequently, only the two most important classes are outlined.

3.2.2 Sub-basin class

Objects of this class simulate the water balance at sub-basin level. At present, only three basic types of land use are distinguished: vegetated soil, impervious surfaces and water. Runoff generation on vegetated areas was computed from the water balance of a homogeneous soil column. Direct runoff (surface runoff, sub-surface stormflow) was assumed to occur on saturated zones only. The extent of those zones was estimated using the method of Zhao et al. (1980) and the rate of direct runoff was computed analytically following Todini (1996, signs corrected). The slower components (interflow, groundwater recharge) were modelled as functions of the soil water content as in LARSIM (Ludwig and Bremicker 2006). The equations to simulate the energy and water balance of a possible snow cover were mostly adopted from Tarboton and Luce (1996).

Potential evapotranspiration (ET) was estimated from short-wave (gross) radiation and temperature data with Makkink’s approach De Bruin (1987). Actual ET was obtained by multiplying the potential rate with a crop factor (see, e.g. Feddes 1987) and a linear term accounting for moisture limitation. The leaf area index was used as a proxy for the crop factor and a linear correlation was assumed between the two quantities. This simple ET model is primarily attractive for the minimum set of predictor variables required.

Storage and concentration of runoff at the sub-basin level were modelled as in LARSIM (Ludwig and Bremicker 2006). Thus, the generated runoff was routed through a single linear reservoir to emulate the distribution of travel times. This was done separately for the individual runoff components. The retention parameters of the linear reservoirs were derived from a concentration time index (CTI) which accounts for the sub-basin’s size, shape and topography. The CTI was then scaled with dimensionless factors identified from recession analysis and subsequent calibration.

The set of state variables considered in the sub-basin class comprise (1) the soil water content, (2) the storage volumes of the linear reservoirs controlling runoff concentration, as well as (3) snow-related quantities not relevant in this study’s context. The required set of meteorological forcings includes precipitation, air temperature, short-wave radiation, air pressure, wind speed, relative humidity and cloudiness. For regions without snow cover, only the first four variables are relevant. To allow for seasonal variation, the leaf-area index was treated as an external forcing too.

3.2.3 Reach class

In HYPSO-RR, a river reach is modelled as a single linear reservoir (LR) whose retention constant is variable. The LR approach was adopted for its attractive analytical solution. The particular solution presented in equation (1(1) $\begin{array}{rl}& v(t_0+\mathrm{\Delta }t)=v(t_0)x+\frac{a(x-1)}{b^2}\\ & +\frac{q_{\mathrm{i}\mathrm{n}}(t_0)(x-1)-a\mathrm{\Delta }t}{b}\end{array}$(1) ) allows for a linear variation of the inflow rate q _in during a discrete time step of length ∆t.

(1) $\begin{array}{rl}& v(t_0+\mathrm{\Delta }t)=v(t_0)x+\frac{a(x-1)}{b^2}\\ & +\frac{q_{\mathrm{i}\mathrm{n}}(t_0)(x-1)-a\mathrm{\Delta }t}{b}\end{array}$(1)

with

$\begin{array}{rl}& a=(q_{\mathrm{i}\mathrm{n}}(t_0+\mathrm{\Delta }t)-q_{\mathrm{i}\mathrm{n}}(t_0))/\mathrm{\Delta }t\\ & b=-1/k\\ & x={\mathrm{e}}^{(-\mathrm{\Delta }t/k)}\end{array}$

In a LR, the outflow rate (q _ex) is linear function of storage (v), and the retention constant (k) is given by k = v/q _ex. For real channel sections this relation is generally nonlinear. However, in order to take advantage of the analytical solution (equation (1(1) $\begin{array}{rl}& v(t_0+\mathrm{\Delta }t)=v(t_0)x+\frac{a(x-1)}{b^2}\\ & +\frac{q_{\mathrm{i}\mathrm{n}}(t_0)(x-1)-a\mathrm{\Delta }t}{b}\end{array}$(1) )), a locally linear behaviour (k ≈ ∆v/∆q) is assumed. The curve v(q) can be inferred from information on the cross-section’s geometry, bed slope and roughness using Manning’s equation; thus, empirical functions k = f(v) and k = f(q) are available. They allow for the selection of an appropriate k for use with equation (1(1) $\begin{array}{rl}& v(t_0+\mathrm{\Delta }t)=v(t_0)x+\frac{a(x-1)}{b^2}\\ & +\frac{q_{\mathrm{i}\mathrm{n}}(t_0)(x-1)-a\mathrm{\Delta }t}{b}\end{array}$(1) ), reflecting the reach’s current state v(t ₀) as well as its inflow q _in at times t ₀ and t ₀ + ∆t.

3.3 Model set-up and calibration

Drainage network and catchment boundaries were derived from the ASTER global digital elevation model (DEM), a product of METI and NASA, using software described in Kneis (2012b). The catchment area (Fig. 1) was divided into sub-basins with an average size of 2 km². Information on land use was gathered from the National Mapping and Resource Information Authority (NAMRIA). Data on the spatial distribution of soil types and average texture information were provided by the Philippines Bureau of Soils and Water Management. The hydraulic properties of river reaches were derived from survey data of 42 cross-sections distributed over the basin. First, power models r(h) and a(h) were fitted for the individual survey cross-sections with h, r and a denoting flow depth, hydraulic radius and wet area, respectively. Second, r(h) and a(h) were calculated for every modelled reach as a weighted average of the respective curves at two survey cross-sections. The particular survey cross-sections were selected according to distance to the reach (to be minimized) and their upstream catchment areas (upstream area of the reach to be within the range of the two cross-sections’ upstream areas). The r(h) and a(h) functions obtained provide the basis for inferring the curve v(q) (see Section 3.2.3) by Manning’s equation.

The hydrological model was calibrated using hourly streamflow data from the Montalban gauge (Fig. 1). Model runs were initialized at least 10 years before start of the actual calibration period and the outputs for this warm-up period were discarded. The initial state for the warm-up period was generated in a long-term simulation repeatedly using meteorological data for the years 2004–2011 as input.

The calibration of nine parameters was carried out by a supervised stochastic approach (Kneis et al. 2012, 2014) characterized by the following steps:

• Sampling ranges are specified for all parameters. Plausible ranges can be deduced from basic soil properties or hydrograph analysis (Kneis 2012a) or previous model applications.

• Random parameter sets are drawn from within these ranges by latin hypercube sampling (Carnell 2012) and the model is run for all sets (typically 200–1000 runs).

• For all runs, the objective function is computed and plotted against the corresponding parameter values (one plot for each parameter).

• The plots are visually inspected and the sampling ranges for sensitive parameters are manually adjusted, i.e. narrowed or shifted.

• Steps 2–4 are repeated until further modification of the sampling ranges no longer results in an improvement of the objective function after reasonable rounding (e.g. considering three significant digits).

In this study, the Nash-Sutcliffe index (Nash and Sutcliffe 1970) was used as the primary objective function. Being a rescaled sum of squared residuals, the Nash-Sutcliffe index naturally puts the focus on peak flows (Krause et al. 2005). Information on the bias was used to confine the sampling range (Step 4 above) in cases where the primary objective function was insensitive to changes in the respective parameter. Parameters that control groundwater recharge and drainage were fitted by analysing the model error during periods of low flow.

The model was calibrated separately for (1) precipitation information from raingauges; and (2) radar-based precipitation estimates (see Section 3.4.2). This was considered necessary since errors in forcings may be compensated partly by parametrization (see, e.g. Gourley and Vieux 2005, Heistermann and Kneis 2011). In the case of gauge-based rainfall, the calibration dataset covered the years 2004–2011; hence, there was a proper distinction between the period of parameter estimation and the period of forecast verification (2012–2013). A Nash-Sutcliffe index of 0.67 was obtained for the entire period 2004–2011. The annual values of Nash-Sutcliffe index ranged from 0.42 to 0.85 and the bias varied between about $-$25 and +25%. The model performed worst in 2008 and 2010 due to underestimation of peaks.

In contrast to raingauge data, radar-based precipitation estimates were available for the years 2012 and 2013 only, making a distinction between calibration and verification periods impossible. Consequently, some artificial skill was to be expected in the radar-based forecasts, resulting in a too optimistic assessment of forecast performance. In order to examine the hypothesised effect, a cross-check was carried out by feeding the model with radar rainfall while using parameters optimized for the gauge-based QPE. Contrary to expectation, the substitution of model parameters did not result in a relevant drop in the Nash-Sutcliffe index. Hence, the acceptable fit of the radar-driven simulations (Nash-Sutcliffe index of 0.79, bias −0.3%) is unlikely to result from parameter (over)fitting and, therefore, artificial skill in the radar-based forecasts appears to be less of a concern.

Simulated flow rates based on either raingauge or radar data are plotted in Fig. 2 together with observations.

Figure 2. Observed and simulated streamflow for the two major events in 2012 and 2013

3.4 Generation of forecasts

3.4.1 Time frame and basic concepts

The study follows the approach of retrospective forecast verification known as hindcasting. That is, runoff forecasts are generated for a series of events in the past and predictions are then compared to the known outcome, i.e. observed flow data.

The verification experiment covered the two rainy seasons 10 June 2012–15 October 2012 and 6 June 2013–10 September 2013 according to the availability of radar data. During these periods, a new forecast was computed every hour unless streamflow data were missing/invalid at the respective time. The total number of verifiable forecasts amounts to 3810.

The quality of hydrological forecasts strongly depends on adequate initialization of the model’s state variables at the time of forecasting. Those initial values cannot be measured because the state variables often do not strictly correspond to observable quantities (e.g. filling of conceptual reservoirs). Even if they do (e.g. in the case of soil moisture), it is difficult to obtain real-time data with sufficient spatial resolution (see Bronstert et al. 2012). Hence, proper initial states must be generated with the model itself by means of a “warm-up” run (leftmost arrow in Fig. 3). The warm-up run is typically a long-term simulation based on observed meteorological forcings. The warm-up period must be chosen long enough to eliminate the effect of a purely guessed initial state.

Figure 3. Outline of a single forecast experiment. Arrows indicate individual runs of the hydrological model. T0 and LT_max denote the forecast initialization time and the maximum lead time, respectively

The set of arrows in the centre of Fig. 3 indicates repeated simulations carried out in the context of sequential data assimilation (DA, see Section 3.4.3 for details). The DA aims to further refine values of state variables obtained from the warm-up run. Specifically, DA attempts to optimally estimate the values of state variables at time T0 taking into account recent observations.

The rightmost arrow in Fig. 3 finally represents the actual forecast run initialized at T0. Since the model is simulating into the future, it must be fed with predicted meteorological forcings in that case (see Section 3.4.2).

3.4.2 Meteorological forcings

In simulation mode, i.e. for all times prior to T0, the hydrological model was forced with meteorological observation data. In the case of rainfall, two alternative variants of QPE were tested. In the first variant, denoted GG, hourly observations from the raingauges shown in Fig. 1 were interpolated to the sub-basin centroids by inverse distance weighting (IDW). Terrain was not used as background information since data did not show a close correlation between rainfall and elevation. In the second variant (hereafter denoted RR), the radar-based precipitation estimates according to the specifications in Section 3.1 were used. Observations of other meteorological variables (Section 3.2.2) were also regionalized by IDW.

The forcings used to drive the hydrological model in forecast mode, i.e. beyond T0, were chosen on a pragmatic basis. Regarding rainfall, replicates of the latest known values (i.e. those corresponding to time T0) were used as estimates for future time intervals. This implies stationarity of the precipitation field as well as persistence with respect to rainfall intensities. The assumption was applied to both variants of QPE, i.e. GG and RR. For all other meteorological variables monthly average values were used as dummy forecasts.

The pragmatic choice with respect to meteorological forecasts is justified by the study’s emphasis on short-term predictions. Consequently, the discussion in Section 5 focuses on lead times being not much greater than the basin’s response time (recall Section 1). Experiments targeted at extending the forecast range by means of rainfall nowcasting are currently in preparation.

3.4.3 Data assimilation

Data assimilation (DA) aims to estimate a system’s true state by combining recent observation data with model-generated information. An introduction to DA in the context of hydrological forecasting can be found in Refsgaard (1997). Recent developments and challenges are reviewed in Liu et al. (2012). Many techniques of DA (e.g. state updating through filtering or optimization methods) are targeted at optimal initialization of the model at the forecast time (T0). In contrast, error prediction (EP) refers to the correction of the forecasts (i.e. model outputs) based on the statistics of recent simulation errors. From a theoretical point of view, EP is clearly inferior to updating techniques since predicted flow rates are related to the catchment’s initial state in a nonlinear way. Nevertheless, EP is attractive for practical applications as it is simple to implement and cheap in terms of computational cost.

In this study, we generally relied on EP to merge model predictions with recent observational data. The correction scheme (equation (2)(2) ${\hat{p}}_k=\left(\begin{array}{cc}{p}_k\frac{o_r}{s_r}& \mathrm{i}\mathrm{f}{p}_k<s_r\\ p_k+\left(o_r-s_r\right)& \mathrm{i}\mathrm{f}{p}_k\ge s_r\end{array}\right.$(2) ) was adopted from an operational model commonly used in Germany (Ludwig and Bremicker 2006, pp. 78–79). In the following equation, ${\hat{p}}_k$ is the corrected forecast at lead time with index k and p_k is the corresponding “raw” forecast; symbols o_r  and s_r denote the observed and simulated flow rates just before forecast initialization:

(2) ${\hat{p}}_k=\left(\begin{array}{cc}{p}_k\frac{o_r}{s_r}& \mathrm{i}\mathrm{f}{p}_k<s_r\\ p_k+\left(o_r-s_r\right)& \mathrm{i}\mathrm{f}{p}_k\ge s_r\end{array}\right.$(2)

The distinction between the two cases of equation (2(2) ${\hat{p}}_k=\left(\begin{array}{cc}{p}_k\frac{o_r}{s_r}& \mathrm{i}\mathrm{f}{p}_k<s_r\\ p_k+\left(o_r-s_r\right)& \mathrm{i}\mathrm{f}{p}_k\ge s_r\end{array}\right.$(2) ) is important. The multiplicative approach (first case) guards against negative values in the case of predicted flow recession. The additive correction (second case) avoids excessive amplification of predicted floods in situations where o_r /s_r is large while the absolute error is small, e.g. during low flow. In its present form, equation (2(2) ${\hat{p}}_k=\left(\begin{array}{cc}{p}_k\frac{o_r}{s_r}& \mathrm{i}\mathrm{f}{p}_k<s_r\\ p_k+\left(o_r-s_r\right)& \mathrm{i}\mathrm{f}{p}_k\ge s_r\end{array}\right.$(2) ) is applicable to short-term predictions only. An extension of the equations may be necessary in the case of medium- and long-term forecasts, so as to allow for a gradual attenuation of the correction terms over the forecast period.

In addition to EP, combined input and state updating of the hydrological model was applied. In this study, rainfall (observed input) as well as storage in the linear reservoirs controlling sub-basin runoff (state variables; see Section 3.2.2) were updated. The approach belongs to the class of non-sequential DA methods; hence, the goal was to minimize the discrepancy between observed (o) and simulated flow rates (s) in a time period prior to T0. This time period is subsequently called the assimilation period (AP; see Fig. 3). Its length was set to 24 hours. To measure the discrepancy between o and s, a weighted sum of squared errors (wSSE, equation (3(3) $\mathrm{w}\mathrm{S}\mathrm{S}\mathrm{E}=\sum _{i=1}^n\left({\left(o_i-s_i\right)}^2{w}_i\right)$(3) )) was adopted as the objective function. The weights (w) were chosen to gradually increase over the AP from w ₁ = 0 to w_n  = 1, in order to give greater weights to observations close to T0 (n: number of hourly observations; equal to 24). The length of the AP and the weights were fixed based on empirical tests.

(3) $\mathrm{w}\mathrm{S}\mathrm{S}\mathrm{E}=\sum _{i=1}^n\left({\left(o_i-s_i\right)}^2{w}_i\right)$(3)

In the given context, rainfall input during the AP and the initial catchment storage at T0 − 24 h represent the parameters to be optimized when minimizing equation (3(3) $\mathrm{w}\mathrm{S}\mathrm{S}\mathrm{E}=\sum _{i=1}^n\left({\left(o_i-s_i\right)}^2{w}_i\right)$(3) ). Varying these quantities directly is rather inconvenient, not least because rainfall is a time series rather than a scalar. Therefore, the optimization was performed with respect to dimensionless multipliers applied to the original values of rainfall and initial catchment storage. Optimization was generally started with an initial guess of 1 for all multipliers; hence, it was assumed that updating of input and states is not actually necessary. Minimization of the objective function with respect to the multipliers was achieved by stepwise line search. The original algorithm (Kuzmin et al. 2008) was adapted to speed up convergence in cases where the objective function is insensitive with respect to a parameter (i.e. a multiplier). To prevent very strong corrections, optimization was constrained. The multipliers’ allowed ranges were set to 0.5–1.5 for rainfall and 0.75–1.25 for initial catchment storage. This is equivalent to maximum corrections by ±50% and ±25%, respectively.

Note that error prediction can easily be combined with model updating. In that case, symbol s_r in equation (2(2) ${\hat{p}}_k=\left(\begin{array}{cc}{p}_k\frac{o_r}{s_r}& \mathrm{i}\mathrm{f}{p}_k<s_r\\ p_k+\left(o_r-s_r\right)& \mathrm{i}\mathrm{f}{p}_k\ge s_r\end{array}\right.$(2) ) represent the optimized hydrograph after updating of states and inputs.

3.5 Verification

Skill scores (equation (4(4) $\mathrm{S}\mathrm{k}\mathrm{i}\mathrm{l}\mathrm{l}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\frac{\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}}-\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}}{\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}}-\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}}$(4) )) rate the quality of a forecast in relation to a reference forecast. They are particularly useful if the performance of alternative forecasts is to be compared.

(4) $\mathrm{S}\mathrm{k}\mathrm{i}\mathrm{l}\mathrm{l}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\frac{\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}}-\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}}{\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}}-\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}}$(4)

For a perfect forecast, the skill score takes a value of 1 (upper bound). A value of 0 indicates “no skill” as compared to the reference. The lower limit depends on the choice of the actual score used to rate the individual forecasts. In this application, the mean squared error (MSE) was chosen, hence the first term in the denominator of equation (4(4) $\mathrm{S}\mathrm{k}\mathrm{i}\mathrm{l}\mathrm{l}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\frac{\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}}-\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}}{\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}}-\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}{\mathrm{e}}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}}$(4) ) was set to zero. The resulting skill score (lower bound −∞) is also known as the reduction of variance (RV, equation (5(5) $\mathrm{R}\mathrm{V}=\frac{\mathrm{M}\mathrm{S}{\mathrm{E}}_{\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}}-\mathrm{M}\mathrm{S}{\mathrm{E}}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}}{-\mathrm{M}\mathrm{S}{\mathrm{E}}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}}$(5) )). Identical values of RV were obtained if MSE was substituted by the Nash-Sutcliffe index and the score of the perfect forecast was set to 1, accordingly.

(5) $\mathrm{R}\mathrm{V}=\frac{\mathrm{M}\mathrm{S}{\mathrm{E}}_{\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}}-\mathrm{M}\mathrm{S}{\mathrm{E}}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}}{-\mathrm{M}\mathrm{S}{\mathrm{E}}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}}$(5)

In this study, persistence was used as the reference forecast. It is the natural first choice for short-term predictions because streamflow data exhibit strong auto-correlation for lag times not much greater than the basin’s response time. In other words, recent observations are already valuable predictors for future conditions.

In order to assess the potential value of the deterministic runoff forecasts in absolute numbers, the probability of detection, POD, and the false alarm ratio, FAR (Jolliffe and Stephenson 2003), were also evaluated. In order to apply POD and FAR to forecasts of continuous variables, the latter need to be converted to binary variables by means of a threshold value X (see, e.g. Bürger et al. 2009). Then, POD represents the fraction of cases where an observed exceedence of X was actually forecast (optimum value 1). FAR represents the fraction of cases where an exceedence of X was predicted but not observed (optimum value 0).

There is often particular interest in POD and FAR for large X corresponding to major events. If X is chosen too large, however, estimates of POD and FAR may be corrupted by the low number of cases (major events are also rare events). In this study, POD and FAR were analysed for a whole range of thresholds (see Fig. 7), but results based on fewer than five cases were omitted from the presentation.

Figure 4. Reduction of variance (equation (5)(5) $\mathrm{R}\mathrm{V}=\frac{\mathrm{M}\mathrm{S}{\mathrm{E}}_{\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}}-\mathrm{M}\mathrm{S}{\mathrm{E}}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}}{-\mathrm{M}\mathrm{S}{\mathrm{E}}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}}$(5) ) as a function of lead time for different forecast variants. GG and RR denote model-based forecasts driven by raingauge or radar data, respectively. Active input/state updating is indicated by the final character “u” while “p” stands for plain forecasts without updating (error prediction was still enabled). A simplistic forecast obtained by linear extrapolation of the observed hydrograph is labelled “Extra”. The vertical line at +3 h marks the basin’s typical response time (RT)

Figure 5. Correspondence between observed (obs) and predicted (prd) flow rates for the “best” forecast according to Fig. 4 (variant “RRu”). Results are shown for lead times of 1, 2 and 3 hours (from left to right). E represents the coefficient of efficiency (Nash-Sutcliffe index)

Figure 6. Similar to Fig. 5 but displaying the increase/decrease of the flow rate (∆Q) with respect to the base value at T0. Positive values indicate rising flow, negative values mark recessions

Figure 7. POD and FAR score for predicted changes of the flow rate (∆Q) with respect to the base value at forecast time (T0). The upper limit of the x-axis was chosen so as to guarantee a minimum number of five cases behind all sections of the graphs

4 Results

All model-based forecasts were found to be skilful for lead times up to 3 h (and beyond) as indicated by a positive reduction of variance (Fig. 4). Forecasts based on radar QPE outperformed forecasts relying on interpolated raingauge data only. Interestingly, the activation of input/state updating did not lead to a consistent improvement in performance. While a slight increase in RV was obtained for the radar-based forecast, updating caused a significant drop in RV for the gauge-based forecasts (see discussion in Section 5). In the longer run, i.e. for lead times greater than the basin’s response time, the difference between the forecasts in terms of RV basically vanishes. It is emphasized that results for lead times beyond +3 h are increasingly affected by the simplistic assumptions regarding future rainfall (recall Section 3.4.2).

For benchmarking purposes, Fig. 4 also contains a graph corresponding to a purely data-based forecast (labelled “Extra”). It was generated by linear extrapolation of the observed hydrograph using just the two most recent flow records. Apparently, this simple approach also yields skilful 1 h forecasts but RV is close to zero or negative for all other lead times.

Figure 4 allows for a convenient comparison of different forecast variants. The practical usefulness of a particular variant, however, can only be assessed by means of additional analyses. In the following, the focus is on the best-performing forecast according to RV, i.e. the model-based forecast using radar QPE and input/state updating (labelled “RRu” in Fig. 4).

Figure 5 illustrates the correspondence between observed and predicted streamflow. For the 1 h forecast, all points lie close to the 1:1 line. Even though the scatter becomes more pronounced at lead times of +2 and +3 h, the results look quite impressive at first glance (E > 0.95). It must be recognized, however, that Fig. 3 does not differentiate between the contributions of the known initial conditions (i.e. streamflow at T0) and the model-based prediction. Therefore, it is informative to look at a persistence-adjusted variant of the original scatter plot (Fig. 5). Instead of absolute values (Q), this figure considers the change in streamflow (∆Q) based on the known initial value at T0. Compared to Fig. 6, the scatter plots in Fig. 4 appear much more dispersed (but note the different scaling of axes). Dispersion is particularly large in the case of positive values of ∆Q, i.e. when flow rates are increasing.

Probability of detection and false alarm rate were also computed with respect to changes, ∆Q, rather than absolute flow rates, Q. According to Fig. 6, POD hardly exceeds 0.5 and, for higher thresholds and lead times, values go down to a level of about 0.3. Thus, the fraction of cases where an observed rise in streamflow was predicted successfully (in terms of threshold exceedence) amounts to 30–50% only. At the same time, a considerable number of the predictions (20–40%) are classified as false alarms.

5 Discussion

The comparison of different forecast variants with respect to the reduction of variance (Fig. 4) indicated clear superiority of the radar-based QPE over the raingauge-based QPE. While this meets our expectations, it is not a trivial finding because other case studies (e.g. Kneis and Heistermann 2009, Gourley et al. 2011) came up with different results. In this specific case, superiority of the radar-based QPE is probably due to the unfavourable configuration of the current raingauge network (Fig. 1). The northeastern part of the basin is essentially ungauged in terms of rainfall. Moreover, most of the gauges are located at low elevations. Thus, rainfall rates in the mountains of the Sierra Madre, which are probably controlled by orographic effects, cannot be well represented. Given the current configuration of the raingauge network, the radar-based QPE should preferably be used for basin-scale hydrological applications. This finding is particularly relevant because, in spite of the unfavourable configuration, MRB is still one of the most densely gauged basins in the entire Philippines.

In the presented case, data assimilation (DA; Section 3.4.3) essentially contributed to forecast quality. Without any DA, i.e. when recent observations of streamflow were ignored, the correspondence between observed and predicted flow rates would deteriorate dramatically (results not shown). The simple error prediction scheme (EP) was effective irrespective of the source of rainfall data. The additional benefit from input and state updating (Section 3.4.3), however, was comparatively small, as indicated by the graphs “RRu” and “RRp” in Fig. 4. This finding also holds for relaxed constraints on the optimized multipliers. In fact, the overall verification statistics deteriorated if stronger corrections of rainfall data and initial storage were allowed. This was due to massive overestimation of catchment wetness in a small number of cases, possibly triggered by non-representative observations of flow, rainfall, or both. The deterioration of forecasts through updating in the case of the gauge-based QPE (graphs “GGu” and “GGp” in Fig. 4) is surprising. This gives rise to the assumption that input updating, i.e. the adjustment of precipitation estimates, may be counterproductive if the original estimate shows major deficits, e.g. with respect to spatial patterns. It appears that the drop in RV is explained by a small number of severe overpredictions while, on average, updating had the intended effect.

The closer analysis of forecast performance revealed substantial deficits with respect to the system’s predictive capabilities. These deficits become evident if the initial conditions (known flow rates at forecast time T0) are considered explicitly in the analysis of forecast errors. While the predicted flow rates correspond rather well with the observed ones (E > 0.95, Fig. 5), the correlation between forecast and observed changes is moderate at best (Fig. 6). The scatter is particularly large in the case of positive ∆Q, i.e. for predictions of rising flow rates. The tendency for underestimation is visible in both Figs 5 and 6.

The mentioned deficits are also reflected in the statistics of POD and FAR with respect to predicted flow changes (Fig. 7). Given the high number of missed events (POD ≤ 50% for positive ∆Q at lead times greater +1 h), the suitability of the forecast system for flood-warning purposes seems questionable. Besides the low POD, a considerable number of false alarms (FAR > 20%) was identified. The problem of overprediction is, however, restricted to small positive ∆Q (cf. Fig. 6).

In view of the identified deficits, several options for improvement of the forecast system need to be considered. These options may be assigned to one of three categories: (1) reduction of errors in deterministic forecasts, (2) quantification of predictive uncertainty through probabilistic forecasts, and (3) general extensions of the forecast system, as discussed in the subsequent paragraphs.

5.1 Reduction of errors

One strategy to minimize forecast errors may be targeted at enhanced data assimilation. At present, the greatest potential is seen in the exploitation of additional data sources such as gauge readings from upstream sites. If the forecast system were to be used operationally, quick and reliable data transmission would be important prerequisites. Moreover, a frequent update of rating curves would be necessary since the assimilation of incorrect observation data may result in predictions being worse than those produced with the model alone.

The optimization procedure used in the context of data assimilation (Section 3.4.3) is a rather simplistic one and alternative algorithms should be explored in future studies. Algorithms eligible for operational applications ideally share the following features: (1) the balance between computational speed and accuracy should be tunable via parameters; (2) support for box constraints is important unless a confrontation with corrupted observation data is definitely impossible; (3) non-finite return values of the objective function, e.g. due to a crashed model, need to be handled; and (4) a deterministic algorithm should be selected if reproducibility of results matters as, e.g., in verification experiments. In the context of this study, the employed line-search procedure proved to be robust and computational speed was no cause for concern. For more demanding applications, requiring data assimilation at multiple gauges, for example, the algorithm could be adapted to allow for parallel processing.

5.2 Quantification of uncertainty

Although a further reduction of the error in deterministic forecasts may be achievable, predictive uncertainty cannot be fully eliminated. Consequently, quantification of uncertainty needs to be addressed whenever the forecast system is to be used in decision making, e.g. flood warning. Many concepts to produce probabilistic forecasts by means of deterministic hydrological modelling have evolved over the past two decades. Usually, different sources of uncertainty need to be distinguished, such as inputs, model structure and parameters (e.g. Montanari and Brath 2004) as well as initial conditions. While input uncertainty is typically accounted for explicitly by means of meteorological ensembles, the contribution of all other sources is often treated in a lumped way, commonly termed “hydrological uncertainty” (Krzysztofowicz 1999, 2001). The approaches used to quantify hydrological uncertainty largely fall into one of two groups. The first group encompasses simulation techniques such as GLUE (Beven and Binley 1992), and more formal Bayesian methods (Engeland and Gottschalk 2002, Liu et al. 2005, Sadegh and Vrugt 2013). The second group consists of statistical methods to predict the error in forecasts produced by the mechanistic hydrological model (Krzysztofowicz 1999, Seo et al. 2006, Todini 2008, 2013, Weerts et al. 2011, Brown and Seo 2013). Statistical post-processing can be applied to both single-valued forecasts or ensembles originating from explicit treatment of input uncertainty. The statistical methods are particularly attractive for operational uses (see, e.g. Demargne et al. 2014) as they can treat uncertainty from various sources at low computational cost.

In this case study, a quite limited number of hindcasts was examined, making it impossible to calibrate an advanced statistical error prediction model. In order to quantify the total uncertainty anyway, a simple but straightforward approach was used to generate a probabilistic forecast of the change in flow rates, ∆Q, (cf. Fig. 6). Similarly to Weerts et al. (2011), an additive model was adopted (equation (6(6) $\mathrm{\Delta }Q=\mathrm{\Delta }{Q}_{\mathrm{p}\mathrm{r}\mathrm{d}}+e$(6) )). Here, ∆Q _prd represents the deterministic prediction and e is the corresponding error:

(6) $\mathrm{\Delta }Q=\mathrm{\Delta }{Q}_{\mathrm{p}\mathrm{r}\mathrm{d}}+e$(6)

Quantile-quantile plots suggest that the transformed variable $\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(e)\cdot \mathrm{l}\mathrm{o}\mathrm{g}(|e|)$ follows a normal distribution reasonably well if cases where $|e|<1$ (not being of practical concern anyway) are dropped (Fig. 8). Hence, confidence limits of ∆Q can be obtained by picking a quantile from that normal distribution at the desired probability, e.g. at p = 0.95, and adding/subtracting the back-transformed quantile to/from the deterministic forecast ∆Q _prd.

Figure 8. Normal quantile-quantile plots of transformed forecast errors for lead times up to +3 hours

Selected values of e for a number of probabilities are reported in Table 3. These are “global” estimates of forecast uncertainty conditioned on the studied events. Attempts to produce more case-specific error estimates by (1) relating e to the value of the deterministic forecast, or (2) exploiting information on the flow rate at T0 or the last known model error were not successful. Correlations between e and these potential predictors were generally weak. Table 3 differentiates between forecasts of rising and falling flow rates at least.

Table 3. Statistical errors of the deterministic forecasts ∆Q _prd for different lead times, probability levels and signs. All values correspond to model set-up “RRu” (best-performing variant according to Fig. 4). Example: Row 2 indicates that predictions of increasing discharge at a lead time of +1 hour are associated with absolute errors of less than about 15 m³ s⁻¹ with a probability of 90%

Lead time (h)	Probability (%)	Forecast error e (m^3 s^−1)
		if ∆Q prd ≥ 0	if ∆Q prd < 0
1	50	±3.0	±2.4
	90	±14.6	±8.1
	99	±66.5	±26.5
2	50	±3.6	±2.7
	90	±23.0	±11.3
	99	±135.2	±44.5
3	50	±3.9	±2.9
	90	±28.0	±14.1
	99	±184.9	±63.3

5.3 Extensions of the forecast system

Our study focused on lead times of up to 3 hours. Within this time frame, runoff predictions are not sensitive to future precipitation. This was also verified experimentally by re-running the verification experiment with alternative assumptions on future rainfall (e.g. zero values). For practical use, e.g. in flood early warning, an extension of the forecast range beyond the basin’s response time may be desirable. Based on these premises, the generation of radar-based, short-term precipitation forecasts is currently under development. Nevertheless, given the deficits in the short-term forecasts (Figs 6 and 7), the possible quality of runoff predictions for lead times greater than +3 h is limited.

There are many options to further develop the hydrological model engine and to refine parameter estimates. For example, one could explore alternative sub-models for runoff generation, distinguish a greater number of land uses, or increase the overall spatial resolution. From the perspective of an operational user, however, it might be more interesting to see whether an empirical model, e.g. based on an artificial neural network, was a competitive substitute for the conceptual hydrological model. Empirical alternatives to a conceptual model are generally attractive for their unbeaten computational performance and low data requirements. The training of empirical models employed in forecasting, however, requires special techniques (Schmitz and Cullmann 2008) to compensate for lack of extrapolation capability.

6 Conclusions

The verification of model-based runoff predictions over two seasons revealed severe deficits in forecast quality. This is true even for short-term predictions entirely controlled by observable rainfall rather than by future precipitation. The performance is particularly limited with regard to the prediction of strong rises in streamflow. The observed deficits are probably rooted in both the insufficient quality of quantitative precipitation estimates (QPE) and the simulation of catchment-scale hydrological processes. Regarding QPE, the case study yields a clear recommendation: given the present configuration of raingauges, runoff forecasting for the study basin should rely on the radar-based QPE.

Storm runoff from the Marikina River Basin has repeatedly caused casualties and high economic losses in recent years. Reliable early warnings are, therefore, of vital interest. Given that, further investment in the forecast system is necessary to lessen the identified deficits. Priorities include the use of additional gauge readings for model updating and further improvements of radar-based QPE, e.g. by enhanced hardware calibration, by applying dynamic raingauge adjustment, and by combining observations from different radar stations in regions of overlap. (See Crisologo et al. (2014) for progress in radar-based QPE for a C-band radar that also covers the Marikina Basin.)

Although the implementation of those measures will reduce the errors in deterministic forecasts, considerable uncertainty is likely to remain. In this context, the parametric uncertainty post-processor outlined in Section 5 should be compared against alternative statistical approaches such as quantile regression. Also, the parameters of any statistical post-processor need to be verified/updated when new hindcasts become available. Finally, more effort must be put into the generation of case-specific estimates of predictive uncertainty from all operationally available information, e.g. by means of generalized linear models.

Extension of the hydrological forecast range via rainfall nowcasting will become a priority task as soon as significant progress has been made at short lead times.

Acknowledgements

All meteorological and hydrological observation data were provided by EFCOS (Effective Flood Control Operation and Services). Special thanks go to the Environment Monitoring Laboratory (EML) members for providing equipment and assistance for field campaigns.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The study was financially supported by the German Federal Ministry of Education and Research (PROGRESS project, work group D2.2, grant no. 03IS2191A). C. Abon was funded through the GeoSim graduate school.