Toward Improving Medium Range Ensemble Streamflow Forecasting—Dissective Analysis of Predictive Skill for the Hydrologic Ensemble Forecast Service (HEFS)

1. Introduction

Improving streamflow forecasting entails reducing total uncertainty in streamflow prediction such that the forecast accuracy is increased while the forecast attributes meet the needs of the user decision making and support. In such efforts, headwater basins are of particular interest as they are generally the largest source of uncertainty. Total uncertainty at headwater locations comprises future input uncertainty and hydrologic uncertainty (Krzysztofowicz 1999, WMO 2025). Future input uncertainty refers to the uncertainty in the future forcings such as precipitation and temperature used as input to the hydrologic models. Hydrologic uncertainty refers to structural, parametric and initial condition (IC) uncertainties in the hydrologic and water management models used as well as anthropogenic uncertainties associated with unknown or partially known human control and alterations of movement and storage of water. Accordingly, streamflow forecast systems and processes include components and elements that are aimed at reducing the above uncertainties individually or collectively (Demargne et al. 2014, Troin et al. 2021).

Ensemble forecasts aim at describing higher-order statistical moments of the predictands beyond central tendency while maximizing predictive skill. Because generating such forecasts requires dynamical-statistical coherence, ensemble forecasting lends itself to holistic design, assessment and improvement of the forecast systems and processes. To that end, the Hydrologic Ensemble Forecast Service (HEFS, Demargne et al. 2014; see Figure 1) of the US National Weather Service (NWS) was designed to address input uncertainty with the Meteorological Forecast Ensemble Processor (MEFP, Schaake et al. 2007, Wu et al. 2011, NWS 2019ab), the IC uncertainty with data assimilation (DA) (NWS 2010) and other hydrologic uncertainties with streamflow post processing (Seo et al. 2006, Alizadeh et al. 2020). To date, however, only the MEFP has been implemented across all River Forecast Centers (RFC), leaving large room for potential improvement (Kim and Seo 2025ab).

An important consideration in the design and implementation of an ensemble streamflow forecast system is the availability and quality of streamflow assimilation. Of all hydrologic uncertainties, only the IC uncertainty is lead time-dependent. Because soil moisture and other state variables are rarely observed at the catchment scale, they are typically inferred from real-time observations of streamflow and forcing variables via hydrologic models. If streamflow assimilation is included in the ensemble forecast system, it may not be necessary to include real-time streamflow observations as a predictor in streamflow post processing. The implication for the HEFS, for example, is that one may use the MEFP or its recent enhancement, the conditional bias penalized regression (CBPR)-aided MEFP (Kim et al. 2025, Kim and Seo 2025a), for streamflow post processing as well. The MEFP combines multiscale probability matching, regression (in normal space), downscaling and the Schaake Shuffle (SS, Clark et al. 2004, Wu et al. 2018) for statistical post processing of precipitation forecasts from numerical weather prediction (NWP) models. The CBPR-aided MEFP uses CBPR in place of ordinary least squares regression (OLSR) to ameliorate under-forecasting of heavy-to-extreme precipitation by reducing Type II error (i.e., false negative) for large events but at some expense of increasing Type I error (i.e., false positive) for small events. The use of a common post processing method for both precipitation and streamflow would not only simplify and streamline statistical modeling but may also improve parsimony for reduced data requirements.

With the above backdrop, this paper aims at addressing two objectives. The first is to describe the full HEFS in which streamflow assimilation and streamflow post processing operate in tandem as designed at inception (see Figure 1). The second is to dissect the performance of the full HEFS toward maximizing the collective predictive skill for medium range (~14 days) ensemble streamflow forecasting. A similar study has recently been reported by Zhang et al. (2025) for prediction of mean daily flow at the outlet of the 15,876 km² upper Huai River Basin in east China using the ensemble Kalman filter for streamflow assimilation, Bayesian Joint Probability method for post processing of precipitation forecast, and Quantile Regression for post processing of streamflow forecast. This study differs significantly from the above in two important ways. First, the spatial scale of the catchments of interest is much smaller in this work, ranging from 190 to 2,400 km² with a median of 985 km². In addition, the predictand of primary interest is instantaneous discharge at a 6-hr timestep in line with the operational practices at the RFCs. Second, the relatively small space-time scale associated with the predictands of precipitation and streamflow means that predictability is rather limited (Kim et al. 2025, Kim and Seo 2026), and that Type II error is also an important source of uncertainty particularly for all-important heavy precipitation and high flows (WMO 2025). Accordingly, we use the CBPR-aided MEFP for post processing of precipitation forecasts from NWP and DA-aided streamflow predictions, and the fixed-lag smoother application of the adaptive conditional bias-penalized ensemble Kalman filter for streamflow assimilation (AEnKF, Shen et al. 2022ab, Seo et al. 2022, Kim and Seo 2025b). We then assess the flow magnitude-dependent contributions of hydrologic modeling, NWP precipitation forecast and its post processing, assimilation of streamflow and streamflow post processing to the collective predictive skill of the resulting ensemble streamflow forecast in reference to climatological ensemble streamflow forecast.

This work advances assessment and understanding of the performance of a full ensemble streamflow forecast system for headwater basins and the contribution of each component to the collective predictive skill in short to medium ranges. Whereas the details are given in the context of the HEFS, the general approach and methodology apply to similar ensemble streamflow forecast systems that partition input and hydrologic uncertainties. This paper is organized as follows. Section 2 describes the study locations, data, models and tools used. Section 3 describes the methodological framework, the simplified streamflow post processor and the evaluation approach. Section 4 presents the results and discussions. Section 5 provides conclusions and future research recommendations.

2. Study Locations, Data, Models and Tools Used

The study area comprises 43 headwater locations in the HEFS Testbed (Lee et al. 2018, 2020) in the service areas of 11 RFCs (see Figure 2 for the map and Supplementary Data (SD) 1 for geographical names and hydrologic attributes). In Figure 2, the forecast points are identified by their 5-character NWS abbreviations. The same locations were used in Kim and Seo (2025a) from which 20-yr ensemble precipitation hindcasts covering 2000 to 2019 are available from the NWS Office of Hydrologic Prediction (OWP). The precipitation hindcasts were generated with the CBPR-aided MEFP forced by the Global Ensemble Forecast System Version 12 (GEFSv12) ensemble mean precipitation hindcasts (Guan et al. 2022, Hamill et al. 2022). The observed mean areal precipitation (MAP) and mean areal temperature (MAT) data are from the RFCs. Climatological mean areal potential evapotranspiration (MAPE) is used for both observed and future MAPE. The observed instantaneous streamflow data are from the US Geological Survey.

