## Abstract

This paper describes a fully automated scheme that has provided calibrated 1–10-day ensemble river discharge forecasts and predictions of severe flooding of the Brahmaputra and Ganges Rivers as they flow into Bangladesh; it has been operational since 2003. The Bangladesh forecasting problem poses unique challenges because of the frequent life-threatening flooding of the country and because of the absence of upstream flow data from India means that the Ganges and Brahmaputra basins must be treated as if they are ungauged. The meteorological–hydrological forecast model is a hydrologic multimodel initialized by NASA and NOAA precipitation products, whose states and fluxes are forecasted forward using calibrated European Centre for Medium-Range Weather Forecasts ensemble prediction system products, and conditionally postprocessed to produce calibrated probabilistic forecasts of river discharge at the entrance points of the Ganges and Brahmaputra into Bangladesh. Forecasts with 1–10-day horizons are presented for the summers of 2003–07. Objective verification shows that the forecast system significantly outperforms both a climatological and persistence forecast at all lead times. All severe flooding events were operationally forecast with significant probability at the 10-day horizon, including the extensive flooding of the Brahmaputra in 2004 and 2007, with the latter providing advanced lead-time warnings for the evacuation of vulnerable residents.

## 1. Introduction

Accurate and timely forecasts of river flow have the potential of providing critical information for water resource management, agriculture practice optimization, and disaster mitigation. Nowhere is the need for such forecasts more urgent than in the Bangladesh delta that lies at the confluence of two of the largest river systems in the world: the Brahmaputra and the Ganges (Fig. 1). Combined, these rivers can possess discharges of well over 100 000 m^{3} s^{−1}. Once these flood-stage waters reach Bangladesh—a small country the size of England (140 000 km^{2}) but with a population of 155 million—they spread out onto a flat and wide floodplain occupying roughly 80% of the country, making Bangladesh part of the world’s largest river delta. The combined catchment area of the Brahmaputra–Ganges system (approximately 1.5 × 10^{6} km^{2}) ranks tenth in the world, with only the Amazon and the Congo surpassing the combined discharge of the two rivers.

Each year, short-lived flooding occurs throughout the summer and early autumn with sufficient irregularity to produce adverse social and agricultural effects. In these normal years, about one-fifth of the country can be inundated for days. The timing of the flooding is critical for the harvesting of the winter rice crop, the planting of summer rice, and the viability of other crops planted during the monsoons. Major flooding events occur in Bangladesh with a return period of 4–5 yr. For example, during the summer of 1998, more than 60% of Bangladesh was inundated for nearly three months (Ahmad and Mirza 2000). Again in July 2004, widespread flooding occurred as the Brahmaputra exceeded dangerous flood levels twice over a two-week period. In July and September 2007, the Brahmaputra severely flooded twice with inundation lasting weeks.

There is considerable utility in producing timely and accurate river flood forecasts for Bangladesh (RTI and EGIS 2000). Such forecasts can impact short-term agricultural and societal decision making, specifically when to evacuate flood-prone regions and ensure protection of food and water supplies and livestock safety. Preemptive actions taken can significantly lower relocation and reconstruction costs following a period of flooding. RTI and EGIS estimates that a 7-day forecast has the potential of reducing postflood costs by as much as 20% in household and agricultural costs compared to 3% with a 2-day forecast.

Predicting river discharge in Bangladesh is little different from that existing in many developing nations but with one very particular problem: the Brahmaputra and Ganges Rivers exit India to flow into Bangladesh, and there has been little consistent river discharge data sharing between these two countries. The only reliable river streamflow data come from what Bangladesh measures once the rivers cross the India–Bangladesh border (see Fig. 1). As a result, Bangladesh forecast lead times have been limited to 2–3 days for the interior of the country and have essentially no lead time in areas close to the border.

In the absence of upstream river discharge data from India, a different modeling approach would need to be taken to extend warning times of severe flooding. The Ganges and Brahmaputra catchments upstream of the India–Bangladesh border would have to be treated as if they were ungauged and effectively naturalized rivers. Can the latter be an accurate approximation for these highly utilized river basins?

The Brahmaputra has yet to have a major hydraulic structure built along its course (Singh et al. 2004) and as such, it can be modeled as a naturalized river. Although there are several irrigation projects along its course, their design discharges are typically small (well below 100 m^{3} s^{−1}). However, the management of the Ganges River presents a further complication, primarily as a result of at least 34 barrages and other structures functional in India and Nepal (Mirza 2004). The actual amount of the water diverted by these structures is not known. However, these structures were designed primarily for use in the dry season and not during the monsoonal flood season, which is the period of interest of this study. For example, the Farakka Barrage, a major diversion located 10 km from the India–Bangladesh border, is diverted primarily to reduce silt in the Calcutta harbor during the dry season. During the monsoon season (June–October), there appears to be little difference between discharge–precipitation relationships before or after the 1974/75 construction of the barrage (Jian et al. 2009). On the basis of the Jian et al. (2009) study and the Flood Forecasting and Warning Centre (FFWC; 2000, personal communication), it is assumed Farakka diversions do not affect discharge into Bangladesh beyond 15 June.

