Improving the Characterization of Initial Condition Improving the Characterization of Initial Condition for Ensemble Streamflow Prediction Using Data for Ensemble Streamflow Prediction Using Data Assimilation Assimilation

. Within the National Weather Service River Fore-cast System, water supply forecasting is performed through Ensemble Streamﬂow Prediction (ESP). ESP relies both on the estimation of initial conditions and historically resam-pled forcing data to produce seasonal volumetric forecasts. In the western US, the accuracy of initial condition estimation is particularly important due to the large quantities of water stored in mountain snowpack. In order to improve the estimation of snow quantities, this study explores the use of ensemble data assimilation. Rather than relying entirely on the model to create single deterministic initial snow water storage, as currently implemented in operational forecasting, this study incorporates SNOTEL data along with model predictions to create an ensemble based probabilistic estimation of snow water storage. This creates a framework to account for initial condition uncertainty in addition to forcing uncertainty. The results presented in this study suggest that data assimilation has the potential to improve ESP for probabilistic volumetric forecasts but is limited by the available observations.


Introduction
Accurate prediction of seasonal streamflow information is essential to effectively managing surface water supply.For this reason, recent studies have examined techniques that have potential to provide skillful predictions of seasonal runoff volume (Kennedy et al., 2009;Moradkhani and Meier, 2010;Regonda et al., 2006;Thirel et al., 2008).It has been well documented that the seasonal volume of runoff is controlled by both the initial water storage of the land surface at the beginning of the season and the future water fluxes into Correspondence to: H. Moradkhani (hamidm@cecs.pdx.edu) and out of the system (Lorenz, 1975).The effect of the initial water storage is particularly important in mountainous regions where melt from seasonal snowpack can dominate the spring and summer streamflows (Pagano et al., 2004).Given that there is accurate quantification of the snow water storage in a specified region, this information can be utilized to improve the accuracy of seasonal streamflow prediction.This idea is implemented in the National Weather Service River Forecast System (NWSRFS) through Ensemble Streamflow Prediction (ESP).
The ESP framework was first introduced by Twedt et al. (1977) and later clarified by Day (1985).In ESP, the initial condition of the water stored in the snowpack and soil are produced by running hydrologic models up to the initial forecast time-step.This provides information to initialize the model for seasonal forecasting with an ensemble of historical forcing data.Since the approach for generating the initial condition does not account for uncertainties in the modeling framework, and the sampled forcing data is not necessarily representative of the future climate, it is advantageous to both constrain the forcing data and account for the uncertainties associated with the initial conditions.Previous work has focused on constraining the forcing to more realistic predictions (Najafi et al., 2011), reducing the error in initial conditions (McGuire et al., 2006;Tang and Lettenmaier, 2010) and examining the effect of uncertainties in the initial conditions (Wood and Lettenmaier, 2008;Li et al., 2009).While both initial condition and forcing uncertainty are important to address, this study focuses on using ensemble based data assimilation to more accurately quantify the uncertainty with respect to snow.In this study, it is hypothesized that a more accurate quantification of uncertainty in model initial conditions will lead to a more reliable ESP forecast.
Modeling of snow accumulation and ablation is subject to a range of uncertainties.These uncertainties stem from errors in observing the forcing data, model structure and parameterization.In addition to modeling errors, observing the initial C. M. DeChant and H. Moradkhani: Ensemble streamflow prediction using data assimilation condition of the snowpack is complicated because of the spatially heterogeneous and complex nature of snowpack.In order to improve estimation of snowpack states, several recent studies have examined combining the modeled and observed states through data assimilation to increase the accuracy of snow estimates and quantify their uncertainty (Andreadis and Lettenmaier, 2006;DeChant and Moradkhani, 2011a;Durand et al., 2009;Leisenring and Moradkhani, 2010;Rodell and Houser, 2004;Slater and Cark, 2006;Sun et al., 2004;Zaitchik and Rodell, 2009).These studies showed that there is potential for using snow data assimilation to improve snow estimates through a variety of observations including in-situ (SNOTEL) and remotely sensed (snow cover, snow water equivalent and passive microwave brightness temperature) measurements.For the purposes of this study, SNOTEL observations are utilized for their simplicity and reliability.While remotely sensed data have the potential to be more effective in quantifying spatial quantities of snow, the use of these observations for data assimilation is still in the development phase and therefore not fit for this study.
This study is organized into 4 subsequent sections.The methods section discusses the SNOW-17 and Sacramento Soil Moisture Accounting (SAC-SMA) models, the study area and the SNOTEL observation data.The third section describes the experimental design, including a description of data assimilation, ESP, with and without data assimilation, and the performance metrics with which the results are verified.This is followed by a results section and the final section contains a brief discussion and conclusion.