The hydrologic models used are SNOW-17 for snowmelt (Anderson 1976), Sacramento for soil moisture accounting (SAC, Burnash et al. 1973) and unit hydrograph for surface runoff routing (UHG, Chow et al. 1988). All models are run at a 6-hr timestep even though some RFCs use shorter timesteps for fast-responding basins. The use of the same timestep for all locations is to allow skill assessment under the same hydrologic assimilation cycle of 6 hours. The data assimilated into SAC and UHG include streamflow, MAP and MAPE where MAP comprises rainfall and snowmelt. For SNOW-17, no DA was performed. Hence, the snow model states, including the snow water equivalent (SWE), are not updated throughout the 20-yr hindcast period which may degrade the quality of snowmelt prediction relative to that attainable operationally via manual DA. To counter the potential degradation of MAP due to lack of snow DA, we assumed clairvoyant future temperature in this work by using MAT in place of temperature forecast. Franz et al. (2014) reports that, for water supply forecasting, the impact of SWE DA relative to DA-less prediction is mixed; whereas DA improved percent bias for ensemble mean forecasts, ensemble forecasts were more skillful without DA. The above finding suggests that the combination of no snow DA and clairvoyant temperature forecast likely inflates the accuracy of snowmelt prediction. Whereas the assumption of clairvoyant temperature forecast is an apparent limitation of this study, it avoids potentially under-representing the predictive skill attainable in the current operations at the RFCs, which is generally not reproducible, and has little impact on the principal findings of this work as shown in the Results and Discussions Section. For medium range streamflow forecasting, SWE and streamflow should ideally be assimilated jointly as a coupled inverse problem. Such coupled DA and the use of temperature forecast warrant separate attention and is left as a future endeavor. Figure 3 describes the dataflow for the hindcast experiment.

Reducing hydrologic uncertainty necessarily depends on the quality of hydrologic modeling, including that of observed forcings. Significant but easily addressable biases or systematic errors in observed forcings or hydrologic model parameters likely inflate the impact of streamflow assimilation and post processing whereas poorly modeled UHG likely deflates the impact of streamflow assimilation (Kim and Seo 2025b). For critical assessment of predictive skill of an ensemble streamflow forecast system, it is hence desirable to reduce such errors as much as possible. To that end, we used AB_OPT (Adjoint-Based OPTimizer, Seo et al. 2009) to correct long-term biases in MAP and MAPE, refine the operational UHG (or newly estimate if necessary) and locally optimize the SAC parameters using the operational settings as starting points. AB_OPT uses adjoint-based optimization to estimate the biases and UHG and the Sequential Line Search (SLS, Kuzmin et al. 2008) to refine the SAC parameters. For the majority of the locations, the magnitude of the adjustment and refinement is relatively small. For some locations, however, the adjustment is rather significant, particularly for UHG due to flow magnitude-dependent variations. Figure 3 shows how AB_OPT was used for the hindcast experiment in this work.

3. Methodology

Ensemble forecasting is a form of Monte Carlo integration for the predictive conditional probability density function (PDF) sought. Hence, an ensemble streamflow forecast system and its composition are determined by this PDF as described below for the HEFS and similar systems.

3.1. Ensemble Streamflow Prediction

For medium range ensemble streamflow forecasting, The HEFS numerically derives the conditional PDF of streamflow valid up to 14 days into the future, given the observed streamflow valid up to and including the current time and the GEFS ensemble mean forecasts of precipitation and temperature valid at the same future times. The above PDF may be written as $f_1\left(q_f|q_o,\overline{p_f}\right)$ where q_f is the future streamflow, q_o is the observed streamflow and $\overline{p_f}$ is the ensemble mean forecasts of precipitation and temperature. Throughout this paper, f( ) denotes a PDF in which only the experimental values of the random variables are indicated for brevity. For example, $f_1\left(q_f|q_o,\overline{p_f}\right)$ abbreviates $f_{{1Q}_f|Q_o,\overline{P_f}}\left(q_f|q_o,\overline{p_f}\right)$ where Q_f, Q_o and $\overline{P_f}$ are the random variables, and q_f, q_o and $\overline{p_f}$ are their experimental values, respectively. In some forecast systems, the predictor for the forcings may be ensemble forecasts from NWP, p_f, rather than their ensemble mean, $\overline{p_f}$, or some other single-valued representation. In such cases, p_f would replace $\overline{p_f}$ with little loss of generality. One may of course choose to model $f_1\left(q_f|q_o,\overline{p_f}\right)$ purely statistically or exclusively via machine learning, but most operational agencies use hydrologic modeling to leverage encapsulated knowledge (Abbott 1999).

Introducing the hydrologic model-predicted streamflow valid at the same future times, s_f, the above conditional PDF may be written via the total probability law as (Krzysztofowicz 1999, Seo et al. 2006):

$$f_1\left(q_f|q_o,\overline{p_f}\right)=\int f_2\left(q_f|q_o,\overline{p_f},s_f\right)f_3\left(s_f|q_o,\overline{p_f}\right)d{s}_f$$

(1)

where s_f is the future streamflow due to the future forcings and the ICs of the hydrologic models. In the above, $f_3\left(s_f|q_o,\overline{p_f}\right)$ may be written in terms of the uncertainty sources modeled:

$$f_3\left(s_f|q_o,\overline{p_f}\right)=\iint f_4\left(s_f|b_f,{i_c,q}_o,\overline{p_f}\right)f_5\left(i_c|b_f,q_o,\overline{p_f}\right)f_6\left(b_f|q_o,\overline{p_f}\right)d{i}_{c}d{b}_f$$

(2)

where b_f is the future input forcings and i_c is the ICs of the hydrologic models. Inserting Eq. (2) into Eq. (1), we have for $f_1\left(q_f|q_o,\overline{p_f}\right)$:

$$f_1\left(q_f|q_o,\overline{p_f}\right)=\int f_2\left(q_f|q_o,\overline{p_f},s_f\right)\iint f_4\left(s_f|b_f,{i_c,q}_o,\overline{p_f}\right)f_5\left(i_c|b_f,q_o,\overline{p_f}\right)f_6\left(b_f|q_o,\overline{p_f}\right)d{i}_{c}d{b}_{f}d{s}_f$$

(3)

Noting that B_f is independent of Q_o and the ICs are independent of future forcings, one may reduce $f_6\left(b_f|q_o,\overline{p_f}\right)$ and $f_5\left(i_c|b_f,q_o,\overline{p_f}\right)$ to $f_6\left(b_f\right|\overline{p_f})$ and f₅(i_c| q_o), respectively. Similarly, noting that i_c and b_f in $f_4\left(s_f|b_f,{i_c,q}_o,\overline{p_f}\right)$ reflect q_o and $\overline{p_f}$ via f₅(i_c| q_o) and $f_6\left(b_f\right|\overline{p_f})$, respectively, one may reduce $f_4\left(s_f|b_f,{i_c,q}_o,\overline{p_f}\right)$ to f₄(s_f| b_f,i_c). Because s_f already reflects q_o and $\overline{p_f}$ via $f_3\left(s_f|q_o,\overline{p_f}\right)$ in Eq. (1), one may reduce $f_2\left(q_f|q_o,\overline{p_f},s_f\right)$ to f₂(q_f| s_f). With the above reductions, Eq. (3) may be rewritten as:

