Combined assimilation of streamflow and snow water equivalent for mid-term ensemble streamflow forecasts in snow-dominated regions

The potential of data assimilation for hydrologic predictions has been demonstrated in many research studies. Watersheds over which multiple observation types are available can potentially further benefit from data assimilation by having multiple updated states from which hydrologic predictions can be generated. However, the magnitude and time span of the impact of the assimilation of an observation varies according not only to its type, but also to the variables included in the state vector. This study examines the impact of multivariate synthetic data assimilation using the ensemble Kalman filter (EnKF) into the spatially distributed hydrologic model CEQUEAU for the mountainous Nechako River located in British Columbia, Canada. Synthetic data include daily snow cover area (SCA), daily measurements of snow water equivalent (SWE) at three different locations and daily streamflow data at the watershed outlet. Results show a large variability of the continuous rank probability skill score over a wide range of prediction horizons (days to weeks) depending on the state vector configuration and the type of observations assimilated. Overall, the variables most closely linearly linked to the observations are the ones worth considering adding to the state vector due to the limitations imposed by the EnKF. The performance of the assimilation of basin-wide SCA, which does not have a decent proxy among potential state variables, does not surpass the open loop for any of the simulated variables. However, the assimilation of streamflow offers major improvements steadily throughout the year, but mainly over the short-term (up to 5 days) forecast horizons, while the impact of the assimilation of SWE gains more importance during the snowmelt period over the mid-term (up to 50 days) forecast horizon compared with open loop. The combined assimilation of streamflow and SWE performs better than their individual counterparts, offering improvements over all forecast horizons considered and throughout the whole year, including the critical period of snowmelt. This highlights the potential benefit of using multivariate data assimilation for streamflow predictions in snow-dominated regions.

Abstract.The potential of data assimilation for hydrologic predictions has been demonstrated in many research studies.Watersheds over which multiple observation types are available can potentially further benefit from data assimilation by having multiple updated states from which hydrologic predictions can be generated.However, the magnitude and time span of the impact of the assimilation of an observation varies according not only to its type, but also to the variables included in the state vector.This study examines the impact of multivariate synthetic data assimilation using the ensemble Kalman filter (EnKF) into the spatially distributed hydrologic model CEQUEAU for the mountainous Nechako River located in British Columbia, Canada.Synthetic data include daily snow cover area (SCA), daily measurements of snow water equivalent (SWE) at three different locations and daily streamflow data at the watershed outlet.Results show a large variability of the continuous rank probability skill score over a wide range of prediction horizons (days to weeks) depending on the state vector configuration and the type of observations assimilated.Overall, the variables most closely linearly linked to the observations are the ones worth considering adding to the state vector due to the limitations imposed by the EnKF.The performance of the assimilation of basin-wide SCA, which does not have a decent proxy among potential state variables, does not surpass the open loop for any of the simulated variables.However, the assimilation of streamflow offers major improvements steadily throughout the year, but mainly over the short-term (up to 5 days) forecast horizons, while the impact of the assimilation of SWE gains more importance during the snowmelt period over the mid-term (up to 50 days) forecast horizon compared with open loop.The combined assimilation of streamflow and SWE performs better than their individual counter-parts, offering improvements over all forecast horizons considered and throughout the whole year, including the critical period of snowmelt.This highlights the potential benefit of using multivariate data assimilation for streamflow predictions in snow-dominated regions.