SNOW-17
The SNOW-17 model is used operationally within the NWS-RFS and is the snow model used in operational ESP forecasts.SNOW-17 is a temperature index model that estimates simplified vertical snow processes (Anderson, 1973).The main processes simulated by SNOW-17 include: form of precipitation (snow or rain), accumulation of snow cover, energy exchange at the snow-air interface, internal states of snow cover (temperature, liquid/frozen water content, density, etc.), transmission of liquid water through the snowpack, and heat transfer at the soil-air interface.With six-hourly inputs of precipitation and air temperature, the model predicts the amount of snow accumulation and melt that occur.In order to account for spatial and elevation heterogeneities of snow, the model is run for two or three separate elevation bands for each sub-basin.Historical forcing data and model parameters for each basin elevation band were provided by the Colorado River Basin River Forecast Center (CBRFC).In the data assimilation portion of this study, the precipitation is assumed to have a log-normal error standard deviation equal to 25 % of the precipitation value and the temperature is assumed to have a normal error with a standard deviation of 3 degrees Celsius.This is based on errors suggested by DeChant and Moradkhani (2011a), which used similar forcing data.

Sacramento soil moisture accounting model
The SAC-SMA model, first introduced by Burnash et al. (1973), is the model used operationally at the NWSRFS to translate snowmelt and rain values into streamflow.The model simulates water storage with two soil moisture zones: an upper and a lower zone.The upper zone accounts for short term storage of water in the soil, while the lower zone models the longer term groundwater storage.Water can move vertically from the upper zone to the lower zone, laterally out of the system depending on the state variables and the parameterization, or vertically out of the system through evapotranspiration.The SAC-SMA model is run with snowmelt information from the SNOW-17 model and the potential evapotranspiration (PET), is linearly interpolated from the monthly PET values for each elevation band provided by the CBRFC for the study basins.The model calculates the water balance for the system and any excess is routed to the basin outlet using the unit hydrograph method.

Study area
This study takes place in the Upper Colorado River basin.Fifteen separate sub-basins were analyzed to determine an average effect that data assimilation has on ESP.These fifteen basins are summarized in Table 1.The locations of the basins within the Upper Colorado River Basin are shown in Fig. 1.These 15 basins were chosen because of the availability and reliability of forcing data, streamflow data and nearby SNOTEL observations.