$$f_1\left(q_f|q_o,\overline{p_f}\right)=\int f_2\left(q_f|s_f\right)\iint f_4\left(s_f|b_f,i_c\right)f_5\left(i_c|q_o\right)f_6\left(b_f\right|\overline{p_f})d{i}_{c}d{b}_f$$

(4)

Eq. (4) captures the 4 scientific components of the full HEFS where f₂(q_f| s_f), f₄(s_f| b_f,i_c), f₅(i_c| q_o) and $f_6\left(b_f\right|\overline{p_f})$ represent streamflow post processing, hydrologic modeling, streamflow assimilation and future input uncertainty processing, respectively (see Figure 1). In most operational practices, hydrologic models do not include uncertainty terms. Hence, f₄(s_f| b_f,i_c) is an impulse function which simply maps i_c and b_f to s_f with a probability of unity. Figure 4 depicts the ensemble generation process associated with Eq. (4) where the order of f₅(i_c| q_o) and $f_6\left(b_f\right|\overline{p_f})$ may be reversed with no loss of generality. If desired, one may include parametric uncertainty or multi-model ensembles in the above. Because a unique combination of a hydrologic model and a parameter set defines a unique state space, multi-model ensembles are algorithmically no different from multi-parameter ensembles. In this work, parametric uncertainty or multi-model ensembles are not considered because such processing is unlikely to be operationalized in the very near future due to the computational cost associated with performing DA for each state space and the challenges associates with modeling parametric uncertainty and its interpretation by human forecasters.

Hydrologic modeling is described in Section 2, the MEFP is described in Schaake et al. (2007), Wu et al. (2011) and NWS (2019ab), the CBPR-aided MFEP is described in Kim et al. (2025) and Kim and Seo (2025a) and streamflow assimilation is described in Shen et al. (2022ab), Seo et al. (2022) and Kim and Seo (2025b). Streamflow assimilation is applied in a strongly-constrained mode in this work, i.e., the hydrologic models act as equality constraints in the optimization. Hence, the DA results reflect, and often reveal, deficiencies in hydrologic modeling which may be difficult to detect and identify in forward simulations alone (Lee et al. 2019, Kim and Seo 2026). Below, we describe post processing of DA-aided streamflow predictions via the CBPR-aided MEFP.

3.2. Post Processing of DA-Aided Streamflow Predictions

Streamflow post processing addresses hydrologic uncertainty only. Hence, one may exclude future input uncertainty by replacing the forecast forcings in Eq. (4) with the observed forcings which are assumed uncertainty-free over the prediction horizon:

$$f_1\left(q_f|q_o\right)=\int f_2\left(q_f|s_f\right)\int f_4\left(s_f|i_c\right)f_5\left(i_c|q_o\right)d{i}_{c}d{s}_f$$

(5)

SD 2 shows the numerical integration of Eq. (5) and hence the ensemble generation process for parameter estimation and verification of streamflow post processing. A practical challenge in modeling the streamflow post processor, f₂(q_f| s_f) in Eq. (5), is that, depending on the relative proximity of the ensemble members to the (unknown) true state in the state space, the predictive skill of s_f may vary from one ensemble member of I_c to the next. Hence, it may be necessary for parameter estimation of $f_2\left(q_f|s_f^k\right),k=1,\dots ,n_i$, where n_i is the ensemble size for streamflow assimilation (see Figure 4), to be ensemble member-specific, making it computationally very expensive. Alternatively, one may assume that the predictive skill of ensemble mean of s_f, or $\overline{s_f}$, is comparable to that of $s_f^k,k=1,\dots ,n_i$, and approximate Eq. (5) with:

$$f_1\left(q_f|q_o\right)\approx f_7\left(q_f|\overline{s_f}\right),{\overline{s_f}=\int s_{f}f}_3\left(s_f|q_o\right)d{s}_f=\int s_f\int f_4\left(s_f|i_c\right)f_5\left(i_c|q_o\right)d{i}_{c}d{s}_f$$

(6)

In other words, collapse the ensemble members of S_f, $s_f^k,k=1,\dots ,n_i$, into ensemble mean of S_f, $\overline{s_f}$, and use it as the sole predictor in streamflow post processing, thereby reducing computational cost for parameter estimation by a factor of n_i. The above simplification does not mean that the DA-updated uncertainty in the ICs is ignored. Rather, the updated IC uncertainty is lumped together with all other hydrologic uncertainties into $f_7\left(q_f|\overline{s_f}\right)$. In this work, a single set of hydrologic models (i.e., SNOW-17, SAC and UHG) and a single set of model parameters are used. Hence, one may expect the ensemble members, $s_f^k,k=1,\dots ,n_i$, to be strongly inter-correlated and most of the information content therein to be explainable by a very small number of orthogonal components (Brown and Seo 2009). In such a case, one may expect Eq. (6) to be a good approximation of Eq. (5). Indeed, limited comparative analysis indicates that the impact of the simplification is negligible. If multiple models or multiple sets of model parameters are used, however, Eq. (6) is not likely to be reasonable. SD 2 shows the ensemble generation processes under the above simplification which is to be compared with Figure 4.

When represented by $f_7\left(q_f|\overline{s_f}\right)$ in Eq. (6), post processing of streamflow is completely analogous to that of precipitation forecast represented by $f_6\left(b_f\right|\overline{p_f})$ in Eq. (4). Hence, in principle, one may use the CBPR-aided MEFP not only for precipitation but also for streamflow. The MEFP combines multiscale probability matching for bias correction, linear regression in normal space for statistical prediction, multiscale downscaling via historical ensemble traces and the SS to stitch together the synthetic ensemble traces across different canonical events (CE) while preserving rank correlation in the historical ensemble traces. Contrary to the name, CEs are not events in the hydrometeorological or probabilistic sense but forecast sub-horizons of varying durations and hence equivalent to support in geostatistics. To apply the CBPR-aided MEFP to streamflow post processing, however, a few modifications are desirable as elaborated below.

The first pertains to distribution modeling. In the CBPR-aided MEFP, $\overline{p_f}$ and the verifying observed precipitation are assumed to follow 2-parameter Weibull or gamma distributions (NWS 2019ab). In streamflow post processing, empirical distributions are used instead. Streamflow is strongly modulated by regime-dependent hydrologic processes, fixed boundary conditions across scale and anthropogenic controls and alterations by varying degrees. Hence, the probability distribution of streamflow varies very widely and, unlike precipitation, is often not well described by additive or multiplicative random processes which lead to known parametric distributions (Benjamin and Cornell 2014). The use of empirical distributions means that it is necessary to model the tails outside of the historically observed ranges. The left tail is of little practical concern because the magnitude of flow involved is very small though it may significantly impact certain verification statistics such as the rank histogram. The right tail, on the other hand, has large impact on large to extreme flows, including record-breaking flows, and hence requires care. In this work, we use the hyperbolic model for the right tail end (Deutsch and Journel 1998):

