Snow cover estimation using blended MODIS and AMSR‐E data for improved watershed‐scale spring streamflow simulation in Quebec, Canada

Estimation of the amount of water stored in snow is a principal source of error for spring streamflow simulations in snow‐dominant regions. Measuring this variable throughout large and often remote areas using snow surveys is an expensive task since they are practically point measurements. Remote sensing is an alternative method, which can cover much larger areas in little time, but further research is required to reduce uncertainties on snow water equivalent (SWE) estimations, especially during the melting period. However, optical‐near infrared (NIR) and passive microwave remote sensing can detect snow cover area (SCA) with greater certainty, which can be used as a proxy for SWE. The two datasets work in complementary ways considering their spatial resolutions and cloud cover limitations. This study developed an SCA product from blended passive microwave (Advanced Microwave Scanning Radiometer ‐ Earth Observing System: AMSR‐E) and optical‐NIR (Moderate Resolution Imaging Spectroradiometer: MODIS) remote sensing data to improve estimates of streamflow caused by snowmelt during the spring period. The blended product was assimilated in a snowmelt model (SPH‐AV) coupled with the MOHYSE hydrological model through a modified direct insertion method. SCA estimated from AMSR‐E data was first compared with in situ snow‐depth measurements and SCA estimated with MODIS. Results showed an agreement of over 95% between AMSR‐E‐derived and cloud‐free MODIS‐derived SCA products in the spring. Comparison with ground stations confirmed the underestimation of snow cover by AMSR‐E. Assimilation of the blended snow product in SPH‐AV coupled with MOHYSE yielded an overall improvement of the Nash–Sutcliffe coefficient comparable with simulations with no updates, which is comparable to results driven by biweekly snow surveys. Assimilation of remotely sensed passive microwave data was also found to have little positive impact on streamflow simulation due to the difficulty of differentiating melting snow from snow‐free surfaces. © 2013 The Authors. Hydrological Processes published by John Wiley & Sons Ltd.


INTRODUCTION
Snow cover monitoring is a major challenge in modeling surface water distributions. Since snow can act as a water reservoir controlling the water supply, snowmelt is an important variable for dam management because it is an excellent source of energy for hydropower stations (Lavallée et al., 2006), but it also poses a flood hazard in some areas (Yang et al., 2003). Spatial and temporal variations of snow can also have important consequences for water availability in snow-dominated regions (Barnett et al., 2005), such as the province of Quebec, Canada. A hydrological forecast system (HFS) based on the combination of a snow model and a hydrological model is thus used by the Centre d'Expertise Hydrique du Québec (CEHQ) in order to simulate and predict incoming streamflows for dam management purposes (Turcotte et al., 2004). Although improvements have been made through the direct assimilation of in situ snow depth and snow water equivalent (SWE) in the HFS, obtaining precise knowledge of snowpack properties, most importantly SWE, over the entire region poses a difficult challenge De Lannoy et al., 2012). Insufficient accuracy of the SWE variable can greatly affect streamflow estimates, especially during peak floods in spring caused mainly by snowmelt.
Satellite remote sensing can help to improve near real-time streamflow simulation by providing complete spatial coverage of the territory and through frequent data availability. Over the last few decades, visible remote sensing has been widely used to monitor snow cover area (SCA) for hydrological applications (see Seidel and Martinec, 2004). However, it has been shown that SCA products based on visible remote sensing used for hydrological purposes are limited by cloud cover . On the other hand, SCA products based on passive microwave remote sensing are practically unaffected by cloud cover, but they are sensitive to the presence of snowmelt and are provided at a much lower spatial resolution (Dietz et al., 2012). For operational applications, it is necessary to use an improved snow cover product capable of providing snow cover data on a daily basis independent of cloud coverage. Hence, previous authors, such as Ramsay (1998) and Foster et al. (2011) have suggested blending visible and passive microwave data to improve snow cover monitoring.
The main objective of this paper is to develop a method to monitor SCA during the melt season through a combination of different remote sensing data for operational watershedscale spring streamflow simulation. More specifically, we aim to (1) blend remote sensing SCA (RS-SCA) from visible and passive microwave data and (2) investigate the impact of assimilating this snow cover product and its individual components into the CEHQ's HFS for near real-time spring streamflow simulation. We first describe the hydrological and snowmelt models used in the HFS as well as the datasets used to map snow cover. The method is then presented and evaluated by means of cross-validation between snow cover products and comparison of HFS output with measured streamflows for three different watersheds in southern Quebec. Our results are discussed and concluding remarks given in the final section of the paper.