SNOTEL
SNOTEL sites are managed by the National Resources Conservation Service (NRCS) and provide in-situ observations of snow depth, snow water equivalent (SWE), precipitation, and temperature.Some enhanced sensors can measure more variables, such as soil moisture.The observation used in the data assimilation portion of this study is SWE.At SNOTEL sites, SWE is measured by a snow pillow.A snow pillow is a pressure sensitive pad that weighs the snowpack, which can be directly translated into the volume of water that would be released if the snowpack was melted.In addition, the snow depth is measured with a sonic sensor, the precipitation is measured with a storage type gage and air temperature is measured with a shielded thermistor (NRCS, http://www.wcc.nrcs.usda.gov/snotel/SNOTEL-brochure.pdf).Each SNOTEL site is chosen as the observation of a given model elevation band based on horizontal (latitude and longitude) and vertical (elevation) Hydrol.Earth Syst.Sci., 15, 3399-3410, 2011 www.hydrol-earth-syst-sci.net/15/3399/2011/    evidence that the middle elevation bands are well represented in terms of elevation in comparison to the low bias of observations for the highest elevation bands and the high bias of the lowest elevation bands.Similar elevations between the model band and the SNOTEL observation are necessary for an accurate update because the timing of peak snowpack accumulation changes dramatically with elevation.As will be described in the results, the elevation of the band also controls the timing of snowmelt, which strongly affects the modeling results.In order to account for the uncertainty in horizontal representativeness of the SNOTEL stations, which is strongly affected by the spatial distribution of forcing, vegetation and wind, a temporally and spatially uncorrelated error, with a standard deviation of 25 % of the observation, is considered for the SNOTEL SWE.This value was found to be adequate by examining the accuracy of the expected value and the reliability of probabilistic streamflow estimation during model hindcasts, using Nash-Sutcliffe Efficiency, Rank Probability Skill Score and Normalized Root Mean Square Error for different error values.Of the error values explored, 25 % was found to produce the best verification scores on average, across all 15 study basins.In addition, uncorrelated errors were assumed for simplicity but the effects of correlated errors will need to be evaluated in future studies.
3 Experimental design

Data assimilation using the particle filter
The term data assimilation refers to a variety of techniques aimed at combining model simulation and observations to account for uncertainties in both state reconstruction techniques.In this study, ensemble based techniques are employed because they directly address the prediction of uncertainty in the desired state.Of the available ensemble data assimilation techniques, the Ensemble Kalman Filter (EnKF) is the most popular in the hydrologic literature.Several stud- ies in the past decade have applied this technique to hydrologic models (Andreadis and Lettenmaier, 2006;Durand et al., 2009;Moradkhani et al., 2005a;Reichle et al., 2002;Roddell and Houser, 2004;Slater and Cark, 2006;Sun et al., 2004;Zaitchik and Rodell, 2009).Though these studies have shown that the EnKF is a valuable tool for data assimilation in many applications, this study focuses on the use of the Particle Filter (PF).According to recent studies, the PF is an effective hydrologic and hydraulics data assimilation method, providing predictive uncertainty in model states, parameters and fluxes (DeChant and Moradkhani, 2011a, b;Matgen et al., 2010;Leisenring and Moradkhani, 2010;Moradkhani et al., 2005b;Montzka et al., 2011;Rings et al., 2010;Weerts and El Serafy, 2006).The PF is not subject to the limitations experienced in the EnKF including the Gaussian assumption of joint distribution of observation and model states and the linear updating of model states (Moradkhani et al., 2005b;Moradkhani and Sorooshian, 2008).In addition, the PF is capable of preserving the dynamical balance in the analysis (i.e. the conservation of mass) because the ensemble members are only resampled as opposed to being linearly adjusted, as in the EnKF (Moradkhani, 2008).Maintaining a water balance is assumed to beneficial in this framework, due to the water balance assumption in the SNOW-17, but the value of maintaining a water balance in a data assimilation framework has not been entirely quantified.For these reasons, the PF is used to assimilate the snow observations into the hydrologic models.
In order to explain data assimilation, the modeling framework must be viewed in the state-space framework.The model is SNOW-17, which estimates the prior snowpack state based on the posterior snowpack states at the previous timestep and forcing data.As the model progresses forward in time, the prior distribution of states is produced according to Eq. (1).
where f is the forward operator (hydrologic model), represents the model predicted (prior) states (e.g.density, depth, liquid water in snowpack), x + i,t−1 represents the updated model states at the previous time-step, u i,t represents the meteorological forcing data (precipitation and temperature), ω i,t is the model error, which is normally distributed with a standard deviation of 25 % of x − i,t , i is the ensemble member and t is the time-step.This model error was determined in conjunction with SWE observation error as described in Sect.2.3.Prior to update of the model states, an observational operator must be applied to transfer the states into the observation space, as in Eq. ( 2).
where y i,t is the prediction (SWE) and ν i,t is the observation error.
Based on the recursive Bayes Law (3), the PF sequentially samples prior states and parameters to create an accurate posterior distribution, at each observation time-step.
Equation ( 3) shows mathematically that a posterior conditional probability distribution of model predicted states (x t ), given all previous observations (y 1:t−1 ) and the current observation (y t ), can be computed sequentially in time.The observation is the SNOTEL observed SWE.In this study, the probability of each particle, p(y t |x t ), is calculated via the normal likelihood Eq. (4).
The normalized likelihood, p y t |x i,t , can easily be calculated by: In the above equations, R k is the variance of the observations, which is 25 % as described in Sect.2.3.This probability is necessary to transform the prior particle weights into the posterior via Eq.( 6).
In the PF with resampling, prior particle weights, w − i,t , are set equal to 1/N p (N p is the number of particles) before moving on to the next time-step.This results in a posterior weight, w + i,t , equal to p y t |x i,t which is the normalized likelihood.In this study, we rely on Sequential Importance Resampling (SIR) method as elaborated in Moradkhani et al. (2005b) and Moradkhani (2008).In the data assimilation portion of this study, 500 ensemble members, or "particles", are used for snow estimation.

Ensemble streamflow prediction (ESP) with and without data assimilation
ESP is a method used by the NWSRFS to create probabilistic forecasts of seasonal streamflow volumes.This provides a prediction of the uncertainty in seasonal streamflow resulting from unknown forcing data and initial model states at the time of forecast, as shown in Fig. 4. Initial model states are produced by running the model with observed input data up to the initial forecast time-step.This is called a "spin-up".
Starting at this point, the model is forced with resampled historical forcing, beginning at the initial forecast date, for each historical observation year.This produces a potential streamflow trace that could occur given the current state of the land surface and a previously observed forcing dataset.Given that the assumptions of seasonal forcing stationarity and accurate model initial states are not violated, ESP can provide a skillful probabilistic prediction of seasonal streamflow.
ESP can be coupled with data assimilation for state initialization, also shown in Fig. 4. Rather than beginning the ESP forecast with a spin-up, as is done in traditional ESP, ESP-DA implements a data assimilation strategy to characterize the initial states.This creates an ensemble of state values that represent a probability density (Fig. 4b).After a sufficient number of time-steps, and assuming the uncertainty with respect to the model and observation are accurately quantified, the state distribution produced by data assimilation will accurately reflect the uncertainty in the state with respect to the model and observation.Similar to the spin-up, initialization through data assimilation requires many time-steps to reach reliable state values.At the forecast date, a given number of ensemble members are sampled from the state PDF, 50 in this study, and the ESP is performed from each of these ensemble members (shown in Fig. 4b with 8 sampled ensemble members for visibility).This propagates the uncertainty from the initial condition through the ensemble forecast, which is hypothesized to create a more accurate estimation of uncertainty in streamflow prediction, as is represented in the forecast PDF in Fig. 4b.
A sample of 50 ensemble members from the state PDF was motivated based on the need to minimize the computation time of the algorithm.Since this computational demand increases dramatically with larger sample sizes, this study attempted to minimize the number of ensemble members necessary for ESP-DA.50 ensemble members was found to be suitable by comparing a histogram of the sampled states to the full state distribution, and by analysis of seasonal streamflow prediction results in a small number of the study basins.Based on the reliability of ESP generated by sampling 50 ensemble members to different sample sizes, 50 ensemble members was found to balance the increase in computational demand, due to large sample size, and loss in reliability of the forecast, due to small sample size.

Performance metrics
Ranked Probability Score (RPS) is a widely used measure for evaluating the quality of probabilistic predictions Wilks (2006).RPS is the sum of squared error of the cumulative probability forecasts averaged over multiple events.In streamflow prediction, the probability forecast is usually expressed using a non-exceedance probability forecast within pre-specified categories (i.e. 1 %, 25 %, 50 %, 75 % and 99 % non-exceedance).The observed value for a given threshold (forecast category) takes on the value of 1 if the observed flow value is less than the threshold for that category.Otherwise, the observed value is 0. The discrete expression of RPS is given as: where F t j is the forecast probability at time t given by P (forecast j < thresh j ) and O t j is the observed probability given by P (observed < thresh j ) where j is the probability category.The Rank Probability Skill Score (RPSS) is also computed as the percentage improvement over a reference score (e.g.climatology) Wilks (2006): where RPS climatology is the rank probability score for the observation.A positive value shows to the percentage of improvement over the reference RPS.
In the analysis of total streamflow volume, a rank histogram and Q-Q plot are used to analyze the accuracy of the uncertainty prediction.In the Rank Histogram, the rank in which the observation falls on each ensemble prediction is presented.The rank is calculated according to the following equation.
Rank n = I y n,i < y n and y n,i+1 > y n (9 where y is the seasonal volumetric prediction and y is the observation, similar to the PF explanation.In addition, n is the index for each separate seasonal forecast.All ranks are then placed into a histogram.Since the ensemble volume is assumed to make up a probability density, a uniform histogram, in which the observation falls equally as often in each rank, indicates accurate representation of uncertainty.For a detailed description of Rank Histogram interpretation, see Hamill (2001) andTalagrand et al. (1997).Similar to the rank histogram, a Q-Q plot provides information about the accuracy of the uncertainty estimation of the ensemble forecast.A Q-Q plot is created by first calculating the normalized rank (z) of each observation z n = I y n,i < y n and y n,i+1 > y n N (10) where N is the total number of seasonal forecasts, which is 135 in this study.The ranks are then sorted and graphed.This graph is compared with a uniform line.A Q-Q plot matching the uniform line indicates a highly reliable forecast ensemble prediction.For details of the interpretation of a Q-Q plot, see Laio and Tamea (2007).While this plot is very useful in determining the effectiveness of the ensemble prediction, supportive quantitative scores are also included with this plot.
Both the reliability and resolution, as described by Renard et al. (2010), can be adapted for this study to analyze the accuracy of probabilistic seasonal volume prediction.Equations (12,13 and 14) describe the two reliability statistics and Eq. ( 15) describes the resolution statistic.
With respect to the reliability of the forecast, α is a measure of the uniformity of the Q-Q plot (worst value is 0 and perfect score is 1) and ε is a measure of the portion of observations inside the ensemble prediction (a value closer to 1 indicates a greater occurrence of observations falling within the ensemble).π is a measure of the resolution of the forecast, which will be referred to as precision.Given that two forecasts have similar reliability, the forecast with greater precision is preferred because it indicates a lower level of uncertainty.In equation 12 through 15, U is the normalized uniform cumulative density.In equation 15, E[] is an expected value operator and SD[] is a standard deviation operator.

Results
Prior to analyzing the streamflow prediction results in this study, it is necessary to compare the initial conditions created by the spin-up and data assimilation.This comparison is made in Fig. 5.In this figure, the initial states for the forecasts beginning on 1 March, 1 April, 1 May and 1 June for 2003 through 2005 are presented.Each sub-plot contains the total water stored (TWS) as snow (sum of all 15 basins) for the spin-up, shown as a single value, and the data assimilation, shown as a distribution of values, which represent Poor representation of the upper and lower elevation bands by SNOTEL has important consequences on the accuracy of data assimilation.Since the observation is most representative of the middle elevation band, the lowest elevation bands will be forced to peak later in the model than in reality, and the highest elevations will be forced to peak earlier in the model than in reality.Overall this causes a trend of miscalculation of SWE in the winter months, when the majority of the flow would be from the lower elevation band, and late spring/summer months, when most of the flow would results from high elevation melt.Since SNOTEL is less representative in the lower and upper elevation bands than the middle elevation band, it may be possible to improve the results by using different error characteristics for these bands.Though this would be a useful contribution to SNOTEL data assimilation in the Upper Colorado River Basin, this study does not attempt to manage these errors due to the complexity of melt timing differences between SNOTEL and the model.Because of the errors early and late in the snow accumulation/ablation season, the ESP-DA results are only expected to improve seasonal forecasts beginning on 1 March, 1 April and 1 May.
In addition to errors relating to the elevation of the observations, the model and observations differ on the quantity of snow accumulation for each year.Figure 5 suggests that the model tends to produce similar snow quantities in each year while the SNOTEL stations have their greatest accumulation in 2003.This leads to a similar prediction of snow quantities in 2003 but a strongly low bias in the following two years.Though elevation differences between the model and SNO-TEL are expected to cause errors in the streamflow prediction, yearly differences are expected to improve the seasonal prediction.With an understanding of the effects that data assimilation has on model states during this time, the ensemble streamflow forecasts can be effectively analyzed.
Three month forecasts beginning in March and April are presented in Fig. 6 and beginning in May and June are presented in Fig. 7.These figures show the cumulative runoff for each day in the forecast period, summed across all 15 study basins.Beginning in March, there is very little cumulative runoff through the first month.This is because the snow accumulation will continue through most of March in a typical year.In April, an increasing amount of flow is observed, but it is not until May when the most rapid increase is observed.During May, most of the area within the 15 basins is experiencing significant snowmelt.These high flows do not taper off until late June.By late June, most of the snowpack has melted and only the highest elevations will have snow to melt.For most forecast dates, the uncertainty added through the initial states by data assimilation would be expected to increase the uncertainty in volumetric runoff.Rather than increasing the uncertainty, these figures suggest that the uncertainty is constrained by ESP-DA in comparison to ESP, through most forecasts.This is due to a generally low bias in the ESP-DA SWE in comparison to the ESP SWE.Since there is a lower boundary condition for the flow (baseflow), and no upper boundary condition, the low bias in the snow states has caused most of the ESP-DA forecasts to have less uncertainty in runoff volume than traditional ESP.During the 1 May and 1 June seasonal predictions for 2003, a larger uncertainty in initial conditions, with only a slight bias, is found to increase the seasonal volumetric prediction uncertainty.With respect to the 1 June prediction, it is important to note that, though the upper 95 % predictive bound of the ESP-DA is lower than that of the traditional ESP, the maximum value of the ESP-DA is actually greater than that of the traditional ESP.This is shown by Fig. 8, which shows the total 3 month volumetric streamflow forecast from all 15 basins for each forecast date.Overall for the forecasts beginning on 1 March, 1 April and 1 May, the biased predictions from ESP-DA, comparing to traditional ESP, appear to have more accurately bounded the observation.This is shown to improve the probabilistic prediction with an average improvement in RPSS of 7.5 in 2004 and 11.5 in 2005.An improvement in RPSS suggests that the cumulative ensemble time series of prediction from ESP-DA is more accurate than traditional ESP.Though an improvement was found in forecasts    prediction with a slight low bias.In addition, the higher α for ESP-DA, in relation to ESP, indicates a closer to uniform Q-Q plot, while a higher ε indicates more observation lie within the ensemble range.This implies a generally more reliable prediction.It is also important to note that the ESP produced a more precise prediction (as indicated by a larger π value) than ESP-DA.This is expected because the ESP technique accounts for fewer sources of uncertainty than ESP-DA.Overall the results suggest that seasonal predictions beginning in March, April and May in the upper Colorado River Basin more accurately characterize the uncertainty when initialized by data assimilation than with a model spin-up.

Discussion and conclusion
This study examined the utility of incorporating ensemble data assimilation techniques to improve state initialization in the ESP framework by allowing for uncertainty in the initial condition.A combined ESP-DA framework was implemented in 15 sub-basins in the upper Colorado River Basin.This was compared against traditional ESP to determine if an improved representation of uncertainty could be achieved through data assimilation.Though positive results were found for ESP-DA, issues were found relating to the snow data assimilation.
In general, the flaws in this study stem from the lack of a representative observation for all basins.Since the SNOTEL stations tend to be in the range of middle elevation bands, the upper and lower elevation bands are known to be inaccurately adjusted.There is potential to improve the snow data assimilation in this study, either through more effective management of errors relating to elevation misrepresentation, use of more realistic error assumptions (e.g.accounting for spatial and temporal error correlations), or use of more representative datasets, but rigorous study is still needed in the field of snow data assimilation to achieve this.Though the assim-  ilation process is known to be flawed, some positive results can still be observed from this study.In the late spring and early summer, when the runoff is dominated by the middle elevation bands, the results presented here suggest that ESP can effectively be initialized through SNOTEL data assimilation.Furthermore, initialization through data assimilation can improve the ability to estimate seasonal runoff volume uncertainty.With this result, it can be inferred that as the understanding of snow data assimilation improves, ESP-DA will become more effective for seasonal streamflow prediction.This highlights the potential for accurate seasonal forecasts in mountainous regions but emphasizes the need for improved snow estimation techniques to achieve a more accurate forecast.
In addition, the value of ensemble data assimilation within an operational framework is often debated.Ensemble data assimilation may become questionable in an operational framework due to the computational complexity.The NWS hydrologists currently implement run time modifications to adjust the state in model manually to improve forecasts.Manual adjustments require less computing time and it has yet to be proven that automatic data assimilation is more effective in state estimation than an experienced forecaster.One clear benefit of automatic data assimilation comes from an accurate accounting of initial condition uncertainty.Since data assimilation can characterize an entire initial condition PDF, not just initial condition adjustments, it is likely that automated data assimilation will provide a more accurate characterization of uncertainty than manual modifications.In addition, the computational demand may be a large obstacle for real-time forecasting at the NWS due to the daily time constraints of the forecasts, but seasonal forecasting relaxes the time constraints on producing the forecast, therefore making ensemble data assimilation more attractive.For this reason, it is assumed in this paper that data assimilation is a viable technique for seasonal streamflow forecasting.
Overall this study found that ESP-DA has potential to improve characterization of uncertainty over ESP in snowmelt dominated basins.As expected, improvements were found during the period that achieved accurate state initialization.The results presented here show the importance of both accurate state initialization and accurate estimation of the uncertainty of the initial states.By more accurately characterizing the uncertainty in the states, the total seasonal flow uncertainty is more accurately represented.In future studies, the ESP-DA techniques should be tested with new methods for assimilating snow information and in different climatic regions to confirm the effectiveness of ESP-DA for different types of climates.ESP-DA can also be easily coupled with advanced techniques to constrain the forcing uncertainty.Therefore state initialization of ESP using data assimilation should be considered a potential tool for improving ESP forecasts in snow dominated basins.

Figure 1 .Fig. 1 .
Figure 1.Map of the Upper Colorado River Basin (CRB) with highlighted study basins and SNOTEL locations proximity.In order to show the representativeness of SNO-TEL stations in relation to the model elevation bands, Figs. 2 and 3 are presented.Figure2shows the vertical proximity of model elevation bands and their respective SNOTEL observation site.It is important to note that model elevation bands above 2800 m and below 3400 m have relatively sim-ilar elevations to the closest SNOTEL observation.Depending on the basin, the upper, middle or lower elevation band may fall within this range.Due to the inconsistency in model elevations and number of bands, all elevation bands between 2800 m and 3400 m will be referred to as the middle elevation bands to simplify future discussion.Figure3provides further

Fig. 2 .
Fig. 2. Scatter-plot comparing the elevation that each model band represents and the elevation of its respective SNOTEL observation site.

29Fig. 3 .
Fig. 3. Mean absolute difference (left) and mean difference (right) between model elevation and its respective SNOTEL elevation.A positive mean difference indicates the model is at a higher elevation than its SNOTEL station.

Fig. 4 .
Fig. 4. Diagram of the ESP and ESP-DA algorithms.(a) shows the uncertainty bounds of the data assimilation and forecast, (b) shows the individual traces generated from each sample.

Fig. 5 .
Fig. 5. Comparison of the total water stored (TWS) as snow produced by the spin-up (SU) and data assimilation (DA) for the four forecast dates summed across all 15 basins.

Fig. 6 .
Fig. 6.Cumulative daily volumetric flow plots summed across all 15 basins.This figure shows the 95 % predective bounds of the cumulative runoff volume from ESP (black lines) and ESP-DA (green lines).The expected value of each is a dashed line and the cumulative observed runoff volume is the red line. 33

Fig. 9 .
Fig. 9. Rank histogram of the seasonal volumetric flow prediction for all 15 basins beginning in March, April and May during 2003, 2004 and 2005. 36

Fig. 10 .
Fig. 10.Q-Q plot of the seasonal volumetric flow predictions for all 15 basins beginning in March, April and May during the years 2003, 2004 and 2005.

Table 1 .
Basin data for all of the study basins.