$$F_{Q\ge }\left(q\right)=1-{\displaystyle \frac{q_{\mathit{max}}^{\overline{\varpi }}(1-p_{\mathit{max}})}{q^{\overline{\varpi }}}}\overline{\varpi }\ge 1$$

(7)

In the above, Q is instantaneous flow (observed or simulated), F_Q≥(q) is the cumulative distribution of Q, q_max is the maximum observed or predicted flow, p_max is the probability of less than or equal to q_max and $\overline{\varpi }$ is the user-specified parameter which controls the length of the tail (the smaller, the longer). In this work, we used $\overline{\varpi }$ = 3.25 for all locations based on limited sensitivity analysis.

The second modification pertains to the modeling of temporal variability of the ensemble traces over the forecast horizon sampled from $f_2\left(q_f|\overline{s_f}\right)$ in Eq. (6). Compared to the modeling of probability distribution or conditional bias, temporal variability modeling is of second-order importance (Wu et al. 2018). The quality of temporal variability, however, impacts the quality of multi-daily flow forecasts and hence is likely important for large events. In the MEFP and CBPR-aided MEFP, the SS matches rank correlation of the ensemble traces sampled from $f_6\left(b_f\right|\overline{p_f})$ to that of the climatological precipitation forecast (for graphical illustration, see Figure 3 of Kim and Seo 2025a). In streamflow post processing, the SS is modified to match rank correlation of the post-processed streamflow prediction to that of the DA-aided. The rationale for the above change is that, for large events, the temporal variability of DA-aided streamflow prediction is likely much more realistic than that of the climatological streamflow forecast. Whereas the DA-aided predictions reflect the current hydrologic and hydrometeorological conditions, climatological streamflow forecasts are bound to include many historical traces of low flow which are overly persistent for large events. Below, we provide a brief algorithmic description of the CBPR-aided MEFP as applied to streamflow post processing after Kim and Seo (2025a).

Given the DA-aided ensemble mean streamflow forecast, $\overline{s_f}$, the streamflow post processor samples synthetic traces from $f_7\left(q_f|\overline{s_f}\right)$ for each CE where $f_7\left(q_f|\overline{s_f}\right)$ is modeled using CBPR (Kim et al. 2025). Figure 5 shows the 77 CEs used for streamflow post processing. Though not pursued in this work, the CE selection may be empirically optimized to maximize skill while minimizing the number of CEs used (see Kim 2024 for an example for precipitation). The streamflow post processor then uses the SS to stitch together the CE-specific synthetic traces of streamflow with the temporal variability of the DA-aided ensemble streamflow forecast in the ordinal sense over the forecast horizon of 14 days. The SS starts with the DA-aided ensemble streamflow forecast at a 6-hr timestep valid over the forecast horizon. For a given CE, the ensemble members are aggregated to the duration of the CE and sorted in the ascending order. The streamflow post processor samples synthetic traces from the CE-specific $f_7\left(q_f|\overline{s_f}\right)$ and sort them in the ascending order. Each of the sorted, or rank-ordered, synthetic traces is then assigned the ensemble member identifier in which the DA-aided streamflow forecast of the same rank occurred.

For example, if the CE covers 0-to-96 hrs of lead time (see Figure 5) and the largest DA-aided streamflow prediction in this 4-day period occurred in the 21^th ensemble member for the basin of interest, the largest synthetic trace sampled from the CE-specific $f_7\left(q_f|\overline{s_f}\right)$ is assigned an ensemble member identifier of “21”. The SS then updates each synthetic trace from the latest shuffling of the ensemble streamflow forecast at a 6-hr timestep (recall that this is initially the DA-aided ensemble streamflow forecast) by multiplicatively adjusting the instantaneous streamflow at a 6-hr timestep uniformly within the CE such that the adjusted mean flow is the same as the synthetic trace sampled from the CE-specific $f_7\left(q_f|\overline{s_f}\right)$. The above process is repeated for all CEs beginning with the CE with the smallest Pearson product-moment correlation, or correlation for short, in the bivariate standard normal space between the positive DA-aided ensemble mean streamflow forecast and the positive observed streamflow and ending with the CE with the largest correlation.

The above sequence of operations cumulatively reflects the temporal patterns of granularity in instantaneous streamflow at a 6-hr timestep gleaned from the successive adjustments during the SS while utilizing all available scale-dependent skill in the input DA-aided ensemble mean streamflow forecast. Once the above steps are completed for all CEs, one has replaced the initial traces of the DA-aided ensemble streamflow forecast over the forecast horizon with the same number of ensemble traces of post-processed streamflow. The resulting ensemble streamflow forecast reflects the multi-scale predictive skill in the DA-aided ensemble mean streamflow forecast over different parts of the forecast horizon and share the same rank correlation and temporal variability in the ordinal sense.

The conditional distributions of observed flow, Y, given the ensemble mean of DA-aided streamflow prediction, X, is modelled identically to that of observed precipitation given the ensemble mean precipitation forecast as described in Wu et al. (2011). The only algorithmic difference is that Pr[Y = 0| X = x,X > 0], denoted as c(x) in Wu et al. (2011), is evaluated empirically in this work rather than analytically due to the use of empirical cumulative distribution functions for positive X and Y. For ephemeral basins, intermittency in streamflow is modeled in the same way as precipitation intermittency (Wu et al. 2011).

Lastly, in streamflow post processing, the CBPR coefficient is estimated without magnitude-dependent optimization (Kim et al. 2025), i.e., by minimizing mean CRPS only once using all available pairs of positive DA-aided prediction and verifying observed flow. In magnitude-dependent optimization, the number of pairs of positive observed flow and positive DA-predicted flow are successively reduced by incrementally increasing the conditioning threshold. The rationale for this simplification is that, unlike precipitation forecasts, magnitude-dependent biases in streamflow prediction tend to arise from regime-dependent errors or routing errors. In such cases, the CBPR coefficient from magnitude-dependent optimization likely lacks robustness necessary for use in operational forecasting.

3.3. Evaluation

Eq. (4) provides the framework for evaluation as well; using climatological streamflow forecast as reference, evaluate the incremental contributions of hydrologic modeling, f₄(s_f| b_f,i_c), the use of NWP precipitation forecast and its post processing, $f_6\left(b_f\right|\overline{p_f})$, streamflow assimilation, f₅(i| q_o), and post processing of the DA-aided streamflow predictions, f₂(q_f| s_f) (or $f_2\left(q_f|\overline{s_f}\right)$ with simplification), to the collective predictive skill. Climatological streamflow forecast, or climatology for short, is constructed by sampling observed streamflow over the same prediction horizon of 14 days as the subject forecast from each of the available historical years except the year of the hindcast. Due to many missing streamflow observations, the ensemble size for climatology varies from location to location and from one prediction time to another. The minimum ensemble size used for climatology is 10 to include as many locations as possible.