Introduction
Water resource management for reservoirs located in snowdominated regions relies on an accurate portrayal of the snow water equivalent (SWE) spatial and temporal distribution in order to make accurate streamflow predictions.Some water resources managers make use of ensemble streamflow prediction (ESP) to plan reservoir operations over various lengths of time.ESPs have the benefit of integrating weather forecast uncertainty, either by making use of weather ensemble predictions (de Roo et al., 2003) or by using historical weather data (Day, 1985) as input in a hydrologic model.However, ESPs depend heavily on the model's initial conditions (Franz et al., 2008).Presently, many water resources managers still use a manual approach to adjust the initial state of the watershed based on available observations and the user's experience (Liu et al., 2012).
Data assimilation (DA) methods, such as the ensemble Kalman filter (EnKF; Evensen, 2003) can improve the estimation of the initial state of the watershed while also providing an uncertainty on this initial state (Liu and Gupta, 2007).Several authors have already shown the added value of DA in snow-dominated watersheds to improve the estimation of the state of the watershed (De Lannoy et al., 2012;Dechant and Moradkhani, 2011;Nagler et al., 2008;Slater and Clark, 2006;Andreadis and Lettenmaier, 2006).Some studies have also integrated DA in ensemble forecast systems for relatively short-term (up to 5-10 days) hydrologic forecasts (Abaza et al., 2014(Abaza et al., , 2015;;He et al., 2012), but studies focusing on longer forecast periods are scarce even though the need exists for water resource managers.
Multivariate DA applications in hydrology are becoming more frequent, but generally focus on streamflow and soil moisture (Samuel et al., 2014;Trudel et al., 2014;Lee et al., 2011), omitting the snow water equivalent.In snowdominated watersheds, the key initial states include not only information about the hydric state, such as soil moisture and streamflow, but also the snow cover state, such as SWE and snow cover area (SCA).To the authors' knowledge, no published studies pertain to the combined assimilation of information about a watershed's hydric and snow state.Since the lasting impact of hydric DA and snow DA can be quite different given the different physical processes driving them, the simultaneous DA of both types of data could yield improvements over a potentially longer length of time.
However, data assimilation performance depends on various factors, such as the choice of variables to be updated by an observation (hereby referred to as the state vector configuration).Abaza et al. (2015) demonstrated this importance when assimilating streamflow in a hydrologic model.Going from univariate to multivariate DA increases the number of degrees of freedom, which increases the complexity of the matter.The importance of state vector configuration when using multivariate DA for hydrological modelling has yet to be investigated.
The study's main objectives are to (1) investigate the potential impact that multivariate data assimilation of hydric (streamflow) and snow-related (SWE and SCA) data can have on short-term (1-5 days) and mid-term (up to 50 days) streamflow forecast, and (2) to explore how this impact varies as a function of the state vector configuration.

Study area description and data
Simulations were conducted in a synthetic setting based on the Nechako watershed located in British Columbia, Canada (Fig. 1).The watershed includes a reservoir, which drains an area of approximately 14 000 km 2 .The reservoir is managed by Rio Tinto mainly for hydroelectricity production purposes.The watershed includes part of the Coast Mountains in the west region, such that the difference in elevation between the highest and lowest point in the watershed reaches about 1700 m.At this latitude and altitude, most (estimated at 53 %) of the precipitation falls as snow.
There are various types of data gathered regularly over the watershed.First are the seven weather stations managed by Rio Tinto, three of which measure daily precipitation and air temperature only (yellow circles).Three others also in-Figure 1.The Nechako watershed and the locations of weather stations, snow pillows and a hydrometric station.All of these contain at least daily weather data.The outlet is considered to be at the spillway, located at the blue triangle.The intake is located at the Tahtsa Intake weather station (westernmost yellow circle).
clude snow pillows (red squares), which measure the snow water equivalent.The northernmost snow pillow is located at Mount Wells, the southernmost at Mount Pondosy and the westernmost at Tahtsa Lake.Maximum seasonal SWE observations average 615, 853 and 1393 mm for the Mount Wells, Mount Pondosy and Tahtsa Lake snow pillows respectively.The distribution of snow on the ground follows a strong east-west gradient such that measurements at Tahtsa Lake typically yield much more snow than Mount Well and Mount Pondosy.The northernmost weather station (blue triangle) is located next to the spillway at Skins Lake and also takes hydrometric measurements.Historical daily water levels can then be converted into natural inflows by also taking into account spilled and turbined flow.Finally, daily SCA data derived from the spaceborne sensor MODIS/Terra are also considered (Hall et al., 2002).Because of its spatial coverage and relatively high temporal resolution, remotely sensed snow data from MODIS have proven to be valuable in a number of hydrologic studies (Bergeron et al., 2014;Roy et al., 2010;Tang and Lettenmaier, 2010;Andreadis and Lettenmaier, 2006;Clark et al., 2006), including one applied to the Nechako watershed (Marcil et al., 2016).
The meteorological observations gathered over a period of 10 years (from 15 August 1990 to 14 August 2000) were used as a basis upon which a synthetic experiment (see below) tested the added value of three types of synthetic observations (streamflow, SWE and SCA) for data assimilation purposes.The only data sets actually used were the meteorological station data, as well as some streamflow observations (owned by Rio Tinto) and MODIS/Terra daily L3 snow cover data (Hall et al., 2002) for the initial model calibration performed prior to the synthetic experiment presented in this manuscript.All the observation data were created synthetically to mimic streamflow, SWE and SCA data that could be measured or estimated using hydrometric, snow pillow and MODIS data, respectively.More details on the cre-  ation of synthetic observations and meteorological input in Sect.3.1.1and 3.1.2respectively.

Model description
The hydrologic model used was the spatially distributed, conceptual model CEQUEAU (Charbonneau et al., 1977).It is currently being used by Rio Tinto to model hydrologic processes including streamflow at the outlet of the Nechako watershed, considered to be the spillway where the hydrometric station is also located.All variables are computed at a daily time step using a set of parameters to calibrate and daily meteorological input consisting of mean air temperature and precipitation.The set of parameters used in this study was the result of a manual calibration performed by Rio Tinto by comparing the simulated streamflow at the outlet with the corresponding real streamflow observations.A summary of the main processes concerning the production and transfer functions is presented here to facilitate the understanding of the state variables used in this study.
CEQUEAU divides the watershed into regular square pixels called "whole squares" over which the production func-tion is computed (Fig. 2).The current version of CEQUEAU uses the snow model presented by the US Army Corps of Engineers (1956) to simulate most snow-related processes.The SWE is actually computed separately for forested and open areas, which have their own set of parameters, but is aggregated here as a weighted sum according to the proportion of forested and open areas within each whole squares.The only variable computed separately (i.e.outside from CEQUEAU) is SCA, which is computed using a depletion curve (Anderson, 1973).The depletion curve used here follows Andreadis and Lettenmaier (2006), which uses a three parameter beta distribution: where SCA i is the resulting snow cover area over a whole square i, SWE i is the simulated snow water equivalent over the same area, SWE max, i is the annual maximum snow water equivalent from the beginning of the accumulation period over the same area, SI represents the value of SWE above which it is assumed there is always 100 % snow cover and α SCA and β SCA are shape parameters for the beta distribution itself.The calibration of those three parameters was conducted using the SCE-UA method (Duan et al., 1992) to minimize the root mean square difference between simulated snow cover area and MODIS/Terra daily L3 snow cover data (Hall et al., 2002) averaged within each whole square within the Nechako watershed.It is important to note that SCA is computed as an output only and is therefore not considered to be a state variable since it has no impact on future simulations if its value is tampered with.CEQUEAU then uses three conceptual reservoirs to simulate various hydrologic processes from the available water resulting from rain or snowmelt.There is an optional lake reservoir, an upper reservoir (called "soil moisture reservoir" in this study) and a lower reservoir (called "groundwater reservoir" in this study).
All in all, the state variables simulated over each whole square include SWE, a snow ripening index (SRI), a snow temperature index (STI), the soil moisture level (SML), the groundwater level (GWL) and the lake water level (LWL) should there be one.There are 644 whole squares in the case of the Nechako watershed.
Each whole square is itself divided into "partial squares" according to the subpixel drainage divide.There are a total of 1082 partial squares in the case of the Nechako watershed.Available water from whole squares (whether from surface runoff or drained from the soil and/or groundwater reservoirs) is divided into these partial squares according to the fraction of whole square area drained by partial squares to form volumes (VOLs).These volumes represent the total amount of water available for transfer from one partial square to the next.The actual amount of water transferred over a given period is called streamflows and is defined as follows: where SF j is the streamflow at partial square j , M j is the number of partial squares directly upstream, ext k is a transfer coefficient and t is the time step.VOL is therefore a state variable, but streamflow, like SCA, is not considered to be a state variable since it has no impact on future simulations if its value is tampered with.

Ensemble Kalman filtering
The EnKF is a data assimilation method developed by Evensen (1994).It is an approach often used in hydrology, mainly due to its ability to consider non-linearities in the model and its relative simplicity to implement.The EnKF is a sequential method, meaning it relies only on current observations to update state variables as opposed to non-sequential approaches such as smoothers (Evensen and van Leeuwen, 2000) and recursive methods (McMillan et al., 2013).The EnKF propagates an ensemble of model runs based on a Monte Carlo implementation to represent model errors.The model covariance matrix (P b t ) at a time t is computed from the state vector (x b t ) holding the N ensemble members and their simulated variables; and the ensemble mean of the state vector (x b t ) therefore implicitly taking the model dynamics into consideration: when an observation is available, it is perturbed to form an ensemble of observations that are used to update each ensemble member.The updating step applies the Kalman gain (K t ), which is computed from observation (R t ) and model covariance matrices as well as an observation operator (H t ), which relates the model states to the observation: The Kalman gain acts as a weighted average between the observation and state vector to yield a post-filter analysis (x a t ) computed as such: The EnKF has practical and theoretical limitations.First, the EnKF relies on an ensemble representation of model and observation errors that are valid in the limit where ensemble sizes approach infinity.This is not feasible in practice, so a finite sample is used instead which aims to be sufficiently large such that sampling errors are negligible while ensuring that computational power and memory limitations are met.The method also makes use of model and observation covariance matrices to compute the gain during the updating process.These covariance matrices assume a linear relationship between variables.The EnKF also assumes normally distributed, bias-free and time-independent errors for both the model and the observations.Since these assumptions are not always met, which means that optimality is not guaranteed, a synthetic experiment is recommended to test the applicability of the method to the specific case being studied.
3 Experimental design

Synthetic experiment
Synthetic experiments, such as the ones done by Xie and Zhang (2010) or Weerts and El Serafy (2006), are test beds used to test the robustness of a data assimilation method or to tune various hyper-parameters.This is because the true state is known since it is initially created from true input that is also known.
For the current study, interpolated meteorological input and a specific set of parameters were used to run CEQUEAU, the output of which was considered to be the true state (step 1) (see Fig. 3).Synthetic observations, which include daily streamflow, SWE and SCA (step 2), and meteorological input, which include daily mean air temperature and precipitation (step 3), were then obtained by applying a perturbation to the true state and true meteorological input respectively.This means that the observation sets described thus far are not directly used, but synthetically generated using known parameters and perturbation.The synthetic observation ensemble (step 4) and meteorological ensemble (step 5) were created by further perturbing the synthetic observations and meteorological input.Ideally, these meteorological and observation ensemble perturbations should reflect the true errors of the synthetic meteorological input and observations for an optimal analysis.An ensemble of hydrologic states (step 6) was then obtained by running CEQUEAU using the synthetic meteorological ensemble.The EnKF was then applied using both the model and observation ensembles to produce an analysis (step 7), which was used as an initial state to produce ESP (step 8) using the true meteorological input.Additional details pertaining to the procedure used in generating perturbations and ensembles are included in the following sections.

Synthetic observation perturbation
Three types of observations were considered; namely streamflow, SWE and SCA.These synthetic observations were generated using a daily time step since their real-world counterparts are usually available on a daily basis.
In order to abide by EnKF assumptions, observations errors should ideally have a normal distribution.However, this is not practical due to the physical limits of the observations.For example, SWE observations cannot be negative and adding a normally distributed perturbation to SWE could result in some values being negative.Raising the negative values to zero or above would introduce a bias.Therefore, other distributions that share similarities with a normal distribution, while ensuring that physical limits are respected, were used to generate synthetic observations.Synthetic watershed-wide SCA were created using perturbations that follow a beta distribution since SCA is bounded between 0 and 1. SCA observations are expressed as y t,j ∼ B −1 Q t,j | α t,j , β t,j , where Q t,j is the cumulative probability of a temporally correlated normal random field with zero mean and unit variance at time t for observation j , and α t,j and β t,j are positively valued shape parameters.The shape parameters may be expressed in terms of the mean µ t,j and variance σ 2 t,j , but it must follow that σ 2 t,j < µ t,j 1 − µ t,j .The variance was arbitrarily set to σ 2 t,j = µ t,j 1 − µ t,j /50, such that the resulting shape parameters are α t,j = 49µ t,j and β t,j = 49 1 − µ t,j .Examples of beta distributions for different means using the same definition of variance as described above are shown in Fig. 4a.The distribution has a null variance and greater deviation from a normal distribution when the snow covers either 0 or 100 % of the watershed, as well as a variance at its greatest and resembling most a normal distribution at 50 % SCA.This approach avoids introducing a systematic bias when assimilating extreme values of SCA.When values are at 0 % (or 100 %), perturbations can only introduce higher (or lower) values in order to remain within the physical limits of the observations.This approach also gives the observations a greater uncertainty during the transition periods when SCA is lower than 100 % and higher than 0 %, which loosely follows the greater uncertainty attributed to MODIS observations over the same periods (Hall and Riggs, 2007).
Synthetic streamflow and SWE observations were created using perturbations that have a lognormal distribution since both observations are bounded to the left at zero and are theoretically unbounded to the right.Observation values y t,j can be expressed as y t,j ∼ ln N −1 Q t,j | µ t,j , σ 2 t,j .The true state was used for µ t,j for both streamflow and SWE, while Gamma, lognormal and normal distributions with µ = 1 the relative variance σ 2 t,j was set to 20 and 10 % respectively.The error distributions using these parameters are shown in Fig. 4b.The exact value of these variances is arbitrary for feasibility purposes.However, since the conclusions of this study will likely be used to help set up real-world applications, the variances chosen should ideally be relatively similar to the error of their corresponding real observations.Since these real observation errors are not known, rough estimates are used.
The use of these distributions is a compromise between the normal distribution of observations required by the EnKF and the physical limits of the observations without introducing a bias.

Meteorological input perturbation
Both the true daily precipitation and temperature values were perturbed using a gamma distribution, which has the benefit of generating positive values exclusively.Perturbations were implemented such that the meteorological input z t,i (precipitation or temperature) at time t over the whole square i is the result of the inverse gamma function given the cumulative probability P t,i of a spatially and temporally correlated normal random field with zero mean and unit variance.This can be expressed mathematically as z t,i ∼ −1 P t,i | κ t,i θ t,i , where κ t,i and θ t,i are shape and scale factor respectively.The shape and scale factors can be expressed in terms of mean µ t,i and variance σ 2 t,i , such that κ t,i = µ 2 t,i / σ 2 t,i and θ t,i = σ 2 t,i / µ t,i .In this study, synthetic precipitations are generated using the value of the true precipitation for µ t,i and a relative variance of 50 %, such that σ 2 t,i = 0.5 • µ t,i .Figure 4 shows the resulting error distribution using these parameters.Similarly, perturbed temperatures use the true temperatures for µ t,i and a standard devia-tion of 1 • C. Within the synthetic study where the feasibility of the approach is tested, the exact value of these perturbations is arbitrary, so long as it is coherent between scenarios.The values used were such that the Nash-Sutcliffe efficiency (Nash and Sutcliffe, 1970) of the simulated streamflow resulting from CEQUEAU using the perturbed meteorological input compared with synthetic streamflow was roughly similar to the performance of the simulated streamflow using real-world meteorological input compared with real streamflow observations.

Ensemble streamflow predictions generation
ESPs were generated using the ensemble of state variables resulting from the EnKF as initial states, true meteorological input and true model parameters.Using the true meteorological input implies that over a sufficiently large forecast horizon, every DA scenario considered in this study is likely to converge to the true state, but at different rates.By comparing the relative gains in performance over the ensemble with no data assimilation (open loop), one can then observe the length of time upon which DA impacts ESPs without having erroneous meteorological input affecting the results.This still generates an ensemble of streamflows since each ensemble member has its own initial states (VOL, SWE, SRI, STI, SML, GWL, LWL).
ESPs were generated everyday over the entire study period (10 years) using a forecast horizon spanning 50 days.

Hyper-parameter tuning
The use of the EnKF requires the tuning of hyper-parameters, such as model and observation errors, and ensemble size.Improper specification of these hyper-parameters could lead to filter divergence (Houtekamer and Mitchell, 1998).

Ensemble size
The ensemble size should ideally approach infinity to reduce the impact of sampling when covariance matrices are computed, but this is not feasible given the limits of computing power and memory.In practice, the ensemble size is chosen such that computing time is more reasonable while ensuring that the sampling error remains small.
Tests were carried out using ensemble sizes of 8, 16, 32, 64 and 128 members.An ensemble size of 64 members was used for this study.This number was chosen as a function of the stability between successive runs and computing resources available.It was found to be a reasonable trade-off between having sufficiently consistent results between simulations, such that the sampling error would be dwarfed in comparison with the impact of the actual data assimilation, without exceeding the computing resources available.

Meteorological ensemble generation
Perturbation factors similar to the ones used to generate synthetic meteorological inputs were used to generate an ensemble spread.This means that at every time step, meteorological ensemble (z t,i ) were generated using an inverse gamma function given the cumulative probability P t,i of a spatially and temporally correlated normal random field with zero mean and unit variance, mathematically expressed as z t,i ∼ −1 P t,i | κ t,i , θ t,i , where the shape factors are defined by κ t,i = µ 2 t,i / σ 2 t,i and θ t,i = σ 2 t,i / µ t,i .The prime symbol is used to distinguish between the ensemble variables/parameters and the synthetic variables/parameters. Precipitation was generated using the value of the synthetic precipitation (z t,i ) for µ t,i and a relative variance of 50 %, such that σ 2 t,i = 0.5 • z t,i , while temperature ensembles were generated using synthetic temperatures for µ t,i and a standard deviation of 1 • C. Using similar perturbation factors between synthetic and ensemble versions of the meteorological input reduced the probability of filter divergence cause by a misrepresentation of the model error.Errors from CEQUEAUspecific parameters were not taken into consideration, such that the parameter set used for the generation of the true state were the same for the ensemble generation.

Observation ensemble generation
As with model error representation, the perturbation factors used to generate an ensemble of observations were similar to the ones used to generate synthetic observations.Streamflow and SWE observation ensembles were created using perturbations that have a lognormal distribution centered around the synthetic observations z t,j with a relative variance of 20 and 10 % respectively.Watershed-wide SCA ensembles were created using a beta distribution centered around the synthetic observation z t,j with a variance of σ 2 t,j = z t,j 1 − z t,j / 50.Using similar perturbation factors avoided problems caused by a misrepresentation of the observation errors.

Covariance localization
The main disadvantage in using a finite sample to compute covariance matrices is that the resulting covariance matrices are not exact.This may result in theoretically zero covariance elements between two theoretically uncorrelated variables to become small, but non-zero, which may deteriorate the performance of the EnKF.
One way to overcome this dilemma is to use covariance localization, where covariances are forced to zero between some variables.One option is the Schur product where a covariance matrix is multiplied element-wise by a distancedependant correlation function (Houtekamer and Mitchell, 2001).However, there are other geophysical characteristics, such as land cover and elevation, which could be considered Hydrol.Earth Syst.Sci., 20, 4375-4389, 2016 www.hydrol-earth-syst-sci.net/20/4375/2016/ in the covariance localization.This would further increase the number of parameters to set and the degree of subjectivity in setting those parameters when the degrees of dependence are unknown.
Another approach was used in this study, which is based on the improvements observed in the state vector.First, the open loop is executed, as well as a data assimilation scenario with one observation and the corresponding spatialized state variable included in the state vector (e.g. 1 snow pillow assimilated and all modelled SWE included in the state vector).Then, the two runs are compared with the true state on a spatial basis.In the case of CEQUEAU, these can be whole or partial squares depending on the variable analyzed.The covariance matrix is localized such that the areas that do not show an improvement for the data assimilation scenario over the open loop are set to zero.This process is repeated for each observation.
While this process remains susceptible to the sampling error from the finite ensemble size, it is a simple approach that exploits the availability of the true state in a synthetic experiment and limits the state vector size according to observed improvements.
In this study, only SWE observations have a corresponding state variable, so covariance localization has only been applied to the SWE variable.

State vector configuration
Though the state vector often comprises only of the variables corresponding to the observations or those judged to be relevant enough by the user, there are potentially many state variables that could benefit from the assimilation of available data if there exists a linear (or approximately linear) relationship between the modelled variables and the observations.
To determine which variable could benefit from being included in the state vector, one could execute multiple scenarios where each possible combination is compared with the true state.However, this could get very laborious even for a relatively small number of state variables.The current approach suggests reducing this number by first adding state variables one at a time.The variables that show a global improvement can then be added to the state vector.Assuming that not all variables are added to the state vector, this reduces the number of combinations to try.

Metrics
Various metrics were used to quantify results.The mean square skill score (MSSS), based on the mean square error (MSE), was used to assess the differences between various data assimilation scenarios and the open loop during the state vector configuration and covariance localization processes.The MSE for a variable of interest x is defined as where N is the number of time steps, x t is the ensemble mean analysis of the state variable of interest at time t and x T t is the corresponding true state.It is often more convenient to express this score as a unitless skill score: The ensemble forecast performance was assessed using the continuous rank probability score (CRPS; Hersbach, 2000) and its associated skill score (CRPSS).For this synthetic study, the CRPS is adapted as follows: where F x f t and F x T t are the cumulative distribution function of the ensemble forecast at a horizon f and the true state, respectively.The CRPS has the same units as the variable of interest and is bounded by [0, +∞].A lower CRPS is a better score.As with the MSE and MSSS, it is often convenient to express the CRPS in its skill score form: 4 Results and discussion

State vector configuration and covariance localization
Before investigating the effect of data assimilation on streamflow forecasts, a state vector configuration analysis was conducted.This was done in order to find out which variables, among the seven listed previously (VOL, SWE, SRI, STI, SML, GWL, LWL), should be included in the state vector for each type of data assimilated in order to reduce the number of comparisons to make.

Streamflow data assimilation
First presented are the results from the case where only streamflow at the outlet was assimilated.Streamflow at the outlet is computed by the model, but it is output only.Therefore, in order for the assimilation of streamflow to have any impact on the modelled states, additional variables needed to be added to the state vector.
Figure 5 shows a box plot of the MSSS computed for each variable on the whole watershed when they are individually included in the state vector using the open-loop scenario as a reference.Values above zero mean there is an improvement for a particular partial (for volumes) or whole (for other state variables) square compared with the open loop.The boxes range between the 25th and 75th percentiles, with a red bar to show the median, and the whiskers range between the maximum and minimum values.Outliers are not shown for visibility purposes.Results for the case where water VOLs are included along with the streamflow at the outlet show an improved score for each partial square on the watershed.This is not entirely surprising given the close relationship between streamflow and volume.This suggests a necessity to include VOL in the state vector when assimilating streamflow at the outlet for streamflow predictions.
Results also show a deterioration of SWE for nearly 75 % of whole squares on the watershed.Although there is some improvement for some whole squares, this suggests that including SWE in the state vector when assimilating streamflow at the outlet is unlikely to be beneficial for streamflow predictions.Although SWE does have an important impact on streamflow, there is a time lag between the snowmelt occurrence and the increase of streamflow at the outlet.Since the EnKF assumes linear relationships between variables, the non-linearity between SWE and streamflow can result in a non-optimal analysis.In this case, the results are actually worse than open loop for most whole squares.Clark et al. (2008) discussed the issue of non-linearities between streamflow and other variables.To overcome this issue, one could use either a recursive approach, which allows for adjustments of previously simulated variables, or a smoother approach to DA, which also uses "future" observations to update current state variables.However, this may not be necessary given the positive impact of streamflow DA on VOL, as well as in a multivariate DA scenario where other variables, such as SWE in this case, are also assimilated.
As for the SRI and STI, the median sits around 0, which means there is no improvement for 50 % of the whole squares.This suggests little change can be obtained in the analysis by including those variables in the state vector.Similar to the case with SWE, there is likely a time lag issue between streamflow and these variables.However, there is also a weaker link between these variables such that a change in SRI or STI is not as strongly linked to an eventual change in streamflow as much as it is for a change in SWE.
Finally, results for the three conceptual reservoirs SML, GWL and LWL show an improvement for over half of the whole squares, with a greater number of whole squares improved for the GWL and slightly above zero median for SML.This suggests that including these three variables in the state vector can potentially yield improvements for streamflow predictions.Though the relationship between the water level in these conceptual reservoirs and streamflow at the outlet is not exactly linear, mainly due to reservoirs having multiple orifices (see Fig. 2) and the time lag before water reaches the outlet, it may be sufficiently near linear such that streamflow DA yields an overall improvement for most whole squares.For example, the median correlation coefficient of a simple linear regression between each reservoir for each whole square and the streamflow at the outlet is 0.12, 0.49 and 0.15 for SML, GWL and LWL respectively.Samuel et al. (2014) and Trudel et al. (2014) found that updating soil moisture with streamflow observations actually deteriorated soil moisture simulation compared with real soil moisture observations.However, there are notable differences between these studies and the present one.Aside from the different features of the study area and model structure, the use of synthetic data instead of real data likely strengthens the link between variables and observations.Since synthetic observations are constructed using the same model and parameters as the model in which the observations are assimilated, there is no difference in scale between observations and modelled variables, which is often an important source of error for studies using real data.
Nonetheless, the inclusion of VOL, SML, GWL and LWL in the state vector were considered during the assimilation of streamflow at the outlet.The impact of each scenario for streamflow predictions are compared in Sect.4.2.1.

SWE data assimilation
The same analysis was performed for SWE data assimilation from synthetic snow pillows.However, unlike streamflow, SWE is a state variable such that any changes made upon it will have repercussions on future simulations.SWE at the location of the snow pillows should therefore be included in the state vector and also potentially whole squares in the vicinity that are correlated with these locations.A spatial analysis was performed first to determine the spatial ex- tent upon which each snow pillow may affect modelled SWE in other whole squares.
Figure 6 shows the MSSS of SWE on the spatial level for the snow pillows located at Mount Wells, Mount Pondosy and Tahtsa Lake, using the open-loop scenario as a reference.For each figure, the whole square that shows the most improvement is the area where the corresponding snow pillow is located.Whole squares that show improvements are mainly located around snow pillows, but the range differs for each snow pillows.Various areas in remote locations also show improvements for each snow pillow.As mentioned earlier, relationship with geophysical factors, such as distance from snow pillow, elevation and land cover, could be used to ex- plain this variation, but a simpler approach was used such that the covariance localization was limited to whole squares showing improvements only.The covariance elements representing all the other whole squares were set to zero.
As for the state vector configuration, Fig. 7 shows the MSSS computed for each variable within the extent of whole squares, which were positively impacted by SWE DA during the covariance localization process.The open-loop scenario was used again as a reference.The results show no significant improvement for any other variable except for SWE itself, which yields only positive MSSS values by design.The lack of overall improvement for water-related variables (VOL, SML, GWL, LWL) is coherent with the time delay with changes in SWE.As for the other snow-related variables (STI, SRI), although there may be a relationship with SWE, it is non-linear (US Army Corps of Engineers, 1956), which is further weakened by the distance separating SWE at a snow pillow from STI or SRI at another location.
Only the inclusion of SWE surrounding a given snow pillow in the state vector is considered during the assimilation of SWE for streamflow predictions in Sect.4.2.2.

SCA data assimilation
Like streamflow, SCA is not a state variable.It is computed in parallel with CEQUEAU without having any direct effect on future simulations.In order to have any impact during the assimilation process, there must exist a linear or sufficiently near-linear correlation between SCA and state variables.The update step should bring improvements to the state variables if the computed correlation also reflects the true correlation.
Figure 8 shows a box plot of the MSSS computed for each variable when they are individually included in the state vector.The open-loop scenario is used as a reference.Results show that global snow cover data assimilation yields little or no improvement for all state variables compared to the open-loop scenario.For most cases, a strong deterioration is observed, suggesting that the way SCA data are used in this study is not well adapted for the current assimilation pur- poses using the EnKF.Marcil et al. (2016) have shown that there exists a relationship between the SCA and the percentage of cumulated streamflow at the outlet, but it is neither linear nor is cumulated streamflow a state variable.The EnKF requirement that relationships between variables be linear and synchronized severely limits the value of global SCA data for the current application.This result is coherent with the findings presented by Clark et al. (2006).Using a model which incorporates snow cover area as a state variable, such as the snowmelt runoff model (SRM; Martinec, 1974) or the Soil and Water Assessment Tool (SWAT; Arnold et al., 1998), could overcome the issue of non-linearities between variables, while using recursive or smoother approaches to data assimilation could help with the time lag issue between observations and state variables.
Given the absence of overall improvement for all the state variables, the impact of SCA DA on streamflow predictions was not considered in this study.

Streamflow forecasts
Aside from granting insight into the sensitivity of the system to the state vector configuration, the analysis in the previous section presented a list of state vector configurations likely to favour streamflow predictions improvements based on the improvement of various state variables.This section presents ensemble streamflow prediction results for each configuration selected for each type of data assimilated.

Streamflow data assimilation
Focusing on the case where only streamflow at the outlet are assimilated, Fig. 9 presents the CRPSS of predicted streamflow at the outlet over a forecast horizon of 50 days using the open loop as a reference.Only the state vector configurations that showed some improvements in the state vector configuration analysis section are shown.
First, high values of CRPSS for short-term forecasts can be observed for the case where only volumes are included in the state vector (blue curve).The CRPSS subsides asymptotically to zero over time, which shows assimilating streamflow to update volumes improves streamflow predictions compared to the open loop only for a few days, after which the impact of streamflow assimilation becomes insignificant.The duration of the impact depends on the residence time of the water stored in the model's partial squares (VOL).A relatively short-lived impact would mean a relatively short residence time.The high initial impact is not surprising given the nearly linear relationship between streamflows and volumes.Assimilation of streamflow observations generate a globally positive update on volumes, as seen in Fig. 5, which in turn strongly affects simulated streamflows.
Second, adding each of the three water reservoirs individually to the volumes yields different results.Even though lake water levels showed improvements over the majority of whole squares, the impact on streamflow predictions (red curve) is marginal compared to the case where only volumes are included in the state vector.This is because the weights attributed to lakes in CEQUEAU are very low for most whole squares.Only about 0.5 % of the entire watershed is modelled using the conceptual lake reservoir and its parameters, unlike the soil moisture and groundwater reservoirs, which are present in every whole square.Adding SML (green curve) or GWL (magenta curve) instead of LWL not only increases the initial CRPSS, but also slows the decrease over time.This is consistent with the improvements observed for the updated water levels for over half of the whole squares compared with the open-loop case, which translate as added improvements over the case where only volumes are included in the state vector.The slower decrease over time is also coherent with the increase in time it takes for water in the reservoirs to reach the outlet compared with water already in the routing system.The groundwater reservoir is shown to have an initially similar, but longer-lasting positive impact than the soil moisture reservoir.The soil moisture reservoir controls mainly the fast-flowing surface runoff, the amount of evapotranspiration leaving the system and the amount of water infiltrating into the groundwater reservoir.The groundwater reservoir has a numerically unlimited capacity, with no way out for the water except through evapotranspiration and the outlets that feed the routing system, making its impact on streamflows last longer than the relatively ephemeral soil moisture reservoir.
Finally, the scenario where all four variables are added to the state vector is analyzed (orange curve).The difference noted with the other curves is mainly caused by the simultaneous inclusion of SML and GWL, since all the other curves already have VOL included, while LWL was already shown to have very little impact on streamflow at the outlet.Comparing the fully combined case with the VOL+GWL case, although the initial improvement is similar between the two, with the former slightly above, the latter has a slower decrease over time.This suggests that the addition of SML interferes with the GWL update.As seen for the state vector configuration analysis (Fig. 5), the assimilation of streamflow at the outlet had a positive impact on a greater number of whole squares for the GWL than the SML.Here, the increased number of deteriorated SML, which infiltrates into the groundwater reservoirs, hinders the GWL updates such that the results show some deterioration compared with the VOL+GWL case, even though it is still an improvement over the VOL only updates.
These results have some similarities and differences with other studies.For example, Abaza et al. (2015) assimilated streamflow at the outlet using the EnKF to update two state variables (soil moisture in the intermediate and deep layers of the hydrological model used in their study) using a time step of 3 h.The resulting gain in CRPS was high for the first time step and decreased quickly as a function of the forecast horizon such that mainly the first 24 h benefited from the data assimilation.Chen et al. (2013) found similar results with multiple performance criteria when assimilating streamflow using a variant of the EnKF (the ensemble square root filter) to update various state variables represented by conceptual reservoirs.The observed improvement duration was even shorter, lasting less than 12 h during flash flood events.The difference in impact duration is likely related to the different water retention time in each watershed.These studies were conducted on much smaller watersheds (all less than 800 km 2 ) than the Nechako watershed (around 14 000 km 2 ), further highlighting that the performance of assimilation techniques is related to watershed characteristics.

SWE data assimilation
Following the same method as with streamflow, this section focuses on the case where only SWE from snow pillows were assimilated.Since the state vector configuration analysis showed only improvements for the SWE variable, it was the only variable added to the state vector for streamflow forecast.However, since there are three observations available, Fig. 10 presents the CRPSS of predicted streamflow at the outlet when assimilating SWE from each snow pillow individually and collectively.
An interesting result is that the impact of each snow pillow on streamflow predictions varies greatly.The impact of the snow pillows located at Mount Wells (green curve) and Mount Pondosy (blue curve) are dwarfed in comparison with the impact of the snow pillow located at Tahtsa Lake (magenta curve).This is coherent with results from Marcil et al. (2016) over the same watershed.The lower impact of the Mount Pondosy snow pillow is explained by the relatively small region of influence observed in Fig. 6b.As for Mount Well, even though it has the largest area of influence (Fig. 6a), it is also the snow pillow affecting the regions with the lowest altitudes and also the least amount of maximum SWE.Although the region affected by the Mount Wells snow pillow contains a mean annual maximum SWE of 410 mm, it is 40 % less than for the region affected by the Tahtsa Lake The forecast horizon varies from 1 to 50 days.Lack of parentheses indicates that the variable is affected by both types of observations.snow pillow, which contains a mean annual maximum SWE of 682 mm.
Nonetheless, assimilating all three snow pillows yields better results for mid-term streamflow forecasts (Fig. 10, red curve).Even though the Tahtsa Lake snow pillow carries the most importance, the other snow pillows have a positive effect on regions that are not reached by the Tahtsa Lake snow pillow area of effect.The assimilation of all three snow pillows does yield short-term forecasts improvements, but a better score is reached over time.This is because the impact of SWE over streamflow at the outlet occurs during snowmelt, which can occur at a much later date than when SWE observations are assimilated.Although the curve gives the impression of a monotonous increase over time, this is only due to the limit imposed on the forecast horizon.At a further horizon, the curve should eventually peak and decrease asymptotically to zero since there is no accumulation of snow from year to year.Over time, the simulation should eventually become indistinguishable from the open-loop scenario.Franz et al. (2014) also assessed the impact of SWE data assimilation on ensemble predictions, but using real observations.Their results showed little improvement of forecast performance through SWE data assimilation, but they highlighted the role of a possible bias in the observations, as well as the difference in scale between the point-scale observations and the basin average SWE simulated by the model.In this synthetic experiment, no bias was specified on observations and there is no difference in scale between modelled and observed SWE, which could explain the differences observed between the two studies.Biased observations and meteorological input were purposely omitted in this study, but may be added in future works to test the robustness of the approach.

Combined streamflow and SWE data assimilation
The focus now shifts to the case where streamflow at the outlet are simultaneously assimilated along with SWE from the three snow pillows.The state vector configuration, which provided the best results from the streamflow data assimilation case, is used (VOL + GWL) along with the best configuration from SWE data assimilation (SWE only).Although these configurations worked best with their respective data assimilation case, they could behave differently when both streamflow and SWE are assimilated together.
Table 1 presents four configurations for the combined assimilation of streamflow and SWE observations.These configurations differ in the overlap of their effect during the update phase such that some configurations allow both observations to simultaneously update the same variable, while others do not.
The performance of these configurations on the CRPSS for predicted streamflow at the outlet is presented in Fig. 11.While all four configurations perform in a very similar way for short-term streamflow predictions, the group forms two pairs that differ in that the blue-magenta group allows SWE observations to update modelled streamflow, while the greenred pair does not.Although allowing SWE data assimilation to update VOL and GWL changes very little, a drop in performance occurs if streamflow assimilation updates modelled SWE.This is coherent with the state vector configuration analysis performed in the previous section (Figs. 5 and 7), where SWE data assimilation is shown to have a weak impact on VOL and a median MSSS around 0 for GWL, while streamflow assimilation deteriorated around 75 % of SWE whole squares when they were included in the state vector.
Overall, the simultaneous assimilation of streamflow and observed SWE yields important improvements over the entire forecast horizon analyzed, with the streamflow data assimilation improving mainly short-term streamflow forecasts and SWE data assimilation improving mainly mid-term streamflow forecasts.CRPSS values for combined assimilation of both streamflow and SWE observations were superior to CRPSS values for individual assimilation of streamflow or SWE over all forecast horizons, with the exception of forecast horizons higher than 45 days, where CRPSS values for SWE DA are slightly higher.This reveals that the updated VOL and GWL by streamflow data assimilation may be very beneficial for short-term forecasts; furthermore, they do not further improve the mid-term forecasts when combined with SWE data assimilation in comparison with the scenario where only SWE data are assimilated.
Moreover, the assimilation of each data type (streamflow and SWE) differs not only by their impact over the forecast horizon, but also over the time of the year.Figures 12 and 13 show the monthly CRPS for the open loop (black curve), the streamflow data assimilation including VOL and GWL (blue curve), the SWE data assimilation of all snow pillows (red curve) and the simultaneous, but separated, streamflow and SWE data assimilation (VG-S; green curve) for short-term (average of horizon from 1 to 5 days) and mid-term (average of horizon from 25 to 50 days) streamflow forecasts.The CRPS is shown for the open loop to show the performance change over the time of the year and the period when improvement is most needed.Recall that the CRPS ranges from zero to infinity, with zero being a perfect forecast.The period from May to July, which corresponds to the melt period, is therefore the period when the CRPS is the highest for shortand mid-term forecasts are the most problematic.The scores for mid-term forecasts are lower than for short-term forecasts because the true weather is used as input for forecasts such that the open loop slowly converges to the true states over time.
Assimilating streamflow results in an improved score over the entire year for short-term forecasts, although little gain is obtained for mid-term forecasts.This steady improvement is to be expected since streamflow here is always non-zero and observations are available all-year round.On the other hand, the impact of the assimilation SWE from snow pillows is limited mainly to the melt period for both short-term and mid-term forecasts.However, this period corresponds to the problematic period when most gain can be obtained.The assimilation of SWE provides a better score than the assimilation of streamflow for the same period and improves going from short-term to mid-term forecasts.SWE assimilation complements streamflow assimilation as observed from the performance of the simultaneous assimilation, which yields both the steady improvements over the year and the important gain during the snowmelt period.

Conclusion
This study investigated the impact that multivariate data assimilation can have on streamflow forecasts using the CE-QUEAU hydrologic model applied over the Nechako watershed in a synthetic experiment.The study also showed the importance of the state vector configuration on streamflow forecasts when using the EnKF.
Streamflow data assimilation was found to improve shortterm streamflow forecast considerably.However, the impact dissipated relatively rapidly as a function of the forecast horizon, which was slowed by adding groundwater conceptual reservoir levels to the state vector.Improvements were observed for all months of the year: low-flow and high-flow periods alike.
On the other hand, the assimilation of snow water equivalent data from synthetic snow pillow data yielded streamflow forecast improvements mainly during the snowmelt period.Although the period lasts approximately 3 months, the impact was found to be greater than streamflow data assimilation over the same period.It was also noted that assimilat-ing each snow pillow data individually yielded different results, with various radii of influence, such that the improvement from assimilating all three snow pillows simultaneously covered most of the watershed and yielded streamflow forecasts which outperformed forecasts from any single snow pillow data assimilation.Over the forecast horizon, the peak of improvement was greater than or equal to the 50-day limit over which forecasts were simulated, which contrasts with the short-lived impact of streamflow data assimilation.
Given their complementarity, streamflow and snow water equivalent data were assimilated simultaneously.The resulting streamflow forecast inherited the strengths from both types of data, having a strong, positive impact for both shortterm and mid-term forecasts.Improvements were obtained for all periods of the year, but mainly during the snowmelt period, which is normally the most problematic.
The assimilation of basin-wide snow cover area failed to improve the simulation of any state variable.The most probable factor was determined to be the absence of snow cover area as a state variable or a proxy with a sufficiently linear relationship with SCA.Suggestions to improve the method to accommodate for snow cover area are to use a model that incorporates snow cover area as a state variable and/or to use a data assimilation approach, which takes into account a time lag between observations and state variables.
The results obtained are conditional to some assumptions and limitations.First, all results depend on the general method and parameters used in creating the synthetic framework.Since this is a synthetic experiment, it is assumed that a real experiment would behave similarly to a simulation using CEQUEAU with a specific set of parameters and inputs.Second, the potential impact of data assimilation on streamflow forecasts observed depended on using the true weather inputs.Using real weather inputs may decrease this impact.Third, it is assumed that the error representations for the model inputs and the observations are known.In this study, they have been generated using specific distributions and variances to compromise between the need for normal distributions and the need to remain within the physical limits of the variables without introducing a bias.Finally, the impact of errors from the model parameters is assumed to be negligible, such that the set of parameters was not altered from the true simulation's set of parameters.
The assumptions and limitations listed reveal several challenges posed by the assimilation of multiple types of observations for streamflow forecasting purposes.Future work includes investigating into the performance of multivariate data assimilation in the presence of biases and unknown errors, as well as the economic impact of streamflow forecasts generated with multivariate data assimilation on real management practices.

Figure 2 .
Figure 2. Diagram of the processes included in CEQUEAU's production function.

Figure 3 .
Figure 3. Flowchart for the production of ensemble streamflow predictions obtained from various data assimilation scenarios.

Figure 4 .
Figure 4. Examples of (a) beta and (b) gamma and lognormal distribution compared with their analogous normal distributions using the same mean (µ) and variance (σ 2 ).
where MSE ref (x) is a mean square error of reference; the open loop in this case.The MSSS is bounded by [−∞, 1] and indicates an improvement as the skill score increases.Values above zero indicate an improvement over the reference (open loop) and a value of one indicates a perfect score: a perfect correspondence between the mean of the analysis and the true state.
where CRPS ref(xf )  is the continuous rank probability score of the open loop used as a reference in this case.Like the MSSS, the CRPSS is bounded by [−∞, 1], with higher values indicating a better score.Values above zero indicate an improvement over the reference and a value of one indicates a perfect score.

Figure 5 .
Figure 5. Box plot of the mean square skill score for each variable when assimilating streamflow at the outlet.The open loop is used as a reference.Outliers are not shown for visibility purposes.

Figure 6 .Figure 7 .
Figure 6.Distribution of the mean square skill score of SWE over the watershed when assimilating SWE located at (a) Mount Wells, (b) Mount Pondosy and (c) Tahtsa Lake.The open loop is used as a reference.Values below −1 are cut off from the legend.

Figure 8 .
Figure 8. Box plot of the mean square skill score for each variable when assimilating basin-wide snow cover area.The open loop is used as a reference.Outliers are not shown for visibility purposes.

Figure 9 .
Figure 9. Continuous rank probability skill score of the streamflow ensemble when assimilating streamflow at the outlet.The open loop is used as a reference.The forecast horizon varies from 1 to 50 days.

Figure 10 .Figure 11 .
Figure 10.Continuous rank probability skill score of the streamflow ensemble when assimilating SWE from all three snow pillow locations.The open loop is used as a reference.The forecast horizon varies from 1 to 50 days.

Figure 12 .Figure 13 .
Figure12.Continuous rank probability score for short-term forecasts (average of forecast horizons 1 through 5 days) of various data assimilation scenarios as a function of the month of the year.

Table 1 .
Overview of multivariate DA scenarios.