The catastrophic Bangladesh flooding in 1998 encouraged the U.S. Agency for International Development (USAID) to fund an exploratory project called the Climate Forecast Applications in Bangladesh (CFAB) in 2000. The primary goal of CFAB was to provide advanced warning of Brahmaputra and Ganges flooding and drought in Bangladesh on time scales ranging from daily to seasonal, in partnership with the FFWC and the European Centre for Medium-Range Weather Forecasts (ECMWF). By utilizing ECMWF ensemble weather forecasts and conditional estimates of hydrologic modeling error, the CFAB discharge forecasts were designed to also provide forecasts of their own uncertainty. Operational (probabilistic) forecast dissemination to Bangladesh began during the monsoon season of 2003, based on ECMWF’s medium- and seasonal-range ensemble forecasts. Table 1 describes the CFAB three-tier system and illustrates the concept of overlapping forecasts to allow a nesting of strategic decisions and subsequent tactical adjustments. In this paper, we focus on the medium-term forecasting tier. (The project is described at http://cfab.eas.gatech.edu, and the three-tier forecasts generated as per Table 1 presented at http://cfab2.eas.gatech.edu.)

Section 2 provides a description of the medium-range forecast scheme. Five components are introduced and described in detail. Examples of the forecasts for 2003, 2004, and 2007^{1} for the Brahmaputra and the Ganges are presented in section 3 along with forecast verifications. Section 4 summarizes the utility of the model and why it is necessary to include meteorological forecasts in a hydrological scheme. Possible extensions of the scheme are also described in section 4.

## 2. The hydrological prediction scheme

Figure 2 provides a flowchart of the Bangladesh multimodel discharge forecasting procedure. Step I describes the initial data inputs from ECMWF, National Aeronautics and Space Administration (NASA) and National Oceanic and Atmospheric Administration (NOAA) satellite and rain gauge products, river discharge measured at the India–Bangladesh border by the FFWC, and updated parameters for the hydrological models. In step II, the ECMWF ensemble predictions of precipitation are statistically modified using the satellite and rain gauge data to remove systematic errors. Together with the discharge data, these data are ingested into a set of hydrological models in step III. Step IV takes into account forecast uncertainties while preserving the useful information of the ensemble weather forecast dispersion, and it makes a final error correction to produce the final probabilistic forecast of river discharge. Step V is the dissemination of the forecasts, which is not discussed in detail in this paper. We note that the complete forecasting procedure shown in Fig. 2 is fully automated.

### a. Step I: Initial data

There are four major sets of initial data in Fig. 2: the ensemble forecasts from ECMWF (box E), satellite and rain gauge data (box S), river discharge data (box Q), and hydrological model parameters (box MP). The arrow at the base of the boxes shows where the product or derivative of the product is used.

#### 1) ECMWF (EPS) forecasts (box E)

A catchment integrates spatially and temporally all weather factors (precipitation, evaporation, etc.), together with hydrological effects (such as aquifer flow and soil moisture changes), into river flow and water storage and eventually river discharge out of the basin. Because of the spatial scale of the Ganges and Brahmaputra basins, the relatively large-grid global circulation weather model output from the ECMWF ensemble prediction scheme (EPS) can be used to generate useful discharge forecasts for these essentially ungauged basins. The ECMWF weather variables used in this study are surface fields of wind, humidity, and precipitation. The weather forecasts are provided as 51 ensemble members for each variable and for each forecast lead time. The model resolution was initially at approximately 80 km × 80 km grid (TL255L62) going out to 240 hours at 6-h increments; however, the current operational EPS is at approximately 50 km × 50 km grid (TL399L62 from 0 to 10 days), with the forecasts horizon also extending to 15 days (days 9–15 at TL255L62; Buizza et al. 2007). All forecast fields are interpolated to the same ½° × ½° grid used by the semidistributed hydrological model (SDM).

The joint CFAB–ECMWF effort was one of the first efforts to utilize ensemble medium-range numerical weather forecast input in the generation of operational short-to-medium-range ensemble discharge forecasts. The ECMWF system provides not only weather variable forecast information but also estimates of the likely forecast error through both 51-member medium-range ensembles (Molteni et al. 1996). The ECMWF ensemble predictions are passed in step II for statistical evaluation and correction.

#### 2) Satellite and rain gauge precipitation estimates (box S)

For forecasting discharge of large spatial-scale catchments, and therefore long hydrological temporal scales, it is essential that the initial state of the hydrologic system (soil moisture content, groundwater storage, in-stream flows, etc.) be assessed as accurately as possible. The importance of setting these initial conditions correctly can be understood through the following simplified argument: *i*-day lead-time discharge forecasts require *i*-day weather forecasts and (*n* − *i*) days of weather observations, for watersheds of time scale *n*, where *n* can be appreciable for large river basins.

Three precipitation observations were integrated in the operational forecast system in 2004. One source is ground-based ½° × ½° gridded daily averaged rainfall estimates from the NOAA Climate Prediction Center (CPC) utilizing World Meteorological Organization (WMO) global telecommunication system (GTS) reporting rain gauges and using Shepard’s inverse-distance algorithm for interpolation (Xie et al. 1996). Note, however, that the reporting gauge density is sparse over the Brahmaputra and Ganges catchments, with roughly 40 gauges reporting each day. The other two estimates are near-real-time ¼° × ¼° 3-hourly satellite-derived precipitation accumulations from the NASA Tropical Rainfall Measuring Mission (TRMM; Huffman et al. 2005, 2007) and the NOAA CPC morphing technique (CMORPH; Joyce et al. 2004) programs, both of which utilize infrared- and microwave-frequency satellite data but distinctly different algorithms.

These datasets are continuously updated as they become available and are interpolated to the same hydrologic computational temporal and spatial grid, with a simple grid-by-grid arithmetic average of the available data taken as the best estimate of the observed precipitation field at the time of daily forecast initialization. This observational estimate is used in both the catchment lumped model (CLM) and in the initialization (soil moisture states and in-stream flows) of the SDM. All other observed weather variables used in SDM evapotranspiration calculations are derived from the ECMWF analysis fields.

#### 3) River discharge data (box Q)

An essential data source for skillful forecast generation are discharge and stage measurements at the near-border gauging stations of Bahadurabad for the Brahmaputra and Hardinge Bridge for the Ganges (Fig. 1b) provided by the FFWC. These stations have observed discharge records derived from standard river cross-section survey techniques completed approximately once per week for more than 13 years and daily river stage records of more than 52 years. These historical discharge records are essential for model calibration, and the near-real-time observations are utilized for hydrologic model data assimilation and error corrections considered later in step III of Fig. 2 and for the final discharge forecast probability distribution function (PDF) generation considered in step IV of Fig. 2.

Daily discharge estimates were derived from fitting rating curves to the FFWC observational data using rating curve functions of form

where *Q _{t}* is the estimated discharge for day

*t*,

*H*is the observed river stage (relative to a fixed datum), and

_{t}*i*represents the stage interval for which the equation applies. The Ganges used a two-part rating curve

*i*∈ [1, 2], and the Brahmaputra, a three-part form. Because sediment transport and bed mobilities of both the Brahmaputra and Ganges are substantial, the parameter sets of these equations were fitted independently for each year over the historical record.

^{2}This was done by weighting the residual error inversely proportional to the year separation and minimizing it using a nested version of the downhill simplex method in multidimensions [Nelder and Mead 1965; A Multidirectional Optimum Ecotope-Based Algorithm (AMOEBA) procedure; Press et al. 1992].

Absolute rating curve error ɛ is estimated to be |ɛ| = 0.10*Q* for Bahadurabad and |ɛ| = 0.07*Q* for Hardinge Bridge. By accounting for additional systematic measurement errors in the FFWC method to measure the observed river discharges of order |ɛ| − 0.1*Q* (rough) and by assuming independence of these sources of error to combine them in quadrature (Taylor 1982), we estimate total error bounds on the observed discharges (shown in Figs. 6, 8, 9, and 10) of |ɛ| ≈ ±0.15*Q*.

We note that both of these rivers have low gradients around the forecasting locations, so it was anticipated that pressure gradient effects due to variable surface water slopes, along with effects due to the filling and draining of overbank storage, may limit the utility of fixed stage–discharge relationships via this potential hysteresis. However, an attempt to minimize the 10% rating curve error just discussed for the Brahmaputra by including an additional hysteresis correction term in Eq. (1) never affected the estimated discharge more than 0.2%.

#### 4) Updated SDM parameters (box MP)

The SDM model parameters are continuously updated through an automated background calibration process that uses a modified version of the root-finding AMOEBA procedure (Nelder and Mead 1965; Press et al. 1992) to minimize the square error between observed and modeled discharge. The search is performed over a reduced nine-parameter set to maintain a small ratio of adjustable parameters to daily discharge calibration points, and thus reduce the risk of overcalibration. These nine parameters are discussed in the SDM model, section 2c. If, at the time of daily discharge forecast initialization, a new optimal root has been found, soil moisture states and in-stream flows are automatically updated before the forecast proceeds.

### b. Step II: Statistical rendering and downscaling of forecasts (dotted box, Fig. 2)

Generally, all practical hydrologic models require some form of calibration arising from reasons such as incomplete knowledge of watershed properties and the necessary parameterization of transport at unresolved scales. Through the hydrologic model parameter calibration process, some of the additional errors in the watershed inputs (rainfall, in particular) can also be implicitly reduced (e.g., runoff errors arising from an overbias in rainfall forcing can be implicitly minimized by increasing the parameterization of evapotranspiration). However, such bias reduction through calibration does not occur if there are relative errors between more than one input of the same physical quantity. In particular, using both weather forecasts and observationally based estimates of precipitation concurrently to generate unbiased discharge forecasts requires that these two data sources maintain their statistical similarity, such that there are no relative biases between their statistical moments (mean, variance, and skewness, in particular).

To minimize systematic differences in hydrologic model inputs leading to systematic hydrologic forecasting errors, a quantile-to-quantile (q-to-q) mapping technique was implemented. A similar approach has been applied to adjust nonsimultaneous radar data to rain gauge observations (Calheiros and Zawadzki 1987). As applied here, the q-to-q technique forces each ECMWF ensemble precipitation forecast at each ½° × ½° grid point (input E, Fig. 2) to be statistically sampled from the cumulative distribution function (CDF) of the associated observational historical record for the same grid point (input S, Fig. 2). In essence, this converts the numerical weather prediction product from a physical variable to a probability forecast by performing a bias correction at each quantile of the forecast CDF. This bias correction approach was implemented for the 2004 monsoon season for catchment lumped values of the ECMWF 40-member seasonal forecasts (Webster et al. 2006) and for each forecast grid of the ECMWF medium-range forecasts (Hopson 2005) the same year. The manner in which the technique is applied to the medium-range forecasts is described below.

Calculate the monsoon season “climatological” empirical CDFs for both the observations and the forecast model at each ½° × ½° grid location *i*. The observational (subscript *o*) climatological CDF for *i* and precipitation amount *y*, *C _{o}^{i}*(

*y*), is derived from the previously discussed rain gauge (available from 1978 onward) and merged gauge–satellite data (available from 2002 onward). The forecast (subscript

*M*) CDF’s at

*f*and

*i*,

*C*

_{M}^{i,f}(

*y*) are derived from stored forecasts (available from 2003 onward), which are first downscaled through simple two-dimensional linear interpolation to the ½° × ½° computational grid of the rainfall observations and hydrologic models, and the CDF computation is performed independently for each forecast

*f*(24, 48 h, etc.). Formally, these CDF’s are given by

where *P*{*α*} denotes the probability of the event *α* and *x* (*y*) is a specific value of the generic random variable *X* (*Y*) representing the observed (forecast) weather variable of interest. Using these CDFs, a direct quantile-by-quantile bias correction is made to each grid point *i*, *f*, and precipitation ensemble member fcst of the forecast, , of the forecast by simply setting (1) equal to (2) and determining the value of *x* = *p*_{adj} that satisfies the equation. In this manner, the “observation space” quantile is matched with the “forecast model space” quantile as shown in Fig. 3, and *p*_{adj} is used in replace of in the hydrometeorological application.

Note that this technique ensures that the forecasts produce the same climatological number of “no rain” events as observed. It also preserves the spatial and temporal (ranked) covariances of the weather variable forecast fields (as generated by the numerical weather prediction model). Figure 4 shows an example of both the spatial- and temporal-preserving ability of this approach. The top panel is a 144-h accumulated field of the original ECMWF 24-h forecasts from 20 to 25 July 2007 preceding the severe Brahmaputra flooding event from 26 July to 6 August. The lower panel shows the same accumulated field after being “q-to-q mapped” to the observational PDF. Note that structures in the original forecast precipitation field have been preserved.

If there is no contemporaneous forecast and observational data, we suggest that the q-to-q technique may be the best one for bringing forecasts into statistical similarity with the observations. However, although the mapping improves the forecast bias at each quantile by steering the forecasts into having the same frequency of events as the observation climatological record, this process still does not necessarily ensure that the forecasts’ skill will be increased (which we discuss elsewhere in detail). Such skill degradation can (arbitrarily) arise because the technique does not directly improve the forecast model’s reliability by neglecting the conditional relationship between observations and forecasts. Although not discussed here, to improve the forecast’s conditional bias and ensemble dispersion, following the gridded q-to-q bias correction, we employ quantile regression to the catchment-averaged values of the precipitation forecasts.

### c. Step III: Hydrological modeling (hashed boxes, Fig. 2)

Next we discuss the two hydrological model approaches we employ along with their relative attributes and drawbacks: the CLM and the SDM. These are combined into a multimodel following the philosophy of Krishnamurti et al. (1999), Krishnamurti et al. (2000), and Georgakakos et al. (2004), who noted collectively that some models are better equipped in certain forecast circumstances than others and that the combined performance of the multimodel generally outperforms any one model. The aim of the multimodel approach applied here is to minimize hydrological modeling error, with total hydrologic forecast uncertainty dealt with separately in step IV. The availability of near-real-time discharge observations allows for autoregressive (AR) error corrections to be applied to the SDM forecasts before the two hydrological models are combined into the multimodel prediction, which is also discussed. A flowchart of each of the hydrological models and the manner in which they are combined is shown in Fig. 2.

#### 1) CLM

The lumped modeling technique employed here is based on the work of Young and Beven (1994) and Beven (2000), whose methodology has been applied successfully in a flood forecasting scheme for the CI6 catchment at Llyn Briane in Wales and for the catchment of the River Nith upstream of the town of Dumfries in Scotland (Lees et al. 1994). This approach formed the basis of the 2003 CFAB operational ensemble forecasts for the Ganges and Brahmaputra.

The CLM employs a multiple linear regression between observed discharge and observed and forecasted precipitation, together with a nonlinear “effective rainfall” filter. The model is based on the idea of linear storage reservoirs, and the model structure is similar to the unit hydrograph approach (Sherman 1932), the transfer function plus noise (TFN) models of Box and Jenkins (1970), and to an autoregressive moving average model with exogenous variables (ARMAX) time series forecasting model (Haltiner and Salas 1988; Marco et al. 1993), where here the “exogenous variable” is a filtered (simultaneous and forecasted) catchment-averaged rainfall time series but the “moving average” process is not being explicitly modeled. The benefits of the CLM are that it has a flexible structure, it can be easily recalibrated daily for different forecast time horizons and for different in situ conditions, and it easily assimilates near-real-time discharge measurements. Drawbacks of the CLM include limited long-time-scale base-flow modeling and heterogeneities in soil moisture and the distribution of rainfall–runoff across the watershed for which they are not accounted. However, the CLM approach can still be quite effective under certain settings, such as when precipitation events are of a similar spatial scale to the catchment itself, or less strictly, when there is a sufficiently large number of precipitation events compared to the catchment’s inherent time scale, and their spatial distribution tends to maintain statistical similarity over time.

The CLM approach, employed operationally, uses the time series of both observed satellite–rain gauge combined product, as discussed earlier, and the daily ECMWF ensemble forecasts of catchment-averaged rainfall. These are adjusted to an estimated fraction of effective rainfall contributing to storm runoff *u _{i}* using the nonlinear filter

where *Q _{i}* is the observed discharge at time

*i*(or taken as the most recently observed discharge in the case of a precipitation forecast),

*Q*

_{max}is an arbitrary normalization constant,

*c*is a power to be calibrated, and

*R*is the catchment-averaged rainfall (observed and forecast). This filter provides a simple way to estimate effective rainfall dependent on watershed discharge used as a surrogate for antecedent soil moisture conditions.

_{i}Using the derived effective rainfall time series of (4) and the time series of daily observed catchment discharge, the linear forecast equation for the future discharge *Q*_{t+f} is given by

where *t* is the current day, *f* is the forecast lead time in days, the sets {*a _{i}*} and {

*b*} are coefficients, and ɛ

_{j}_{t+f}is the model error. For model calibration, all

*u*are derived from observed precipitation. However, in forecasting, the set {

_{i}*u*

_{t+f},

*u*

_{t+f−1}, … ,

*u*

_{t+1}} is replaced by individual members of the 51-member ECMWF ensemble precipitation forecasts, which preserves the temporal covariability of the forecasted precipitation. This then produces a 51-member ensemble of discharge forecasts whose dispersion represents the effects on discharge uncertainty of input forecast precipitation uncertainty.

The number and value of coefficients {*a _{i}*} and {

*b*}, the power

_{j}*c*in (4), along with the length of a boxcar smoothing window performed on the time series of

*u*, are recalibrated before each daily forecast initialization using updated observed discharge data. This calibration is done independently for each forecast lead time, with the best model determined by having the minimum Bayesian information criterion (BIC; Schwarz 1978; Katz 1981). The daily calibrated models are then used to generate forecasts, as discussed earlier, and also to calculate hindcasts (box CLM, Fig. 2), which are used in the multimodel selection process (box M, Fig. 2). Note that observed precipitation is used in generating the hindcasts to isolate model uncertainty from weather variable forecast uncertainty in the multimodel selection, discussed further in step IV. Both the forecasts and hindcasts are then passed to the multimodel step.

_{i}#### 2) SDM

In 2004, an SDM approach was included into the operational system (box SDM, Fig. 2) under the multimodel framework with the CLM. In contrast to the CLM, the strengths of the SDM approach are that water-balance physics (evapotranspiration, soil storage, etc.) are more explicitly modeled locally and the time delays of spatially distributed precipitation events are taken into account. The drawbacks are that model recalibration is more costly and therefore cannot be easily tuned to the current conditions of the watershed. Similarly, data assimilation of observed discharge into the forecast structure is not as directly applied as with the CLM.

The SDM employed here is a simple nine-parameter two-layer representative area soil model parameterizing most of the physical processes of water storage and transport with a slow and a fast time-scale response. Water balance calculations are carried out on a 0.5° × 0.5° computational grid, with each grid treated as a subcatchment. At this scale, heterogeneous water storage and transport processes are homogenized, with the model providing only conceptual representations of these processes. This model was designed to be of intermediate complexity but with a limited number of parameters, because both the large spatial scale and the lack of available data in the Ganges and Brahmaputra catchments for both discharge and soil and vegetation data limits the benefits of employing more complex models. Furthermore, the use of more complicated models increases the inherent risk of model overcalibration, as was discussed previously. The coupled nonlinear differential equations for the two-layer model are

where *S*_{1} and *S*_{2} are the amounts (length) of water in the upper and lower soil layers, respectively; *r*_{s1} and *r*_{s2} are the sizes (length) of the upper and lower soil layer reservoirs, respectively; *t*_{s1} and *t*_{s2} are the time constants for lateral outflow from the upper and lower soil layers, respectively; *t _{p}* is the time constant for percolation from the upper to the lower soil layer;

*P*is the precipitation rate (length/time); and

*E*is the potential evaporation rate (length/time), which is described below. However, if

_{p}*S*

_{2}>

*r*

_{2}, then the two-layer model simplifies to

The outflow of water (length/time) is given by

The discharge calculations for each grid are calculated as a daily rate using Eq. (10) along with integrations of the two rate equations, Eqs. (6) and (7). These integrations are performed on 6-h intervals using a fourth-order Runge–Kutta method to generate the final daily discharge value.

The potential evaporation *E _{p}* is calculated using the aerodynamic method (Thornthwaite and Holzman 1939) using the ECMWF analysis and control (deterministic member) forecast fields and with an additional calibration coefficient. Although the aerodynamic method is derived for conditions where moisture supply is not limited, the evaporation losses here are made to depend linearly on the moisture availability of the upper soil surface through the factor (

*S*

_{1}/

*r*

_{s1}) in (6).

To route the discharge out of each subcatchment and down to the forecasting location and to simulate a typical watershed discharge response, the subcatchment is treated as a series of identical linear reservoirs (e.g., Nash 1958; Chow et al. 1988). For our application, this is done using an *n*-reservoir instantaneous unit hydrograph with its own decay constant *t _{c}* (both

*n*and

*t*are calibrated). For the

_{c}*j*th subcatchment, this is given by

where Γ(*n*) is the gamma function, *m* is the number of days prior to *t* presumed to generate outflow (operationally taken as 40), and *T*_{lag}*d ^{j}* represents the subcatchment travel-time delay from the subcatchment outlet to our forecasting location, where

*d*is the (normalized) distance from the subcatchment outlet to our forecasting location, and

^{j}*T*

_{lag}is a calibration coefficient (prescribed over the whole watershed). Summing over all subcatchments,

*j*then gives the total discharge at the forecasting location.

Taken all together, there are five parameters of the two-layer model (*r*_{s1}, *r*_{s2}, *t*_{s1}, *t*_{s2}, *t _{p}*), one evapotranspiration parameter, and the three routing parameters (

*n*,

*t*, and

_{c}*T*

_{lag}), comprising nine calibration parameters for the whole watershed. The calibration of these parameters is discussed in section 2a(4). For each subcatchment, there are the two state variables

*S*

_{1}and

*S*

_{2}and four forcing variables (precipitation, wind, humidity, and temperature).

Note that snowmelt is not explicitly modeled in the SDM as presented here. Although this is a model deficiency, monsoon rainfall is by far the dominant source of river discharge during the summer months. For the Ganges, snow meltwater provides much less than 1% of the annual river flow volume (Central Water Commission 1988). However, for the Brahmaputra, roughly 10% of the monsoon runoff is derived from snowmelt. If one estimates natural variability in snowpack on the order of 10%, then only a few percent of the seasonal discharge variability is neglected without including snow runoff processes.

For each forecast initialization, discharge hindcasts are generated using the updated model parameters of box MP in Fig. 2 and using the ECMWF weather variable analysis fields. Precipitation, however, comes from the hourly updated satellite and rain gauge observations discussed earlier (box S, Fig. 2). Forecasts are then calculated forward from the updated current hydrologic state using the adjusted ECMWF EPS 51-member precipitation ensemble (output from step II, Fig. 2) and the deterministic forecast for all other weather variables.

As with the CLM, ensemble precipitation forecasts are transformed into ensemble discharge forecasts via the hydrologic modeling process, and the effect of the temporal covariances of the precipitation fields is preserved in the resultant ensemble discharge forecasts by tracking single-ensemble-member precipitation forecast fields through the hydrologic model from 1 to 10 days. Both the discharge forecasts and hindcasts are then passed to the error correction step discussed in the next section, which assimilates the FFWC discharge observations at Bahadurabad and Hardinge Bridge (box Q, Fig. 2), available at the time of forecast initialization.

#### 3) Autoregressive error corrections

Young (2002) noted that the accuracy of operational hydrological forecasts could be improved if real-time corrections are applied to improve model performance utilizing recent discharge observations and hydrologic model errors. For our SDM discharge forecasts, we employ an AR data assimilation technique to forecast combined hydrologic model and weather observation error. This model is refit daily for each forecast lead time using SDM hindcasts in which estimated satellite–rain gauge precipitation fields are used instead of weather forecasts (discussed later). Using estimated precipitation fields isolates hydrologic forecasting errors due to deficiencies in hydrologic modeling and estimated precipitation from errors due to weather forecast uncertainty (which will be combined in step IV). The technique is similar to the stepwise AR technique discussed by Granger and Newbold (1986, their section 5.4) and is given by

where *R _{f}*(

*t*+

*f*) =

*Q*(

_{o}*t*+

*f*) −

*Q*(

_{f}*t*) is the hydrologic model residual (error) at time (day)

*t*+

*f*for a modeled discharge forecast

*Q*(

_{f}*t*) initialized at time

*t*with a lead time of

*f*days in advance, where

*Q*(

_{o}*t*+

*f*) is the contemporaneous observed discharge at time

*t*+

*f*,

*μ*is the residual time series mean, and ɛ

_{t+f}is the leftover error from the AR fit.

Note that for the AR model (14) to provide the best linear unbiased estimator of *R _{f}*(

*t*+

*f*), the errors ɛ

_{t+f}need to be homoscedastic and (generally) independent (see, e.g., Helsel and Hirsch 1992). However, stationarity in the errors is generally not valid (not shown), because when the modeled (and observed) discharges are larger (smaller), the variance of ɛ

_{t+f}is generally also larger (smaller). Although a Box–Cox transform could be applied in this case to remove such heteroscedastic behavior, the utility of this is not always effective (see, e.g., Chatfield 1996). And more importantly, because (14) is updated and automatically applied to the operational forecasts each day, our motivation is only to make an initial improvement to the SDM, and then in a later step, we utilize a nonparametric nearest-neighbor approach to account for any remaining nonoptimal, correlative, and heteroscedastic behavior in the final hydrologic (multi)model time series forecasts, as discussed in step IV (Fig. 2). This fitting of (14) is carried out using stepwise regression using forward selection based on the minimization of the BIC. In addition, the selected AR model is also compared with up to fifth-order no-fit Taylor series models.

The optimal error correction model for each forecast lead time then “error corrects” the SDM hindcasts exactly as if they were actual forecasts (but driven with observed precipitation and analysis weather variable fields), as well as forecasting the error of (and then correcting) the current ensemble hydrologic forecasts

Both the error-corrected hindcasts and corrected ensemble forecasts *Q _{f}^{c}*(

*f*) are then sent to the multimodel forecasting step (box M, Fig. 2). It should be noted that because the error correction is applied to the hindcast time series driven by observed fields, it is correcting for the combined hydrologic model and weather observational error and not for the hydrologic error due to forecasted weather variable uncertainty.

#### 4) Multimodel approach

A natural extension of the earlier incorporation of the CLM in 2003 and the SDM in 2004 was to combine the two models under a multimodel framework, capitalizing on the strengths and weaknesses of both approaches as a function of forecast lead time based on the two models’ past performance.

The multimodel approach is automatically recalibrated daily utilizing both the updated CLM and SDM (AR error corrected) hydrologic model hindcasts as follows:

Rank models based on their total hindcast residual error (square difference between modeled and observed discharges).

Evaluate resultant residuals of the highest-ranked model and the multimodel using the BIC to avoid overcalibration. If BIC of the multimodel is minimized, then use the resultant regression coefficients {

*a*_{1},*a*_{2}} to generate “multimodel” forecasts and hindcasts; otherwise, use the highest-ranked model.

It should be noted that by using hindcasts derived from observed weather fields, as opposed to the forecasted fields, we are trying to isolate and minimize hydrologic model error in the multimodel process. Multimodel discharge forecasts are generated by applying the regression equations to the same ensemble member of each hydrologic model achieved by tracking each ECMWF ensemble member through both CLM and SDM hydrologic models and then individually combining them into a multimodel ensemble member.

The multimodel regression coefficients are essentially weighting coefficients and, as such, can be used to evaluate the skill of each hydrologic forecast model relative to one another. Figure 5 plots the actual regression coefficients (or weights) applied to each of the two hydrologic models (CLM, red; SDM, blue) as a function of forecast lead time for the Ganges (upper panel) and the Brahmaputra (lower panel) for a single day (15 October 2003). Because the Ganges catchment has a longer rainfall-response time scale than the Brahmaputra,^{3} for the particular meteorological and hydrological state of 15 October 2003, the CLM outperforms the SDM for the Ganges for short lead times. However, the performance reverses at longer lead times as the importance of modeling the longer time-scale groundwater response captured in the physics of the SDM begins to dominate. For the Brahmaputra (lower panel), the CLM outperforms the SDM at all lead times, requiring further investigation. Figure 4 gives an example of one argument why the CLM has shown good skill for severe flooding events. As discussed previously, the CLM should become more effective as large-magnitude precipitation events’ spatial scale approaches that of the catchment itself, as seen preceding the Brahmaputra 26 July–6 August flooding event, especially evident in the ECMWF forecasts of the top panel.

### d. Step IV: Generation of probabilistic forecasts (gridded box, Fig. 2)

We have argued that by driving the hydrologic multimodel with calibrated ECMWF ensemble weather forecasts, we can project the effect of weather forecast uncertainty onto hydrologic forecast uncertainty. How important are the remaining sources of uncertainty?

Figures 6a and 6b show the operational multimodel discharge forecasts using all 51 ECMWF EPS members (colored lines) compared to observations (solid black lines) for the Brahmaputra (year 2004) and Ganges catchments (year 2006), respectively. Notice how in the day-1 lead-time forecasts for both of these catchments [Figs. 6a(i) and 6b(i), respectively], it appears that there is a single forecast, even though all 51 ensemble members are plotted, showing that numerical weather prediction (NWP) uncertainty provides a negligible source of discharge forecast uncertainty. However, as the lead times increase to 4, 7, and 10 days, the ensemble bundle widens further, showing that NWP uncertainty contributes more and more uncertainty to the hydrologic forecasts for longer and longer lead times. Contrasting the same lead-time Brahmaputra and Ganges ensembles dispersions, notice how the ensemble bundle is much narrower for the Ganges as a result of its longer transportation time for water to migrate through the catchment.

If all sources of uncertainty were properly accounted for, then the spread of the discharge ensemble would provide a representation of forecast uncertainty, and the observed discharge would generally fall within the ensemble bundle.^{4} However, in Fig. 6, the observed discharge is noticeably outside the ensemble bundle for a large fraction of days for both river basins, in particular for the Ganges. This points out how NWP forecast ensemble spread is not sufficient to account for discharge forecast uncertainty; other hydrologic modeling errors are of significance (e.g., watershed initial state uncertainties, hydrologic model deficiencies) and clearly they need to be accounted for, and their relative significance compared to NWP uncertainty varies for different lead times and catchments.

Our approach to optimally account for all sources of forecast uncertainty is to separate the effect of weather forecast uncertainty on hydrologic forecast uncertainty from that of other hydrologic modeling and initial state errors. This retains the useful information in the dynamical weather forecast dispersion. Hindcast errors of the remaining sources of uncertainty are then generated and conditionally sampled based on hydrologic conditions at the time of forecast initialization, which also provide an error correction of any residual forecast bias. These two sets of uncertainties are then combined to create a weighted ensemble set of hydrologic forecasts that provide a representation of the expected total hydrologic forecast error. These steps are represented schematically in Fig. 7, which will be discussed in more detail later.

For a given forecast day with a given current observed discharge, state of the catchment, and weather forecast, symbolically what we seek is a final CDF of the future discharge *Q _{f}* , at a given forecast lead time

*f*, that will tell us what the probability is that the river’s discharge will exceed a severe flooding level:

where *F _{E}*(

*Q*≤

_{f}*q*) gives the (cumulative) probability that the future discharge

_{f}*Q*will be less than the specific value

_{f}*q*,

_{f}*F*

_{1}(

*Q*≤

_{f}*q*|

_{f}*q*) represents the (cumulative) probability that

_{m}*Q*≤

_{f}*q*given the discharge forecast model forecasts the value

_{f}*q*with no weather forecast uncertainty, and

_{m}*F*

_{2}(

*Q*≤

_{m}*q*|

_{m}*p*) gives the (cumulative) probability that the modeled discharge

_{f}*Q*will be less than the specific value of

_{m}*q*given the NWP model forecasts the specific weather (precipitation) fields

_{m}*p*. In this representation, we have then decomposed the future discharge total uncertainty into two parts: one part due to hydrologic uncertainty (

_{f}*F*

_{1}) and second part of the total error due to weather (precipitation) forecast uncertainty (

*F*

_{2}).

Although the uncertainty calculation given by Eq. (16) in theory represents integrations over continuous functions, in practice, the determinations of both *F*_{1} and *F*_{2} are performed here by first discretely representing their associated PDFs ( *f*_{1} and *f*_{2}, respectively) by (weighted) ensembles and then discretely convolving the two PDF’s together:

where *π* represents a parameter set for which *f*_{1} is dependent on (which includes *q _{m}* itself),

*q*in the last equality is the discharge derived by running each weather forecast ensemble

_{mi}*p*through the hydrologic model, and is a set of conditionally selected hindcast errors (ɛ

_{fi}*=*

_{jj}*q*−

_{oj}*q*) between observed

_{Hj}*q*and model hindcast discharge

_{o}*q*(derived using observed weather fields). The derivation of the final CDF

_{H}*F*(represented by the integration ) is done by a simple linear interpolation of the cumulative probability between each ensemble member. The specific details of this process are described next with reference to Fig. 7, with further details in Hopson (2005).

_{E}I: The first step is to generate the PDF

*F*_{2}of Eq. (17) that represents the total error due to forecast weather uncertainty, which is done by running the adjusted 51-member ECMWF ensemble weather forecasts through the discharge (multi)model. The transformed ensembles represent the discretized PDF version of*F*_{2}that we seek, separating hydrologic uncertainty from the uncertainty due to probabilistic weather forecasts. These ensembles are shown schematically in the first block of Fig. 7, where each discharge ensemble member*Q*is equally weighted with probability (although only seven ensembles are shown here). This discretized PDF is saved and is combined with_{p}*F*_{1}in Fig. 7 (IV) described later.IIa: The first step in producing

*F*_{1}is to generate an historic simulated discharge time series using the current updated hydrologic multimodel for each forecast lead time in an identical fashion to how the forecasts are made, except using “observed” satellite–rain gauge precipitation in place of weather forecasts. By doing so, we separate weather forecast variable uncertainty from the hydrologic uncertainty. This multimodel time series is then differenced with the historic rating curve–derived discharge observation record to generate the desired multimodel residual time series shown. Under the assumption that these historic errors (a nonparametric error “prior”) are statistically representative of future errors,^{5}this set of residual errors empirically accounts for all forecasting errors except for errors in the weather variable forecasts.^{6}IIb: Because of this dependence of hydrologic uncertainty on a number of measurable factors, the next step is to make a refinement on the hydrologic error prior by making a conditional selection of the hindcasted error ensembles that occurred within a similar region of “state space” as that for the system at the time of forecast initiation.

To make this conditional selection of past error, we first must specify the “state” of the system (i.e., hydrologic error). In this application, we determine this state by a simple set of selector variables: the last observed hydrologic error

*e*(observed minus hindcast discharge), along with the slope (with respect to time)_{H}*s*(*e*), and curvature_{H}*c*(*e*) of this error, at the time of forecast (hindcast) initiation (the slope and curvature are discretely calculated)._{H}^{7}With this set of state variables, steps Figs. 7(IIc)–7(IIf) below describe the conditional selection procedure, in reference to Fig. 7(II), where the selection (feature vector) in this figure is derived from the hindcasted hydrologic error time series.

- IIe: The next step is to compare the current (and forecasted) values of the feature vector (
**V**) with the feature vectors in the historic record (_{E}**V**), and those vectors (i.e., days_{i}*i*) most similar to the current conditions for which they are selected. The measure used in the selection process is the Mahalanobis distance*M*(see Davis 1986; Yates et al. 2003 for other applications). Specifically, we calculate for each day_{i}*i*, selecting out the subset of days for which*M*is lowest, where 𝗦_{i}^{−1}is the inverse of the covariance matrix of the selector variables themselves. Use of*M*diminishes the weighting of the selector variables based on the value of their correlation._{i} IIf: The subset of historic days for which

*M*are lowest is then selected out. The number taken is_{i}*n*=*N*, where*N*is the total number of historic days available. Referring to the visual schematic of Fig. 7(II), here the most recent residual string (the piece of the residual time series that precedes the solid vertical line of the forecast horizon) is compared with the rest of the historical residual error time series via their value, slope, and curvature, and “*n*nearest neighbors” for which they are selected. The arrows highlight the residual strings (dashed) that follow favorable comparisons in the historical record. These dashed residual time series then represent*forecasted*error estimates of the hydrologic error of the*current*hydrologic forecast, and as such, they forecast the hydrologic error into the forecast horizon, as shown, where the dashed forecasted model residuals are appended to the most recent residual error at the point of the forecast horizon, providing an ensemble estimate of hydrologic model error.IIIb: From this weighted set of hindcast initialization days

*i*, the associated*hindcasted*hydrologic error that occurs at day*i*+*f*days is then assigned the same weighting. The result is that we now have constructed a*forecast*of the estimated PDF of the hydrologic error that will occur at day*t*+*f*.IV: The final step is to convolve the forecasted PDF of hydrologic error (steps II and III) with the discharge forecasts derived from weather forecast uncertainty (step I). This is done by mapping the

*inverse*of the hydrologic error PDF onto the PDF of forecast input uncertainty (*f*_{2}); this is done by*subtracting*the value of each hydrologic error ensemble from each weather forecast ensemble–derived discharge ensemble*Q*, and we then end up with a new set of ensembles with size [51 ×_{p}*n*], and each new ensemble has a probability given by [(1/51) ×*W*]. This set of ensembles represents a combined PDF of both weather forecast uncertainty along with hydrologic error uncertainty, and it also has simultaneously provided an additional error correction to the forecast._{i}

To see how this last step provides a correction, recall that the hydrologic error ensembles represent error forecasts. So, we want to remove this error forecast from the final discharge forecast, and in essence, the hydrologic error PDF then represents a final hydrologic error correction. From this PDF, the final discharge forecast CDF *F _{E}* is calculated.

## 3. Results

Figures 8 and 9 and the online supplemental figures (available at http://dx.doi.org/10.1175/2009JHM1006.s1) display discharge forecasts for the Brahmaputra at Bahadurabad and the Ganges at Hardinge Bridge based on the forecast system described in section 2. All forecasts shown are operational forecasts disseminated within the respective years with the following caveats: 1) all rating curves have been objectively recomputed and the forecasts rerun using the new rating curve–derived discharges; and 2) the operational forecasts for 2003 and part of 2004 were based only on the CLM model, they did not utilize the TRMM and CMORPH precipitation estimates, and they did not include the final error accounting algorithm of step IV, but they have been rerun to include these features as shown in the figures.^{8}

Of the years of operational forecast dissemination discussed in this paper (2003–07^{9}), the Brahmaputra reached flood stage five times, in July 2003 (3–6; 9–17), July 2004 (10–26), and two distinct periods in 2007 (26 July 26–6 August; 8–17 September), with the latter shown in Fig. 8. The Ganges, on the other hand, only briefly exceeded critical flood level once in 2003 (19–21 September), shown in Fig. 9, but extensive flooding never occurred. All other operational years are given in the figures in the online supplement. In all these figures, the left column (I) of each panel summarizes the 5-day discharge forecasts, and the right (II) shows the 10-day discharge forecasts (Note that 1- to 10-day forecasts were operationally produced and disseminated but only 5- and 10-day forecasts are shown here for brevity). For presentation purposes, the forecasts have been sorted, removing the useful temporal covariability of the ensemble traces. The dates along the abscissa denote the forecast day. The bold black line shows the observed discharge verification curve. The distance between each plot’s right edge and the vertical black line is a metric of the forecasting lead time.

In Figs. 8 and 9, and those in the electronic appendix, the discharge ensembles [row (i)] were used to calculate the forecast 95% (red) and 50% (gold) confidence intervals [row (ii)]. Perhaps the most important metric is the probability that the rivers’ flow into Bangladesh will exceed its respective danger level (corresponding to discharges of approximately 62 000 and 55 000 m^{3} s^{−1} for the Brahmaputra and Ganges, respectively). These 5- and 10-day lead-time forecasted danger-level exceedance probabilities are shown in row (iii) for each river. These curves correspond in essence to the last step of the technical forecasting process (step V of Fig. 2).

To highlight further the error growth with lead time, in Fig. 10 we zoom in on the 10-day lead-time forecasts during the 2004 severe flooding event and isolate two days for which operational forecasts were issued (although forecasts were issued every day) during periods of extreme flood onset and recession. Shown in Figs. 10(v) and 10(vi) are plume plots of the ensemble forecasts and derived danger-level exceedance probabilities from two separate issuing dates out to the end lead time of 10 days. Vertical black lines bound the severe flooding periods. The two July-issued forecasts shown were chosen specifically for days when the observed hydrograph was decreasing on the issuing date (for which a simple “trend” forecast would have no skill). The first forecast [Fig. 10(v)] was issued 6 days before the discharge crossed the critical threshold and was forecasting 31% likelihood of exceedance for that date, but it was forecasting 66% likelihood at 7 days, and 84%–87% at 8–10 days. The second plume forecast shown [Fig. 10(vi)] points outs the forecast system’s ability to estimate flood duration. Again, as the observed hydrograph was descending, the forecasts instead warned of a second flood peak (maintaining severe flood probabilities of 81%–100% during this above-critical-level period).

Overall, both rivers’ discharge was well forecast at the 5-day and 10-day limit for the Brahmaputra but appeared to be conservative in forecasting the precise date of crossing danger-level thresholds. For the Brahmaputra, the flood forecasts of 2003 (refer to the electronic appendix) showed only moderate skill. However, the major flooding in 2004 was well forecast at 10 days, including the double peak and the duration of the forecast (Fig. 10). In particular, using the 10-day forecast lead-time metric (the gap between the vertical line and the plots’ right edge) in Fig. 10, we note that the first July 2004 flood wave was being forecast as the hydrograph was rapidly descending (continuously for eight days), implying a strictly discharge time series forecast model would not have captured this event. No floods were forecast during 2006 and none occurred. For the first major flood wave of 2007 (Fig. 8), the forecast amplitude was lower than occurred, but the timing and forecast critical-level probabilities were skillful. The second 2007 flood wave was generally well forecast. Overall, except for moderate skill in 2003, all major flooding events were well forecast at 10-day lead times with high probability of occurrence and with well-captured peak amplitudes and dates of flood recession. Notably, during the 5-yr period, there were no false positives of major flooding for either basin.

For the Ganges, the forecasts were generally less skillful than the Brahmaputra, although the modeling system has not been tested yet for severe extended flooding events. Low probabilities of short-duration flooding occurred throughout the forecast record, but no high probability false-positive flooding events were forecast. Overall though, forecast skill begins to lose utility past five days and as discussed previously, one reason for this may be that the forecast system presented here does not explicitly account for the Ganges’s water management (unlike the generally “naturalized” Brahmaputra) despite the Jian et al. (2009) study. Current improvements to the forecast system are being implemented to account for this management, a topic of future work.

All sources of discharge forecasting uncertainty were taken into account in each of the forecasts shown in Figs. 8 –10. Modeling the effect of weather forecast uncertainty, each of the 51-member ECMWF ensemble precipitation forecasts was run through the multimodel (using 1-, 2-, … , 5-day precipitation forecasts for the 5-day discharge forecasts and using 1-, 2-, … , 10-day precipitation forecasts for the 10-day discharge forecasts) and then merged with hydrologic and weather observation uncertainty, as discussed in step IV. Although not shown, there is a great deal of uncertainty—as depicted by the ensemble spread—in even the 1-day precipitation forecasts. However, after the transformation through the modeled Brahmaputra and Ganges watersheds with their basin time scales of order 5 days and 13 days, respectively, 1-day forecasts of precipitation only weakly affect the 1-day discharge forecasts. However, this is not the case for the 10-day discharge forecasts, where the weather forecast’s uncertainty clearly affects the discharge uncertainty. The probabilistic forecasts, so produced, provide the necessary information to decide if there is “sufficient risk” upon which to take action.

Objectively, how successful was the forecast system at providing skillful, calibrated forecasts? The rank histograms shown in Fig. 11 (Brahmaputra only; Ganges not shown, but they give similar results) provide one way to assess the calibration of the ensemble forecasts’ own uncertainty estimates. These panels were generated for the final calibrated 51-member ensemble (with 52 intervals) 5- and 10-day forecasts. The red dashed lines are 95% confidence intervals (CIs) on the expected counts in each bin for a perfect model. For 52 intervals, we would expect roughly 2–3 interval counts to fall outside the CI; instead, roughly 5–6 intervals for both 5- and 10-day forecasts for both basins lie outside the CI (with similar results for other forecast lead times). Although the forecast calibration is not that of a perfect model (at 95% confidence), overall the rank histograms are close to being statistically uniform.

The top panel of Fig. 12 provides an objective Brier scores (BSs) assessment of forecasting skill for the event of Brahmaputra discharge exceeding critical level, corresponding to a discharge of 62 400 m^{3} s^{−1}, or the 94th percentile of the climatological record (*q*_{0.94}). (BSs of the Ganges are not provided because only three days in the forecasting period exceeded critical level.) The results are presented as a function of forecast lead time for the forecast BSs (solid black) and 95% CI (black dash), along with a persistence forecast BS CI (red dash) and a climatological forecast BS CI (green dash). A persistence forecast refers to using the river discharge in day *t* to forecast the river discharge in day *t* + *f*. For the BS, a climatological forecast is given by the probability to exceed the river stage danger level on any given day during the monsoon season (calculated to be 0.06 for the Brahmaputra). The CIs are generated by “bootstrapping” (Huber 1981) the days (and the associated forecasts and observations) in the forecast period (years 2003–07) but preserving the temporal covariance of the time series (which widens the CI). Forecast skill analysis is performed on each bootstrap sample (here, performed 2000 times) and the results are ranked to determine the values corresponding to the 2.5 and 97.5 percentiles (giving the 95% CI). The discharge forecasts are superior to both climatology and persistence at all lead times, with only minimal-to-moderate overlap of their CI. [If one calculates the conditional ranking of the BS for each bootstrap sample (not shown), then there is more than 99% certainty the forecasts are superior to both.]

We further examine these comparisons through Brier skill scores (BSS) in the bottom panel of Fig. 12, where the BSS is calculated as

where BS* _{f}* is the Brier score for the forecasts, BS

_{ref}is for a “reference forecast” [here, either “climatology” (green) or “persistence” (red)], and BS

_{perf}for a “perfect” forecast. The most likely BSSs (solid lines) show the forecasts outcompete both persistence (red solid and dashed CI) and climatology (green solid and dashed CI) at all lead times at >60% improvement over these respective reference forecasts.

Further skill score assessments were done for both the Brahmaputra and Ganges using other measures. Figures 13 and 14 show skill scores for the root-mean-square error of the ensemble mean forecasts [where RMSE is substituted for BS in Eq. (23)] and continuous rank probability skill score of the full ensemble [CRPSS substituted for BS in Eq. (23)] of the Brahmaputra (top panels) and Ganges (bottom panels) river discharge forecasts, respectively. The persistence forecast referred to in these figures is as defined earlier. The RMSE climatological forecast is the mean of the monsoon season historical discharge values, whereas the CRPS climatological forecast is given by the whole PDF of the same (i.e., not just the mean value). Again, the Brahmaputra forecasts (top panels of both figures) significantly outperform both climatology (45%–98% improvements) and persistence (55%–80%). The Ganges forecasts (bottom panels, both figures) also show significant skill gains over climatology (60%–98% improvements) and persistence, except for the RMSE of the 24-h Ganges forecast where we see only slight improvements over persistence. [However, the conditional ranking of the RMSE for each bootstrap sample (not shown) gives >67% certainty of the 24-h forecast outperforming persistence.]

## 4. Conclusions

In this work, we have presented an approach to produce probabilistic river discharge forecasts in the medium range for two large-scale river basins. The ability to produce accurate and timely forecasts at these spatial scales depends primarily on three factors: the availability of quality ensemble weather predictions [produced at a few centers, such as ECMWF, NOAA/National Centers for Environmental Prediction (NCEP), etc.], the availability of a range of satellite precipitation products (TRMM, CMORPH, etc.) to initialize modeled watershed states and to reduce systematic error in the meteorological forecasts, and a method of incorporating these data sources into a hydrological modeling system. In addition, an infrastructure to disseminate and utilize the forecasts is essential. The CFAB methodology presents an operational system that bridges the gap between the production of meteorological observations and forecasts, and the delivery of a useable hydrological product.

The CFAB system was developed in its present configuration because of the absence of upstream streamflow data beyond the borders of Bangladesh. The system depends critically on current and historical discharge data at the India–Bangladesh border for calibration and forecast data assimilation. The quality of the forecasts shown in Figs. 8 –10 may speak to the utility of similar schemes applied to other regions where upstream data may also be inaccessible or just not measured. However, if river flow measurements higher up in the catchment were available and could be routed downriver to the forecast location, errors in rainfall–runoff modeling (due to imprecise satellite-derived precipitation estimates and hydrologic modeling) could be reduced, potentially leading to significant improvements in forecasts at shorter lead times. However, hydrologic runoff process modeling coupled to statistically rendered weather forecast information (as per step II in Fig. 2) would play an increasingly greater role (and in situ river data less so) as the forecast lead times increase, or as catchment flow-through time scales decrease.

CFAB began producing operational 1–10-day advance operational flood forecasts in 2003, with the system presented here completed in 2004. During July 2004, the system produced skillful real-time forecasts of the severe Brahmaputra flooding event, motivating the establishment of a pilot program in 2006 to disseminate the forecasts down to vulnerable citizens in five selected districts with a total population of approximately 110 000 (Fig. 1, right panel; CEGIS 2006; Webster et al. 2010). During 2007, the Brahmaputra flooded severely twice, and the dissemination network was activated to warn vulnerable citizens, leading to the early evacuation of many thousands of people and livestock. As tragic as these floods were to the country as a whole, having the forecast system described here be of service to the people of Bangladesh was tremendously gratifying.

## Acknowledgments

ECMWF has enthusiastically supported CFAB since its inception. In particular, we thank Philippe Bougeault, Roberto Buizza, the late Tony Hollingsworth, Dominique Marbouty, Martin Miller, and Tim Palmer for specific assistance. The Asian Disaster Preparedness Centre has directed the in-country work, and we are especially indebted to Sri. A. R. Subbiah and Sri. S. Ramasamy for their inspirational enthusiasm and dedication in helping the poorest of the poor in South Asia. We are also appreciative of the partnership with the Flood Forecasting and Warning Centre and the Bangladesh Meteorological Department, along with the Center for Environmental and Geographic Information Services, the Institute for Water Modelling, and the many other stakeholder institutions within Bangladesh for supporting this project. We would also like to thank Keith Beven for his very valuable assistance in early 2003 as the project went operational, and thanks to the inspiration of Gilbert White, whose life’s work focused on scientific research meeting the needs of the societal greater good. Initial support for the project was given by USAID under Grant USAID/OFDA AOT-A-00-00-00262-00 (2000–03) and most recently by CARE Bangladesh (2005–07). Basic science has been supported by National Science Foundation Grant ATM-0531771. We also appreciate the support of Georgia Tech’s School of Earth and Atmospheric Sciences, which provided the means for the project to continue during 2005, and the Advanced Study Program and the Research Applications Laboratory of the National Center for Atmospheric Research.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**(

**,**

**,**

**,**

## Footnotes

+ Current affiliation: Advanced Study Program and Research Applications Laboratory, National Center for Atmospheric Research, Boulder, Colorado

*Corresponding author address:* Thomas M. Hopson, NCAR, P.O. Box 3000, Boulder, CO 80307-3000. Email: hopson@ucar.edu

^{1}

Years for which severe flooding occurred in one or both basins.

^{2}

However, Mirza (2003) concluded for the same gauging locations presented here that year-to-year variation in river stage for any common threshold discharge was minimal.

^{3}

The (calibrated) Ganges SDM upper-layer (*t*_{s1}), lower-layer (*t*_{s2}), and unit hydrograph (*t _{c}*) time constants are 4.2, 20.7, and 1.9 days; for the Brahmaputra, they are 0.3, 8.5, and 0.3 days.

^{4}

Specifically, the observed discharge should have equal probability (here, ) of falling between any two ensembles or above (below) the largest (smallest) ensemble.

^{5}

Because such residual errors are often generated from a “tuned” discharge model, strictly speaking the error estimates should come from cross validation or be inflated by calculated prediction intervals (see, e.g., Helsel and Hirsch 1992).

^{6}

Although not discussed in detail here, we argue this prior also includes an implicit (although under) estimate of discharge rating curve error.

^{7}

Note that there are many other hydrologic states that could also be utilized here, such as forecasted or observed catchment-averaged precipitation, soil moisture, forecast or observed discharge, among others.

^{8}

Also, the Akaike Information Criteria (Akaike 1974) has been replaced with the BIC.

^{9}

Operational forecasts were not made for 2005 as a result of a lapse in funding.

* Supplemental information related to this paper is available at the Journals Online Web site: http://dx.doi.org/10.1175/2009JHM1006.s1.