The primary measure used to assess the predictive skill of the subject forecast in reference to climatology is the CRPSS. The CRPSS is equivalent to the percent reduction in mean CRPS by the subject forecast over the reference forecast:

$${\mathit{CRPSS}}_F={\displaystyle \frac{\overline{{\mathit{CRPS}}_R}-\overline{{\mathit{CRPS}}_F}}{\overline{{\mathit{CRPS}}_R}}}=1-{\displaystyle \frac{\overline{{\mathit{CRPS}}_F}}{\overline{{\mathit{CRPS}}_R}}}$$

(8)

where $\overline{{\mathit{CRPS}}_R}$ and $\overline{{\mathit{CRPS}}_F}$ are the mean CRPS (Hersbach 2000) of the reference and subject forecasts, respectively. The CRPS is a continuous integration of the Brier Score (BS, Brier 1951) weighted by the magnitude of the predictand. Hence, mean CRPS reflects reliability, resolution and uncertainty via the calibration-refinement decomposition of the BS and Type II conditional bias, discrimination and sharpness via the likelihood base rate decomposition (Murphy and Wrinkler 1987). Because reliability and Type II conditional bias are competing attributes over the tails of the distribution, it is necessary to assess at minimum reliability and Type II conditional bias or, alternatively, discrimination (WMO 2025). In this work, we assess reliability using the reliability diagram and discrimination using the relative operating characteristic (ROC) which is insensitive to reliability (WMO 2025).

Large-sample hindcasting of the full HEFS for many locations is expensive. Given that flood and water supply forecasting is often the most important service for the producers of ensemble streamflow forecasts, we limit ourselves to the wettest consecutive 3 months for each location in this work. This means that the predictive skill for the entire year is likely larger than that seen in this work. Comparative verification of the MEFP and CBPR-aided MEFP using 20-yr precipitation hindcasts has shown that 10-fold cross validation results are very similar to the dependent validation results (Kim and Seo 2025a). This is an indication that parametric uncertainty due to interannual variability of precipitation is relatively small for the CBPR-aided MEFP owing to the parsimony of the statistical modeling. Accordingly, we used the 20-yr precipitation ensemble hindcasts generated with the CBPR-aided MEFP in a dependent validation mode in this work available from the NWS OWP. For streamflow post processing, 2-parameter Weibull and gamma distributions are replaced by a combination of empirical distributions and hyperbolic modeling of the upper tail as described above. Hence, streamflow post processing is likely subject to larger sampling uncertainties than precipitation post processing. For this reason, we carried out leave-one-year-out cross validation for streamflow post processing.

3.4. Results and Discussion

Before presenting the results for dissection of predictive skill, we first assess the performance of the streamflow post processor described above under clairvoyant input forcings. In this way, we may validate its performance without being obscured by future input uncertainty. The following questions are to be addressed. Does the streamflow post processor consistently improve accuracy over streamflow assimilation alone? How significant is the margin of improvement? Are the post-processed streamflow ensembles reliable? How does streamflow post processing impact the discriminatory skill of the DA-aided predictions? What are the factors impacting the performance?

3.5. Assessment of Streamflow Post Processing

Figure 6a shows the CRPSS averaged over all 43 locations vs. lead time of the base (i.e., DA-unaided) prediction, or CRPSS_BASE (black), the DA-aided prediction, or CRPSS_DA (red), the post-processed DA-aided prediction, or CRPSS_PP (blue), and the post-processed DA-aided prediction under perfect weak-sense reliability, or ${\mathit{CRPSS}}_{\mathit{PP}}^{\mathit{POT}}$ (purple), for all ranges of observed flow. The reference forecast is climatology. The number parenthesized in the title of each panel indicates the number of locations included. For base and DA-aided predictions, the hydrologic models are initialized with warm states generated by running the models with observed forcings over the entire period of record prior to hindcasting. For the base prediction, the warm states are perturbed based on the observational uncertainties for MAP and MAPE over the assimilation window used for DA. For the DA-aided prediction, the perturbed warm states are updated via AEnKF based on the observational uncertainties for MAP, MAPE and streamflow (Kim and Seo 2026). In this way, we do not unfairly penalize the base predictions for lack of hedging vs. the DA-aided predictions. ${\mathit{CRPSS}}_{\mathit{PP}}^{\mathit{POT}}$ is calculated using the potential mean CRPS, $\overline{{\mathit{CRPS}}_{\mathit{PP}}^{\mathit{POT}}}$, via the mean CRPS decomposition (Hersbach 2000):

$$\overline{\mathit{CRPS}}=\mathit{REL}-\mathit{RES}+\mathit{UNC}=\mathit{REL}+\overline{{\mathit{CRPS}}^{\mathit{pot}}}$$

(9)

where REL, RES and UNC are reliability, resolution and uncertainty, respectively. Figs 6bcd show the results for verifying observed flow exceeding the 90^th, 95^th and 99^th percentiles, respectively.

The CRPSS results vary greatly from one location to another in reflection of the widely varying hydrologic predictability (Kim and Seo 2026). To provide a sense of the variations across all locations, SD 3 shows the location-specific CRPSS overlaid in single plots where the worst-performing locations for CRPSS_BASE are identified by name. Before summarizing the results, two comments are in order. First, the largest conditioning threshold of the 99^th percentile is generally much smaller than flood flow. Hence, the 99^th percentile results provide only a glimpse into how streamflow post processing might perform for flood flows. Second, the conditional CRPSS shown in Figs 6bcd and SD 3 is not a proper score. If used by itself, such scores will lead to frequent poor decisions. Potential for misinterpretation or misuse, however, may be avoided by accompanying them with the unconditional results such as Figure 6a and the relationship between the two (WMO 2025). For example, one may write mean CRPS as:

E [CRPS] = E[CRPS| obs > q_c]Pr[obs > q_c] + E[CRPS| obs ≤ q_c]Pr[obs ≤ q_c]

(10)

where obs is the observed flow, q_c is the threshold flow of choice and Pr[∙] is the probability of occurrence of the event bracketed. Note in Eq. (10) that Pr[obs > q_c] is the familiar exceedance probability. Depending on one’s verification interest, one may use any conditioning event as long as the two conditioning events in Eq. (10) are mutually exclusive and collectively exhaustive.

Figure 6 indicates the following. On average, streamflow post processing improves over streamflow assimilation beyond the first day or so of lead time. Over very short lead times, however, streamflow post processing adds little for most locations. This is because DA-aided predictions are highly accurate where the IC memory is very strong. Statistical correction via post processing, on the other hand, improves accuracy at longer lead times where the IC memory fades and other hydrologic errors take hold. When all ranges of flow are considered, CRPSS_PP is very close to ${\mathit{CRPSS}}_{\mathit{PP}}^{\mathit{POT}}$ (see Figure 6a), an indication that the post-processed streamflow ensembles are generally very reliable in the weak sense. This is because low flows, which dominate the sample size, are highly predictable. For higher flows, however, the post-processed streamflow ensembles are increasingly less reliable as predictability quickly decreases (see Figs 6bcd). Figure 6 also shows that CRPSS_BASE for high flows is significantly larger than that for all ranges of flow. This is because hydrologic modeling is generally more skillful for high flows than low flows relative to climatology which is very poor for high flows. The location-specific results (see SD 3) indicate that the impact of streamflow post processing is mostly positive for most locations but the magnitude of the impact varies very significantly from location to location in reflection of the wide variations in hydrologic predictability. The margin of improvement by streamflow DA and streamflow post processing is particularly large for basins where snowmelt is important due to a combination of longer hydrologic memory associated with snowmelt (Kim and Seo 2026) and increased error in streamflow simulation arising from increased error in MAP which comprises rainfall and (DA-unaided) simulated snowmelt.