METHODS AND DATA
This section provides a general description of the models, the study area, the data, and the methods for processing the data that were used in this study.

Model description
This study uses the MOHYSE (Modèle Hydrologique Simplifié à l'Extrême) hydrological model (Fortin and Turcotte, 2006), which was adapted to use water from snowmelt and rainwater as input that is not retained by the snowpack as calculated by the snow model SPH-AV (Système de Prevision Hydrologique -Apports Verticaux).

SPH-AV.
Based on the hydrological model HYDROTEL (Fortin et al., 2001), SPH-AV is an empirical model based on a mixed degree-days-energy-balance algorithm for spatiotemporal simulations of a single-layered snowpack. The simulated snowpack simulates five state variables using precipitation and air temperature data as input, namely SWE, heat deficit, albedo, snow depth, and liquid water retention within the snowpack. Computations of these variables are made on 120 by 200 grids with pixels of dimensions 0.1˚by 0.1˚for three superimposed tiles corresponding to evergreen forests, deciduous forests, and open areas. The parameters for these layers (temperature threshold at which melting occurs, melt rate coefficient at the air-snow interface, melt rate coefficient at the soil-snow interface, maximum snow density, and snow compaction coefficient) have been calibrated to reflect these three dominant types of land surfaces in the southern part of Quebec. The modelled SWE and snow depth are also updated through direct assimilation of biweekly snow survey data. Precipitation and temperature data are measured at various permanent stations and interpolated on the model grid. Additional details on the model and the calibration can be found elsewhere in Turcotte et al. (2007).
MOHYSE. MOHYSE is a lumped, conceptual hydrological model that uses only the mean water from the snowpack calculated by SPH-AV and mean air temperature over the watershed as input data (Fortin and Turcotte, 2006). The model includes the following processes: evaporation, transpiration, infiltration, surface runoff, subsurface runoff, base-flow, and the transfer between the vadose and the aquifer reservoirs. The model contains several parameters for the following functions: potential evapotranspiration adjustment, transpiration adjustment, maximal infiltration rate, transfer rate from vadose to aquifer reservoirs, transfer rate out of the vadose reservoir, transfer rate out of the aquifer reservoir, a watershed shape parameter, and a watershed scale parameter. These eight parameters were calibrated for the watersheds in this study using the SCE-UA (Shuffled Complex Evolution-University of Arizona) method (Duan et al., 1992) to optimize the Nash-Sutcliffe coefficient (Nash and Sutcliffe, 1970) over the period from 1 January 2004 to 1 January 2006 using SPH-AV output with snow survey updates. Additional details on the model and the calibration can be found in Fortier-Fillion (2011) and Roy et al. (2010).

Study area
The study sites were three watersheds in southern Quebec: aux Écorces (1110 km 2 , 47.9˚N), which is a predominantly evergreen-forested sub-watershed of the Kénogami River; the northern part of the mainly deciduous-forested du Nord watershed (1170 km 2 , 45.8˚N) located northeast of Montreal; and au Saumon (738 km 2 , 45.4˚N), which is a watershed composed mainly of deciduous forest (Figure 1). In these areas, the snow season typically begins around early December and ends in April. Annual precipitation is approximately 1000 mm of which about a third is solid precipitation. A gauging station was located at each site and at least one snow survey station was also present; however, these stations were typically located in open forest areas. These watersheds were chosen because of past difficulty in accurately modelling the streamflow peaks at these sites.

Datasets
The datasets presented in this section were collected at the study sites over the 2004-2007 period. Unless specified otherwise, all the results presented in the Results section covered this period.
MODIS. Numerous SCA products from satellite data are made accessible through the National Snow and Ice Data Center (NSIDC). Two of these that are used in this study are daily products derived from Moderate Resolution Imaging Spectroradiometer (MODIS) data, a payload scientific instrument on board both the Terra and Aqua Earth Observing Systems (EOS) with 36 spectral bands ranging from visible to shortwave infrared wavelengths. These two products are the MODIS/Terra Snow Cover Daily L3 Global 500 m Grid Version 5 (MODT) and the MODIS/Aqua Snow Cover Daily L3 Global 500 m Grid Version 5 (MODA). These consist of 1200 km by 1200 km tiles of approximately 500 m resolution SCA data gridded in a sinusoidal map projection in a compressed Hierarchical Data Format -Earth Observing System (HDF-EOS). A resampling was required to match the resolution of SPH-AV, thereby aggregating not only SCA pixels but also cloudy pixels. Andreadis and Lettenmaier (2006) use 20% of fractional cloud cover as a threshold to decide whether or not to use the observation while stating that fractional cloud cover as high as 50% adds little information and does not significantly change the gridded SCA mean and variance. Roy et al. (2010) use a 30% threshold, which falls within the 20-50% range, to maintain the representation of cloud-free pixels while maximizing the number of resampled 'cloudfree' pixels. For our study, we chose a 30% threshold for similar reasons. The resampled cloud-free pixels represent a snow cover fraction (SCF), which was then converted back to binary data for data assimilation purposes (see section 'Accuracy assessment' below). A threshold of 10% was used for SCF to extend the snow cover duration (SCD), or the number of days where a pixel is considered snow covered, in spring to compensate for the general underestimation of modelled SCD . The imperfect overall accuracy of the MODIS SCA products, as shown by Hall and Riggs (2007), Parajka andBlöschl (2008), andSimic et al. (2004) where a 5-10% overall error is reported, was also considered in choosing the 10% threshold. The selection of a 10% threshold thus aims to minimize the threshold in order to increase SCD while keeping outside MODIS error range. Though not shown here, it was found that the precise value of this threshold played a minor role on results.
AMSR-E. This study also used satellite passive microwave data, which is also accessible through the NSIDC in HDF-EOS format. The Advanced Microwave Scanning Radiometer-Earth Observing System (AMSR-E) provides daily brightness temperatures (TB) at 6.9, 10.7, 18.7, 23.8, 36.5, and 89.0 GHz frequencies in both horizontal and vertical polarizations. These data are included in the AMSR-E/Aqua Daily L3 12.5 km Brightness Temperature, Sea Ice Concentration, and Snow Depth Polar Grids product on a 12.5 × 12.5 km grid on a daily basis. This product is itself created from AMSR-E/Aqua L2A Global Swath Spatially Resampled Brightness Temperatures data.
The microwave TB emitted by a snowpack at higher frequencies (19 and 37 GHz) is related to many physical characteristics, which are mainly SWE, snow density, grain size, and temperature (Grody, 2008). In order to extract snow cover from passive microwave satellite data, we used a threshold on 19 GHz and 37 GHz TB differences. This simple method has long been used to detect snow on the ground (Chang et al., 1987). The basis of the method is that the scattering effect of snow at 37 GHz is greater than at 19 GHz because of the shorter wavelength at 37 GHz (Mätzler, 1987). A significant difference in TB (ΔTB) can thus be interpreted as being caused by the presence of surface snow. Our method applied the threshold on ΔTB that minimized the total error or disagreement of SCA derived from AMSR-E data with MODT data for cloud-free days over the watersheds during spring periods in 2004 and 2005. The remaining springs of 2006 and 2007 were used for validation with data from MODT, MODA, and Environment Canada stations (section 'Environment Canada stations'). We chose to use MODT for optimization purposes because its potential for snow cover mapping during springmelt events has been previously demonstrated . Also, in order to minimize the effect of melting snow during the day, only the nighttime data (descending pass) was used.
As with MODIS SCA maps, the newly created SCA maps obtained from TB differences (AMSR-E SCA ) required a resampling to match the resolution of SPH-AV; here, the method was simply a nearest-neighbour resampling since the AMSR-E TB resolution is slightly lower than the SPH-AV grid.
Environment Canada stations. Environment Canada (EC) manages various weather stations throughout the province of Quebec, which includes the monitoring of snow thickness within the study zones. Historical data were extracted online from EC's National Climate Data and Information Archive website (http://climate.weatheroffice. gc.ca/prods_servs/index_e.html). Snow-depth measurements from 40 stations located in open areas were converted to binary snow/no snow data through the application of a 1 cm threshold. This dataset was used as a complementary validation of the AMSR-E SCA product. The main advantage of this dataset over the MODIS snow cover dataset is the continuous temporal coverage of snow measurements regardless of cloud cover, although it has been shown that SCD measured at these stations is shorter than that of its surroundings .

Blending RS-SCA datasets
Among the benefits of having two available MODIS instruments is the possibility to combine their data to potentially decrease the number of cloudy pixels. Different combination methods exist, which may yield different results since a given pixel for a given day may contain different information between the two MODIS products. Wang et al. (2009) use a priority approach on pixel information while Parajka and Blöschl (2008) use MODT data to substitute cloudy pixels in MODA. This last method assumes that no snowmelt or snowfall occurs within the time shift between local overpass times. We used two combination methods: (1) MODA data was used to replace cloudy pixels in MODT (henceforth TAC, Figure 2a) and (2) Wang et al.'s method (henceforth TACW) with the following priority scheme: snow > no snow > clouds ( Figure 2b). The preference of MODT over MODA in TAC was a matter of the reliability of its data. Reliability in this case refers to the frequency of performance gain of hydrological simulation when assimilating MODT compared with MODA. As will be seen in the Results section, MODT and MODA yield different results. Also, due to most of the Aqua MODIS band 6 detectors being non-functional, MODIS band 7 replaces band 6 in the snow-mapping algorithm, which has been shown to lead to a lower overall classification accuracy (Wang et al., 2009). Knowing this before the decision was taken on how to combine MODT and MODA, we decided to try one case where MODT is preferred over MODA (TAC) and one case with a different bias, where the information is prioritized instead (TACW). More on the reliability matter is presented in the Discussion section.
To further increase the number of pixels with surface snow information (snow/no snow), we blended the combined MODIS SCA products with AMSR-E SCA using a method similar to Gao et al. (2010). The cloudy pixels left in the combined MODIS SCA products were simply replaced with the corresponding AMSR-E SCA pixels to produce a cloud-free MODIS and AMSR-E combined SCA product (MAC SCA ). However, due to the sensitivity of AMSR-E SCA to melting snow, because the emissivity of wet snow drastically increases (Dupont et al., 2012), only AMSR-E SCA pixels reporting snow were used such that the no data (or 'cloudy') class was not entirely removed. MAC SCA created using TAC and TACW are henceforth referred as MAC and MACW, respectively. The upper part of the flowchart presented in Figure 3 shows the MAC blending.

Accuracy assessment
The MODIS products aggregated onto the hydrological model's 0.1˚× 0.1˚grid were compared over all three watersheds, for a total of 45 pixels, over the springs of 2004-2007. The standard error matrix with producer's and user's accuracies, as defined in Congalton (1991), was used for the relevant classes (snow and no snow). The omission and commission errors are the errors associated with the producer's and user's accuracies, respectively. The same accuracy assessment method was used to compute error

Data assimilation
Earth observation data assimilation applied to hydrology is not used to its full potential, mainly because of the lack of dedicated models designed to accommodate these data (Nagler et al., 2008). Various data assimilation methods exist to update time-dependent model states through observations, but their implementation in models is often not trivial and the computational effort can be very demanding (Moran et al., 2004). In an operational setting where 3-h forecasts are made and updated frequently, simplicity is an asset. The simplest data assimilation method and easiest to implement in a model is the direct insertion method, which assumes there is no error in measurements and that there is no correlation between measurements and the model states. The measured values simply replace the modelled values. However, since RS-SCA products contain only SCA data, which are not directly a modelled variable, we used the method in Roy et al. (2010) that is itself based on McGuire et al. (2006). Using this method, the hypothesis was that snow cover cannot be detected using remote sensing below a certain SWE threshold. This SWE threshold hypothesis is based on the decrease of the surface albedo through discontinuous snow cover and snow impurity observed in spring. The decrease in surface albedo can be enough to make the distinction from snow-free surfaces difficult. A similar assumption can be made that snow cover detected using remote sensing must contain at least as much SWE as the threshold. Using the direct insertion approach, if the modelled SWE was lower than the SWE threshold and RS-SCA detected snow, the modelled SWE was increased to match the SWE threshold. In the opposite case where the modelled SWE was greater than the threshold and RS-SCA detected a snow-free surface, the modelled SWE was decreased to the same amount. The model and RS-SCA were considered to be in agreement under the remaining scenarios and therefore no changes were made (Figure 3). In its present mode, this method does not preserve the water balance, as neither does the direct insertion of snow surveys.
A simple alternative to the direct insertion method is what will be referred to as the empirical method, which uses an empirical factor proportional to the difference between a modelled variable and the observed one to obtain a compromise between the two. In this approach, a weight is applied to the modelled states and observations (SWE threshold), where a 100% weight on observations corresponds to the direct insertion method and a 0% weight on observations corresponds to the modelled state without any update. This data assimilation method was applied to a few selected RS-SCA products to evaluate its performance.

MODIS SCA product intercomparison
Although there are two similar MODIS-based SCA products from the two MODIS instruments, they may differ due to differences in overpass time and the algorithm used to estimate snow presence.   2004 to 2008 for the three watersheds. The two products are in agreement over 50% of the time for all three classes (i.e. snow, no snow, and cloud). A significant difference in the cloud-covered pixels is noted, where MODA reports 75.4% of pixels to be cloud covered compared to MODT's 70.5%. The number of times that MODA reports clouds when MODT reports snow is over twice as high as when MODT reports no snow. The number of times that MODA reports snow when MODT reports no snow is, on average, over five times higher than when MODA reports no snow when MODT reports snow. Comparing the watersheds, this factor changes from 4.5 for the du Nord watershed, to 6.7 for the aux Écorces watershed, and 10 for the au Saumon watershed. These differences between the two SCA products were expected to lead to noticeable differences in streamflow estimates when assimilated into the HFS. Removing cloud-covered pixels from the analysis revealed agreement between MODT and MODA ranging from 86 to 97.2% for aggregated pixels reported as snow covered and 91.8 to 98.5% for snow-free surfaces, depending on whether MODT or MODA was considered as the reference.
The merging of the MODT and MODA images was carried out on a pixel by pixel basis. The merging reduced the percentage of cloud-covered pixels to 64.9%, providing 19.3% and 43.1% more usable pixels than MODT and MODA, respectively. As with their individual components, the combined products TAC and TACW also showed differences. While TAC had 2725 snowcovered pixel-days, TACW had 9% more snow due to the prioritization of snow in its algorithm.

Accuracy assessment of AMSR-E-based SCA
The spatial coverage and relatively high accuracy of MODT made it a good candidate on which to base the design of the AMSR-E-based SCA product. The ΔTB threshold used to create these maps was such that the total disagreement of AMSR-E SCA with MODT data was minimized. The optimal threshold calculated was found to be À2.5 K and it was applied in order to produce an AMSR-E SCA product, which was then compared with MODT, MODA, and ground-based EC station data (Table II). Some differences are noticed between the accuracies produced with the two MODIS SCA products that may be caused by the different overpass time of the satellites and the algorithms used to process MODT and MODA (see section 'MODIS'). Above all, an important difference is noticed in the accuracies (both user's and producer's) of the snow class, likely due to undersampling in MODIS/Aqua's case (162 snow pixels versus 681 for MODIS/Terra). Also, section 'MODIS SCA product intercomparison' already showed that MODIS/Aqua's snow class has lower user's accuracy when comparing it to MODIS/Terra's producer's accuracy. As for the comparison with EC stations, the overall agreement dropped to about 83%, but the underestimation of snow was the dominant source of disagreement, which was likely due to the drastic increase in snow emissivities at 19 and 37 GHz when it was wet.

Time series analysis of watershed-scale SCA
Results in the previous section showed good overall agreement between AMSR-E SCA and MODIS-based SCA products. However, the timing of the discrepancies can have an important impact on real-time hydrological simulation or forecast applications. An erroneous snow melt before the true snow melt can generate erroneous simulated streamflow peaks, which generate misleading information for dam managers. The time series analysis of SCA for the three watersheds ( Figure 4) shows AMSR-E SCA occasionally detected sudden snow-free surfaces for entire watersheds that are days and even weeks earlier than TAC. TAC often detected more gradual depletion curves as are generally seen for mid-latitude watersheds of several hundreds of square kilometres. Incorrect snowmelt detection can greatly affect the accuracy of streamflow simulation and forecast. However, using only AMSR-E SCA snow data, as opposed to the full AMSR-E SCA product including the nosnow information, in the blending algorithm eliminated these anomalies. As a consequence, this change greatly modified the contribution of AMSR-E SCA to the blended product such that the overall number of cloudy pixels was reduced to 40% instead of 0% since AMSR-E SCA often reported pixels to be snow free when TAC and TACW reported clouds.
Another interesting point to the time series analysis is that the snow cover transitions observed with the MODIS product occur on a shorter time scale for aux Écorces watershed. Hall and Riggs (2007) showed that MODIS snow cover mapping accuracy is lowest for evergreen forests because they obscure the underlying surface. Since the aux Écorces watershed is mainly composed of evergreen forests, the quick transitions between snowcovered and snow-free pixels observed with MODIS are likely caused by the faster transitions between snowcovered and snow-free canopies rather than the slower transition between snow-covered and snow-free surfaces underlying the canopies.

SWE threshold optimization
RS-SCA data were assimilated in the SPH-AV snow model through the direct insertion method described previously and the resulting water from the snowpack was introduced into the MOHYSE hydrological model used to simulate streamflows. The SWE thresholds used in the direct insertion of RS-SCA data resulted in different values of the Nash-Sutcliffe coefficient (Nash coefficient), and bias for various RS-SCA products and watersheds ( Figure 5). The simulations included only the spring melt event (25 March to 25 May) from 2004 to 2007 for the SWE threshold optimization period.
The analysis of the Nash coefficient and bias curves produced by assimilation of the same RS-SCA product, but different SWE thresholds, revealed many differences between watersheds. Many RS-SCA assimilation simulations performed better than the snow survey assimilation's Nash coefficient for the du Nord and aux Écorces watersheds. However, most did not provide much improvement compared to simulations with no updates, if at all, for the au Saumon watershed. The curves also varied between simulations assimilating MODT and MODA individually, which should be similar if the two products were equivalent. It is interesting to note that MODT was the only product that could produce Nash coefficients better or comparable to snow survey assimilations for all watersheds, while MODA yielded systematically better, and also higher, streamflow biases compared to MODT. MODIS-combined SCA products also yielded distinct curves. The SWE threshold optimum for the TACW data assimilation was systematically lower than TAC's while its streamflow bias was systematically higher than TAC's. This was consistent with its prioritization of snow over no snow in its algorithm, as if the lower threshold was, in a way, compensating for the more frequent snow detection. While the Nash coefficient for both products performed adequately for the du Nord and aux Écorces watersheds given a certain range on the SWE threshold spectrum, the opposite was true over the same range for the au Saumon watershed.
AMSR-E SCA data assimilation behaved differently, which was partly because only pixels reporting snow were used due to its sensitivity to snow melt. The Nash coefficient and streamflow bias coefficients were steadily improved as a function of the SWE threshold when compared to simulations with no updates, with the exception of the au Saumon watershed. However, these improvements remained marginal compared to assimilation of snow surveys or the other RS-SCA products. AMSR-E SCA 's small impact on streamflow simulations carries over to blended RS-SCA data assimilation since both MAC and MACW data assimilation changed negligibly from their unblended MODIS SCA component, namely TAC and TACW, respectively. It is important to note that while increasing the SWE threshold improved the streamflow bias in nearly all RS-SCA assimilation cases, it was not necessarily true for the Nash coefficient. Adding a sufficient amount of SWE did not translate into matching streamflow curves. This is clearly observed for the au Saumon watershed (Figure 5a and 5b). Figure 6, Figure 7, and Figure 8 show the hydrographs simulated by direct insertion of MAC using its optimal SWE threshold of 6 cm and how it compared with observations and simulations with no updates and with snow survey assimilation for the du Nord, aux Écorces, and au Saumon watersheds, respectively. Little variation is noted for the du Nord watershed with the exception of the year 2007, when the added snow from MAC resulted in a more accurate estimate of the main streamflow peak. For the aux Écorces and au Saumon watersheds, the variations between simulations were much greater. These variations were caused by disagreements between modelled SWE and remotely sensed snow observations via the SWE threshold. Some variations allowed for a more accurate simulation of streamflow, as shown for the aux Écorces watershed in 2005 (Figure 7). The model overestimated the SWE causing the peak streamflow to be overestimated for simulations without remote sensing, while MAC observed no snow causing the simulations with remote sensing to adjust the SWE. Other variations yielded less accurate streamflow simulations, as shown for the au Saumon watershed in 2006 (Figure 8). Around the first of May, MAC observed some snow covering the watershed causing the model SWE to be adjusted upward. This resulted in a significant increase in streamflow which was neither measured nor modelled without remote sensing. This was likely a result of an erroneous observation for the MAC product. Table III summarizes the comparison of streamflow simulation results obtained through RS-SCA data assimilation using SWE thresholds individually to those obtained through optimizing the Nash coefficient for each product. Nash coefficients show that assimilation of RS-SCA based on MODIS was comparable to biweekly snow survey assimilation overall and well above AMSR-E SCA data assimilation or simulations with no updates. In each watershed, streamflows were consistently underestimated, but the overall streamflow bias was improved in all cases using RS-SCA data assimilation with the exception of simulations using MODT, where the underestimation of streamflow for the aux Écorces watershed was important enough to result in the largest overall streamflow bias.

RS-SCA assimilation results
Another analysis evaluated the estimates of the main seasonal streamflow peaks. Table IV shows observations of the main seasonal streamflow peaks compared to simulation results obtained with RS-SCA and snow survey assimilation. The amplitude (difference) and timing (lag) of observed seasonal steamflow maxima were compared with the simulated ones. As with seasonal streamflows, maxi- mum streamflow amplitudes were underestimated overall. While snow survey assimilation yielded the best results for both the maximum difference and lag, the assimilation of MODT and TAC was not far behind and was better than no updates. AMSR-E SCA data assimilation yielded little to no improvement over simulations with no updates, and blended Although this last analysis was useful for evaluating seasonal maximum streamflow simulations, it did not consider the false identification of important peaks. A notable example of these false peaks was evident in the au Saumon watershed between 2004 and 2006 ( Figure 5). However, the Nash coefficients included the error of these false identifications.
Though only the direct insertion method was used so far, as a starting point to assimilate RS-SCA, the empirical method was also applied, but only to some promising candidates (TAC and TACW). We used the same SWE thresholds as obtained through direct insertion in order to compare this method to direct insertion. A weight of 50% was given to the RS-based updates through the SWE threshold, and the other 50% was kept for the modeled SWE. This weight was chosen such that the Nash coefficient was optimized for all three watersheds over the springs of 2004-2007 using discrete values of 0, 20, 50, 80, and 100%. Table III shows that this method essentially closes the gap between the results obtained through TAC and TACW assimilation. Though the streamflow bias was increased to some extent above the values obtained through direct insertion of the same RS-SCA products and SWE thresholds, the overall values remained lower than the simulations with no updates or snow survey assimilation. The mean Nash coefficients were increased above previous results and the variation across watersheds was reduced, although the Nash coefficient obtained with TAC remained somewhat below the Nash coefficient obtained with no updates and with snow surveys. As for peak streamflow estimation, using the empirical method resulted in noticeable bias differences in the overall amplitude difference and lag (Table IV). The RMSE was not significantly changed, however, with the exception of an improved peak difference for the assimilation of TACW and increased negative lag for the assimilation of TAC.

DISCUSSION
Intercomparisons of MODIS pixels aggregated on a 0.1˚× 0.1˚grid revealed that MODA reported cloudy pixels significantly more often than MODT on average. This was consistent with previous studies carried out with MODIS pixels that were not resampled (Parajka and Blöschl, 2008;Wang et al., 2009). This difference mainly occurred when MODT reported snow on the surface. Hall and Riggs (2007) explain a potential difference in snow mapping accuracy as a result of the substitution of the MODIS band 6 used in the snow-mapping algorithm by band 7 due to non-functioning detectors in band 6. The comparison also showed a difference in the type of disagreement between MODT and MODA pixels reported as snow or no snow. Using MODT as a reference, the commission error for snow (additional snow by MODA) was much higher than the omission error (additional snow by MODT), which may have been the cause of the systematically higher streamflow bias simulated when assimilating MODA compared to MODT. This also explains the lower optimal SWE threshold for MODA to compensate for more frequent snow detection. While the higher streamflow bias in MODA generally improved the streamflow bias, the Nash coefficient was generally improved with MODT, which suggested a larger uncertainty in the melting rate and timing of the added snow in MODA. Nonetheless, it was found that direct insertion of these products could significantly improve streamflow estimates compared to modelled simulations with no updates.
The idea of combining existing MODIS SCA products to remove cloudy pixels has been developed and applied in numerous studies (Parajka and Blöschl, 2008;Wang et al., 2009;Hall et al., 2010). Because of the differences in the two MODIS SCA products, the method used to combine these products can affect the end results. Indeed, direct insertion of TAC and TACW performed differently, thus yielding different results. Due to the prioritization of MODT in TAC, and snow, most often reported by MODA, in TACW, direct insertion of these two combined products has some similarities with MODT and MODA, respectively. Overall, direct insertion of these combined products yielded results comparable to the assimilation of biweekly snow surveys. The reduction of cloudy pixels in these combined products was likely the driving factor behind improvements in the estimates of streamflow peaks.
Estimated snow cover derived from AMSR-E brightness temperature differences was used to partly fill gaps in the data due to the remaining cloudy pixels. The snow cover mapping was in overall good agreement with cloud-free aggregated MODT and MODA, and in moderate agreement with EC stations. However, the overall underestimation of snow cover, likely due to its sensitivity to snow melt, significantly shortened the SCD. Choosing to validate AMSR-E SCA with other RS-SCA products, in this case MODT and MODA, was based mainly on spatial coverage. However, this validation method was limited by cloud cover. This is the reason why data from EC stations were also used. But, there is a scaling issue in comparing point data with 12.5 km × 12.5 km surfaces. The problem of the spatial representation of ground stations to validate RS-SCA or related data products has been examined in previous studies (Blöschl, 1999;Romanov et al., 2002;Simic et al., 2004;Chang et al., 2005). Ground measurements are sometimes treated as random samples within a region rather than as representative samples in order to bypass this issue (Brubaker et al., 2005); however, these measurements are generally located in open areas where the effect of the canopy and its projected shadow, which was particularly important for our study sites, on SCD and SWE is not represented (Musselman et al., 2008). Knowing that the stations generally underestimate snow cover in forested areas such as those used in this study  and that AMSR-E SCA underestimates snow cover compared with these stations, such that the errors add instead of cancelling, we can conclude that AMSR-E SCA underestimates snow cover.
Direct insertion of RS-SCA products yielded streamflows that were sometimes more and sometimes less accurate than simulations with no RS. The effect of these erroneous RS-SCA observations could be lowered if their error could be estimated and taken into account. Indeed, by simply adding a 50% weighting factor to updates resulted in better Nash coefficients overall using TAC and TACW data with the SWE thresholds optimized using direct insertion. However, it had a larger negative bias due to smaller increments of modelled SWE per update. The smaller variations of Nash coefficients between watersheds using either TAC or TACW demonstrated the method's flexibility in assimilating different SCA products for different areas. However, the increased lag in the streamflow peak estimates was to be expected since it can take several days before the updated modelled SWE is adjusted approximately to the SWE threshold due to the weighting factor, and the melting snow and cloud cover between updates. This lag was nonetheless shorter compared to using no updates at all.

CONCLUSION AND RECOMMENDATIONS
The aim of this study was to investigate assimilating the blended RS-SCA product and its components in a simple and easy way into the CEHQ's HFS for near real-time spring streamflow simulation applications. MODIS SCA products were first compared and combined and then blended with a SCA product based on the AMSR-E brightness temperature difference. The AMSR-E-based SCA product was validated using ground-based measurements and MODIS SCA products. Each individual, combined, and blended RS-SCA product was assimilated in the CEHQ's HFS through direct insertion, and the empirical method for a few selected products.
Although in overall good agreement with cloud-free optical-and NIR-based RS-SCA products, the SCA product based on passive microwave data was not found to significantly improve streamflow estimates. The increase in SCA data availability from AMSR-E SCA data resulted in a minor improvement of the peak and overall streamflow simulation when compared with simulations with no updates, but it did not show an improvement over the assimilation of snow surveys or MODIS SCA data. This effect was carried over to the combined SCA products, where results resembled MODIS-only products such that the contribution from AMSR-E SCA data was essentially ineffective.
Though direct insertion of SCA derived from AMSR-E data may possibly have a stronger impact for larger watersheds using the proposed method, the study of the scale at which SCA derived from passive microwave becomes significantly useful, if at all, was not part in this study. However, it may be an interesting avenue for future works.
Direct insertion of the combined MODIS products yielded results comparable to the results using direct insertion of biweekly snow surveys for overall streamflow simulation, but the latter was more exact for the main seasonal streamflow peaks. Another suggestion for future works would be to study the potential synergy between in situ and remote sensing data pertaining to snow for operational hydrological simulation purposes.
The comparison of simulated streamflow through assimilation of various RS-SCA products was mainly based on the analysis of the Nash-Sutcliffe coefficient. However, the significance of this coefficient, which varies from one year to another and between watersheds, is difficult to define precisely because it appears to be quite variable. Overall, the Nash criterion for the simulations with assimilation of RS-SCA data was similar to the one using snow surveys assimilation for the cases analysed. However, given the variability between basins, it is difficult to generalize this conclusion over other areas. Therefore, the 6 cm threshold found to be optimal in this study is unlikely to be optimal for other regions, but the method may be applied elsewhere to obtain such a threshold.
Finally, the empirical method was found to be better suited than direct insertion to assimilate RS-SCA products in the HFS. This result suggests that measurement errors must be taken into account. Incorporating more sophisticated assimilation methods, such as the Ensemble Kalman filter (Andreadis and Lettenmaier, 2006) or using a method that conserves mass balance during updates , could improve the precision of streamflow simulations.