Climatological streamflow forecast is free from parametric, structural or anthropogenic uncertainties but has fully saturated IC and future input uncertainties. If the hydrologic uncertainty in streamflow prediction under clairvoyant future forcings is larger than the total uncertainty in climatological streamflow forecast, it is a clear indication that hydrologic modeling warrants attention (Kim and Seo 2026). One may identify such locations by examining location-specific CRPSS_BASE. SD 3 indicates that MFNW1 and ORTW1 in the NWRFC, GTBM3 in the NERFC and CWAC1 and FTJC1 in the CNRFC are such skill-deficient locations. Scatter plots of the DA-aided predictions vs. observed flow for the above locations (see SD 4) show large variations in the quality of the predictions and pronounced sensitivity to the ICs. SD 4 also illustrates how regime-dependent errors may be particularly challenging for streamflow post processing. The first 3 of the above 5 locations have snow fractions over 20% and the first 2 have glaciers (see SD 1). It is very likely that the uncertainties in snowmelt simulation and the lack of glacier modeling contributed to the deficient skill at these locations. Recall that snowmelt is a part of the observed precipitation input to SAC. Hence, even a small but consistent error in snowmelt simulation may greatly deteriorate skill for streamflow prediction particularly for low flows. In addition, snowmelt is strongly seasonal and often very predictable which makes climatological streamflow forecast highly competitive. For the above locations, hydrologic modeling should ideally be addressed first before streamflow DA and streamflow post processing may be applied.

Figure 7 shows the reliability diagrams of the post-processed DA-aided predictions of instantaneous flow at selected lead times of the first 6-hr period in Day 1, the last 6-hr period in Day 4 and the last 6-hr period in Day 14. The thresholds used are the 50^th, 90^th, 95^th and 99^th percentiles of observed flow. In each plot, there are 43 curves corresponding to the 43 locations. Recall that the streamflow post processor samples ensembles from $f_7\left(q_f|\overline{s_f}\right)$ (see Eq. (6)) which is equivalent to the conditional PDF in the calibration-refinement factorization (Murphy and Winkler 1987). Hence, if the assumptions behind the post processor are met, one may expect the post-processed ensemble to be reliable in the strong sense. Figure 7 thus informs how well the assumptions for the multiscale probability matching, regression, downscaling and the SS may be met given the limits of hydrologic predictability.

In Figure 7, the slanted “S” shape in the reliability curves at the shortest lead time is due to heteroscedasticity of the error in the bivariate normal space as explained below. In the MEFP or CBPR-aided MEFP, the predicted flow and the verifying observed flow resulting from normal transformation are assumed to be jointly normal. At very short lead times, the errors in the DA-aided predictions primarily reflect flow-dependent structural and parametric uncertainties particularly when DA is strongly constrained. Such errors tend to vary significantly across different flow regimes in the joint normal space and are often larger over the tails than in the middle of the distribution relative to the ellipsoidal joint normal PDFs. Because the CBPR or OLSR does not model heteroscedasticity, the predicted probability tends to over- and under-prescribe small and large probabilities, resulting in over- and under-spread ensembles in the middle and over the tails of the distribution, respectively. As the lead time increases, the errors become more random and the pattern largely disappears.

Figure 7 indicates that the post-processed streamflow ensembles are generally reliable. The degree of reliability, however, varies significantly from location to location and tends to deteriorate at the lowest and highest thresholds. There are several locations where reliability is noticeably poor. Such results are due primarily to deficient or very challenging hydrologic modeling and cannot, in general, be addressed effectively by streamflow assimilation, which depends on the fidelity of hydrologic modeling, or post processing, which depends on the statistical attributes of the DA-aided predictions. Whereas a more complex streamflow post processor might improve reliability, it is likely to increase the degrees of freedom and hence data requirements for parameter estimation.

Reliability, or Type I conditional bias, competes with Type II conditional bias over the tails of the distribution of the predictand. Recall that AEnKF used for streamflow assimilation addresses Type II error jointly with Type I error. Streamflow post processing without magnitude-dependent optimization of the CBPR coefficient, on the other hand, mainly addresses reliability. Hence, it is necessary to assess the impact of streamflow post processing on the discriminatory skill of the DA-aided predictions. For this purpose, we use the area under the curve (AUC) of the ROC. Figure 8 shows the AUC averaged over all available locations of the base predictions (black), DA-aided predictions (red) and post-processed DA-aided predictions (blue) under clairvoyant future forcings. The conditioning events are the verifying observed instantaneous flow exceeding the 50^th, 90^th, 95^th and 99^th percentiles for the respective locations. Note that some locations could not be included at the highest threshold due to lack of data points. Predictions with perfect discrimination have an AUC of unity. Predictions with no discriminatory skill have an AUC of 0.5.

Figure 8 shows that the AUC of the post-processed DA-aided predictions is nearly identical to that of the DA-aided (i.e., red and blue lines overlap), an indication that streamflow post processing has little impact on the discriminatory skill of the DA-aided prediction. This is not very surprising in that streamflow post processing without magnitude-dependent optimization largely functions as a reliability operator, and that ROC, and hence AUC, is insensitive to reliability. Figure 8 shows that streamflow assimilation using AEnKF significantly improves discrimination over the base prediction for all conditioning thresholds owing to the conditional bias-aware assimilation of real-time observations. To provide a sense of the variations across all locations, SD 5 shows the location-specific results overlaid in single plots. Note that some of the reliability-deficient locations in SD 3 reappear in SD 5, but that there are new locations that are deficient in discrimination, i.e., LNDK1 in MBRFC and CATA1 in SERFC. Similarly to the reliability-deficient locations, discrimination-deficient locations indicate regime-dependent errors in the DA-aided predictions (see SD 4). For these locations, hydrologic modeling should be improved first before streamflow assimilation and streamflow post processing may be applied. Collectively, the results presented in this subsection indicate that the streamflow post processor improves accuracy over streamflow assimilation, renders the DA-aided predictions reliable for most locations and retains the discriminatory skill of the DA-aided predictions. For a number of locations, however, performance is noticeably poor due to deficient or very challenging hydrologic modeling.

3.6. Component-Specific Contributions to Predictive Skill

With the performance of the simplified streamflow post processor assessed, we are now in a position to assess the component-specific contributions to the overall predictive skill in reference to climatology. Table 1 summarizes the forecasts assessed.

In Table 1, CLIMP refers to the ensemble streamflow forecast generated by forcing the hydrologic models with climatological ensemble precipitation forecast. The hydrologic models are initialized with warm states generated by running the models with observed forcings over the entire period of record prior to hindcasting. The warm states are perturbed in the same way as described in Section 4.1. Hence CLIMP represents an imperfect translation of climatological forcings to future streamflow using imperfect hydrologic models with skillful initialization and is equivalent to the long-standing NWS operation Ensemble Streamflow Prediction (ESP, Day 1985) but with perturbed ICs. Climatological ensemble streamflow forecast, on the other hand, represents a perfect translation of climatological forcings to future streamflow using perfect hydrologic models but with skill-less ICs. Climatological ensemble precipitation forecast is sampled from each year of the historical record except the year of the hindcast over the same forecast horizon of 14 days as the subject forecast. The forecast GEFS_CM refers to the streamflow prediction forced by the post-processed NWP precipitation forecast where the post processor is the CBPR-aided MEFP and the NWP forecast is the GEFS ensemble mean. The forecast DA refers to the DA-aided streamflow prediction forced by GEFS_CM. The forecast PP refers to the post-processed DA-aided streamflow prediction.

Figure 9a shows the CRPSS of CLIMP (green), GEFS_CM (black), DA (red) and PP (blue) for all ranges of flow vs. lead time averaged over all 43 locations. The gray line, which is to be compared with GEFS_CM (black) only, represents the CRPSS of the DA-less streamflow prediction forced by the ensemble precipitation hindcasts from the MEFP, rather than the CBPR-aided MEFP, with the same GEFS ensemble mean precipitation forecast as input. The CRPSS results vary greatly from location to location in reflection of the widely-varying predictability of precipitation and streamflow. To provide a sense of the variations of the CRPSS, SD 6 shows the location-specific results and identifies skill-deficient locations. It is important to note that the same 5 skill-deficient locations in SD 3 (see Section 5.3), of which 3 have snow fractions over 20%, reappear in SD 6. Recall that SD 3 reflects hydrologic uncertainty only whereas SD 6 reflects both input and hydrologic uncertainties. The above result suggests that, for those basins where snowmelt is important, the hydrologic uncertainty associated with snowmelt modeling is likely a larger modulator of predictive skill seen in Figure 9 than the input uncertainty associated with precipitation and temperature forecasts. The CRPSS of CLIMP at the lead time of 6 hrs closely reflects the skill of DA-unaided streamflow analysis (i.e., streamflow prediction at the lead time of 0 hrs) and hence reflects the quality of hydrologic modeling, including that of observed forcings (Kim and Seo 2026). The CRPSS of CLIMP with perfect hydrologic modeling and perfectly skillful initialization and forcings would be unity whereas that with perfect hydrologic modeling but with skill-less initialization and forcings (such as climatology) would be zero. The asymptotic CRPSS of CLIMP at long lead times also reflects the quality of hydrologic modeling. If hydrologic modeling is perfect, the future forcings are skill-less and the ICs are skillful, the CRPSS of CLIMP would be significantly positive at short lead times and decrease toward zero as the lead time increases. Because hydrologic modeling is never perfect in reality, the CRPSS of CLIMP typically turns negative as the lead time increases by varying margins.

The rate of decay of the CRPSS of CLIMP reflects the hydrologic memory of the catchment for streamflow generation captured by hydrologic modeling. The more slowly or quickly the CRPSS of CLIMP decays with increasing lead time, the longer or shorter hydrologic memory is, respectively, in reflection of the residence time of the streamflow-generating runoff. Hence, one may consider the area between the CRPSS of CLIMP (green) and the line of zero CRPSS as representing the skill due to a combination of hydrologic modeling, initialization and hydrologic memory as represented by the hydrologic models. A negative area due to negative CRPSS contributes negatively. The difference in CRPSS between GEFS_CM and CLIMP (i.e., the area between the black and green lines) represents the skill added by the NWP precipitation forecast and its post processing. The difference in CRPSS between DA and GEFS_CM (i.e., the area between the red and black lines) represents the skill added by streamflow assimilation. The difference in CRPSS between PP and DA (i.e., the area between the blue and red lines) represents the skill added by streamflow post processing.

In Figure 9a, the contributing areas associated with CLIMP, GEFS_CM, DA and PP are 0.22, 1.60, 1.18 and 0.75 (days), respectively. Hence, on average over all 43 locations, the 4 components of f₄(s_f| b_f,i_c), $f_6\left(b_f\right|\overline{p_f})$, f₅(i_c| q_o) and $f_7\left(q_f|\overline{s_f}\right)$ in Eq. (4) contribute 6.0, 42.6, 31.5 and 19.9% to the collective predictive skill, respectively, when all ranges of flow are considered. The rather low contribution by f₄(s_f| b_f,i_c) is due to the fact that its CRPSS usually falls below zero at large lead times which contributes negatively. Figs 9bcd shows the results for observed flow exceeding the 90, 95^th and 99^th percentiles, respectively. It is immediately clear that the conditional results of Figs 9bcd are very different from the unconditional results of Figure 9a. For verifying observed flow exceeding the 99^th percentile, streamflow post processing contributes little and the contribution of streamflow DA is significant only over the first two days or so of lead time where the IC memory is strong. For Figure 9d, the contributing areas associated with CLIMP, GEFS_CM, DA and PP are 0.58, 3.32, 0.30 and 0.01 (days), respectively, and the 4 components of f₄(s_f| b_f,i_c), $f_6\left(b_f\right|\overline{p_f})$, f₅(i_c| q_o) and $f_7\left(q_f|\overline{s_f}\right)$ in Eq. (4) contribute on average over all 43 locations 13.8, 78.9, 7.1 and 0.2% to the collective predictive skill, respectively.

The difference between the gray line and the black line (GEFS_CM) represents the gain or loss in CRPSS due to CBPR replacing OLSR in the MEFP. Figure 9 shows that, for all ranges of flow, using OLSR is slightly better than using CBPR. For the 99^th conditional percentile, however, CBPR significantly improves over OLSR with a contributing area of 0.51 (days) which approaches that of CLIMP of 0.58 (days). The above result is due to actively trading smaller Type I error for small precipitation events for greatly reduced Type II error for large events (Kim et al. 2025). Figure 9 indicates that, given the currently available NWP and hydrologic modeling capabilities, the largest positive impact on the accuracy of medium range ensemble streamflow forecasts is likely to come from reducing future input uncertainty by improving post processing of NWP precipitation forecast for large events and reducing hydrologic uncertainty with a combination of streamflow assimilation for shorter lead times and streamflow post processing for longer lead times.

Figure 10 shows the reliability diagrams for PP for various thresholds for all 43 locations for selected lead times. Figure 10 may be compared with Figure 7 which does not reflect future input uncertainty. Note that, in addition to future input uncertainty, Figure 10 also reflects additional sampling uncertainty relative to Figure 7 due to the increased dimensionality in numerical integration. Expectedly, ensemble predictions that reflect both future input and hydrologic uncertainties are not as reliable as those that reflect hydrologic uncertainties only. The deterioration in reliability is particularly significant at longer lead times and at the 99^th percentile. Figure 10 indicates that, for prediction of instantaneous flow, the post-processed DA-aided predictions are reliable in the strong sense only up to the first 4 days or so of lead time and the 95^th percentile threshold for most locations (see also the CRPSS_PP vs. ${\mathit{CRPSS}}_{\mathit{PP}}^{\mathit{POT}}$ in SD 6).

Figure 11 shows the AUC of PP, DA and GEFS_CM vs. lead time for verifying observed flow exceeding the 50^th, 90^th, 95^th and 99^th percentiles, respectively. SD 7 shows the variations across all 43 locations for various conditioning percentiles. Figure 11 and SD 7 may be compared with Figure 8 and SD 5, respectively. Figure 11 indicates the following. Streamflow assimilation consistently improves discrimination over GEFS_CM across all lead times and conditioning thresholds. As in Figure 8, streamflow post processing has little impact on the discriminatory skill of DA. The higher the conditioning threshold is, the more quickly the discriminatory skill decreases as the lead time increases. At the 99^th percentile, the discriminatory skill is almost completely lost by Day 14 for most locations. Collectively, Figs 9 through 11 delineate the limits of predictive skill in the mean sense for the study locations for instantaneous flow and illustrate how such verification information for the full HEFS may be used for user decision making and support.

All results presented above pertain to instantaneous flow. In many applications, prediction of multi-daily flows is of particular interest. To assess the variations of predictive skill with temporal averaging scale, we examined the contributions of the 4 components of the full HEFS for predictions of daily, 3-daily, weekly and biweekly mean flows. Figure 12 shows the selected CRPSS results averaged over the 6 CN- and 4 OHRFC locations for all ranges of verifying observed flow. These two regions are chosen for their contrasting predictability. Figs 12ab show the average CRPSS vs. lead time of CLIMP, GEFS_CM, DA and PP for the CN- and OHRFC locations, respectively. Figs 12cd show the average CRPSS vs. lead time of mean 6-hourly, daily, 3-daily, weekly and biweekly PP for the CN- and OHRFC locations, respectively. The relative contributions of the 4 components do not vary significantly with the temporal averaging scale. Hence, only the instantaneous flow results are shown in Figs 12ab. Figure 12 indicates the following. The collective predictive skill varies greatly from region to region (see also SD 6) and, as expected, increases as the averaging scale increases. For CNRFC, the CRPSS of CLIMP indicates longer hydrologic memory due in part to snowmelt as seen in the slower decay vs. lead time. For OHRFC, the CRPSS of CLIMP indicates shorter hydrologic memory as seen in the faster decay vs. lead time. The contrasting hydrologic predictability between the two regions partly explains why DA and PP provide large improvement for the CNRFC locations but not for the OHRFC locations. For the CNRFC results, the relatively low CRPSS for streamflow analysis for GEFS_CM and the negative CRPSS beyond Day 8 suggest that there is significant room for potential improvement in hydrologic modeling.

All results presented above are expressed in terms of verification metrics. To assess the real-world significance of the differences in skill, it is useful to visualize the results in physical units. Figure 13 shows the scatter- (crosses) and quantile-quantile (QQ, red solid line) plots of ensemble mean forecasts of 0-6, 0-24, 0-72, 0-168 and 0-336 hr mean flow vs. verifying observed mean flow. All plots combine all 43 locations for simplicity albeit being dominated visually by high discharge locations. Because ensemble mean forecasts linearly reflect all ensemble members, they are sensitive to outlying members while preserving mass in the mean sense. In addition, ensemble mean forecasts are perfectly sharp and hence do not hedge. For these reasons, ensemble mean forecasts offer a discerning check on the quality of ensemble forecasts. In Figure 13, the tighter the scatter is around the one-to-one line, the more accurate the ensemble mean prediction is and, the closer the QQ plot is to the one-to-one line, the more similar the marginal distribution of the ensemble mean prediction is to that of the verifying observation. Also shown in each plot are multiplicative bias (BIAS), root mean squared error (RMSE) and correlation coefficient (CORR) where BIAS is the ratio of the sum of the forecast flow to that of the verifying observed flow. Figure 13 shows that DA tightens the scatter significantly compared to GEFS_CM particularly over the lower ranges of observed flow, and that PP aligns the QQ plot more closely to the diagonal than DA while tightening the scatter further albeit modestly as reflected in the RMSE values.

4. Conclusions and future research recommendations

The main findings and conclusions are as follows. It is shown that the 4 scientific components of the full HEFS serve specific and complementary roles toward maximizing the collective predictive skill for medium range (~14 days) ensemble streamflow forecasting as designed at inception (Demargne et al. 2014). When all ranges of observed flow are considered, hydrologic modeling, the use of NWP precipitation forecast and its post processing, streamflow assimilation and post processing of DA-aided streamflow predictions contribute on average over 43 headwater locations 6, 43, 31 and 20% to the collective predictive skill in reference to climatology, respectively. In the above, hydrologic modeling is represented by initializing the operational models using perturbed warm states and forcing them with skill-less future input. Hence, the hydrologic modeling component is equivalent to ESP but with perturbed ICs. The rather low average contribution by hydrologic modeling is due to the negative skill at large lead times vs. climatology at most locations.

The component-specific contributions of skill vary greatly from location to location and with the magnitude of the verifying observed flow in reflection of the widely varying predictability of precipitation and streamflow. When all ranges of observed flow are considered, streamflow assimilation and post processing contribute significantly to the collective skill at shorter and longer lead times for locations where hydrologic memory is strong and streamflow is highly predictable, respectively. When only high flows are considered, the use of NWP precipitation forecast and its post processing is by far the largest contributor. For several locations, the collective predictive skill is noticeably poor due to deficient or very challenging hydrologic modeling particularly where snowmelt is important.

Given the limits of predictability of the input forcings and streamflow, the effectiveness of NWP precipitation forecast and its post processing, streamflow assimilation and streamflow post processing depends by varying degrees on the goodness of hydrologic modeling. For example, improving hydrologic modeling would in turn improve the efficacy of DA and post processing. The results strongly suggest that, for holistic improvement of medium range ensemble streamflow forecasting, efforts should be directed at improving hydrologic modeling particularly at skill-challenged locations in addition to improving future input forecasts particularly for large events. Finally, the full HEFS implemented in this work provides a reference against which competing ensemble streamflow forecasting methodologies for headwater basins may be assessed. Work is under way to comparatively assess artificial neural network-based techniques vs. the full HEFS for diverse locations of predictability and hydroclimatology, including skill-challenged locations, and the results will be reported in the near future.