A Hydrograph-Based Approach to Improve Satellite-Derived Snow Water Equivalent at the Watershed Scale

: For the past few decades, remote sensing has been a valuable tool for deriving global information on snow water equivalent (SWE), where products derived from space-borne passive microwave radiometers are favoured as they respond to snow depth, an important component of SWE. GlobSnow, a novel SWE product, has increased the accuracy of global-scale SWE estimates by combining remotely sensed radiometric data with other physiographic characteristics, such as snow depth, as quantiﬁed by climatic stations. However, research has demonstrated that passive microwaves algorithms tend to underestimate SWE for deep snowpack. Approaches were proposed to correct for such underestimation; however, they are computer intensive and complex to implement at the watershed scale. In this study, SWE max information from the near real time 5-km GlobSnow product, provided by Copernicus and the European Space Agency (ESA) and GlobSnow product at 25 km resolution were corrected using a simple bias correction approach for watershed scale applications. This method, referred to as the Watershed Scale Correction (WSC) approach, estimates the bias based on the direct runoff that occurs during the spring melt season. Direct runoff is estimated on the one hand from SWE max information as main input. Inﬁltration is also considered in computing direct runoff. An independent estimation of direct runoff from gauged stations is also performed. Discrepancy between these estimates allows for estimating the bias correction factor. This approach is advantageous as it exploits data that commonly exists i.e., ﬂow at gauged stations and remotely sensed/reanalysis data such as snow cover and precipitation. The WSC approach was applied to watersheds located in Eastern Canada. It was found that the average bias moved from 33.5% with existing GlobSnow product to 18% with the corrected product, using the recommended recursive ﬁlter coefﬁcient β of 0.925 for baseﬂow separation. Results show the usefulness of integrating direct runoff for bias correction of existing GlobSnow product at the watershed scale. In addition, potential beneﬁts are offered using the recursive ﬁlter approach for baseﬂow separation of watersheds with limited in situ SWE measurements, to further reduce overall uncertainties and bias. The WSC approach should be appealing for poorly monitored watersheds where SWE measurements are critical for hydropower production and where snowmelt can pose serious ﬂood-related damages.


Introduction
Hydrological processes for watersheds located in northern latitudes and at high elevations are greatly affected by seasonal snowfall. Spring snowmelt is arguably one of the most important processes for operational hydrology of these territories, as it is a critical component for determining flood risk [1]. However, the prediction of spring snowmelt is contingent on having accurate estimates of the amount of liquid water stored in the snowpack, i.e., snow water equivalent (SWE), and more particularly, the maximum snow water equivalent (SWE max ) before the onset of snowmelt. This is especially true for effective management of water resources and for predicting spring floods over large watersheds where spring snowmelt may last for extended periods and where there is a delay between the occurrence of SWE max and peak flow. Additionally, accurate estimates of SWE are proposed to reduce the bias with some success, such as explicitly considering snowpack stratification in the microwave emission model [11,14] and assimilating PM data into a microwave emission model using snow grain size simulated by a physically based snowpack model [14]. However, these approaches are complex to implement operationally and, therefore, there is a need for a simple bias correction approach that can provide SWE and SWE max estimates at the watershed scale. The Watershed Scale Correction (WSC) approach proposed here aims to correct the underestimation of spatially varying SWE, including SWE max , produced by GlobSnow, at the watershed scale. This approach innovatively applies historical streamflow measurements and, thus, can be applied even if there are no available ground-based snow measurements. To do so, runoff water volumes retrieved from streamflow gauge stations are compared with the volumes resulting from SWE max of GlobSnow, and then corrected based on the volumetric difference. The functionality of the approach is demonstrated by comparing the corrected SWE values with those obtained from snow course measurements at selected watersheds of varying sizes in the province of Québec, Canada.

WSC Overview
The WSC approach aims to improve the accuracy of the SWE max of regional databases, such as GlobSnow V3 [11] using streamflow information. The approach was coded in MATLAB.
The motivation behind the WSC approach is that the watershed surface runoff originating from snowmelt, along with any precipitation during the spring snowmelt period, will eventually leave the watershed through its outlet. This approach focuses on direct runoff (i.e., total runoff minus baseflow) instead of looking at total runoff. According to the mass balance equation applied at the watershed scale, the direct runoff calculated from the amount of snow over the watershed, referred to as DRG (Direct Runoff using GlobSnow), should be equal to DRH, the direct runoff hydrograph retrieved from measured streamflow, which is: However, given that SWE max from GlobSnow is typically underestimated, it is likely that DRG be less than DRH. To remove this bias, an adjustment is applied to GlobSnow's SWE max , according to the following equation: where, SWE Glob max is the SWE max estimated from GlobSnow, SWE corr max is the corrected SWE max , and CF is the correction factor value. Note that a single CF value is applicable to the entire watershed, while SWE max (both corrected and GlobSnow) is spatially distributed.
The estimation of DRG and DRH can be done using the following equations: where, RO total is the total runoff estimated from the streamflow hydrograph, B is the baseflow, I is the infiltration, P is the total precipitation occurring during the snowmelt period, and A is the watershed surface area.  (5) Given that SWE Glob max is typically underestimated, DRG is likely less than DRH, and therefore CF will be larger than 1. Recall that a single CF value is applicable to the entire watershed, while SWE max (either corrected or GlobSnow) is spatially distributed. A CF value should be computed for each year, where a spring runoff hydrograph is available. Therefore, it is expected that no single CF value will characterize a watershed for the following reasons: 1.
SWE Glob max is derived using the HUT model. GlobSnow V1 and V2 [12] assume that the snowpack is made of a single homogeneous layer, while V3 [11] represents the snowpack as a stacked system of snow layers. Despite this improvement over the original HUT model, uncertainties in SWE max estimation remain, especially when the snowpack undergoes multiple freeze-thaw and rain-on snow events, which will be reflected by changing CF values.

2.
Because of the limited penetration depth of the PM signal, it is expected that the CF will be larger when deeper snowpack conditions are encountered. 3.
Flow measurements also have errors related to the rating curve used to convert stream water level into streamflow. 4.
Finally, other sources of uncertainty related to precipitation, baseflow, and infiltration estimates will also affect CF.
Therefore, the proposed WSC approach is likely to provide a range of CF values. Although this might raise questions regarding how to select the 'best' CF value for correcting SWE, when the approach is applied to flood forecasting, the range of CF values is compatible and may be seen as reflecting the various sources of uncertainties affecting SWE estimation.
The WSC method can be divided into three main blocks as depicted in Figure 1. The first block relates to the organisation and retrieval of data, based on the watershed's location. Each watershed is characterized by its own hydrometeorological regime and state variables (e.g., streamflow, precipitation, temperature, and SWE).
Given that SWE max Glob is typically underestimated, DRG is likely less than DRH therefore CF will be larger than 1.
Recall that a single CF value is applicable to the entire watershed, while SWEm ther corrected or GlobSnow) is spatially distributed. A CF value should be compute each year, where a spring runoff hydrograph is available. Therefore, it is expected th single CF value will characterize a watershed for the following reasons: Glob is derived using the HUT model. GlobSnow V1 and V2 [12] assume tha snowpack is made of a single homogeneous layer, while V3 [11] represents the s pack as a stacked system of snow layers. Despite this improvement over the ori HUT model, uncertainties in SWEmax estimation remain, especially when the s pack undergoes multiple freeze-thaw and rain-on snow events, which will b flected by changing CF values. 2. Because of the limited penetration depth of the PM signal, it is expected that th will be larger when deeper snowpack conditions are encountered. 3. Flow measurements also have errors related to the rating curve used to co stream water level into streamflow. 4. Finally, other sources of uncertainty related to precipitation, baseflow, and inf tion estimates will also affect CF.
Therefore, the proposed WSC approach is likely to provide a range of CF va Although this might raise questions regarding how to select the 'best' CF value for recting SWE, when the approach is applied to flood forecasting, the range of CF valu compatible and may be seen as reflecting the various sources of uncertainties affe SWE estimation.
The WSC method can be divided into three main blocks as depicted in Figure 1 first block relates to the organisation and retrieval of data, based on the watershed's tion. Each watershed is characterized by its own hydrometeorological regime and variables (e.g., streamflow, precipitation, temperature, and SWE). The second block uses the information obtained in Block 1 to estimate watershed variables and hydrological processes required to compute DRG and DRH, and ultim the correction factor (as shown in Block 3). Watershed state variables include SWEmax The second block uses the information obtained in Block 1 to estimate watershed state variables and hydrological processes required to compute DRG and DRH, and ultimately the correction factor (as shown in Block 3). Watershed state variables include SWE max from the GlobSnow database and the duration of the snowmelt period, defined as the difference between the time SWE max occurs and when the snow has disappeared. Both state variables can be retrieved from time series of SWE provided by the GlobSnow database. Alternatively, snow cover area obtained from ERA5 land cover product, based on information from MODIS datasets, can be used to establish the last day where there is snow in a watershed, and therefore to estimate the duration of the snowmelt period. Watershed hydrological processes include precipitation during the snowmelt period, watershed infiltration and baseflow. Once the duration of the snowmelt period is established, the total precipitation (rain and snow) occurring during this period can be retrieved. Baseflow separation and watershed infiltration methods used here are briefly described below.
Baseflow separation. Baseflow separation was performed using the HydRun tool box [15], a nested group of functions used for rapid and flexible rainfall-runoff analysis. The package, written in MATLAB, allows computing baseflow (B t ) and direct runoff (q t ) from total runoff Q t , see Equations (6) and (7) below. The baseflow is computed using the recursive digital filter technique proposed by Nathan and McMahon [16]: where B t is the baseflow at time t, q t is the direct runoff, Q t represents the original streamflow and β, is a filter coefficient which affects the magnitude of the baseflow values. For example, increasing β produce a lower baseflow. Researchers used prescribed values for β to obtain optimal results [17,18]. Gan and Zuo [17] set the β value between the range 0.90-0.95, whereas Mau and Winter [18] used 0.85 for best results. Nathan and McMahon [16] suggested that the optimal β value should be 0.925. It is advised to thoroughly investigate the β value depending on the hydro-geological conditions in the catchment. To facilitate baseflow separation, the software also offers options to preprocess the hydrograph, such as hydrograph filtering establishing the maximum tolerable difference of discharge between the start and finish of an event and discarding peaks lower than a threshold value, allowing small fluctuations to be ignored while computing B t . The toolbox is effective for analysis of rainfall-generated hydrographs but presents some uncertainty as it remains to be tested for hydrographs generated by snowmelt [15]. However, other studies have applied the digital recursive method of separation to areas with snow cover [17,[19][20][21]. While Hammond and Kampf [19] and Fan et al. [21] used β = 0.925 as the optimal filtering parameter value as suggested by previous studies [16,20], baseflow separation is highly sensitive to the filtering parameter and heuristic baseflow separation should be done with care, as true values of baseflow in snow-dominated systems are seldom available for optimizing the filtering parameter value [20].
As an alternative to observed baseflow values, Ref. [17] simulated baseflow in glacier melt dominated basins using the SWAT hydrological model and compared these 'observations' with baseflow produced according to the recursive filter method with β = 0.925. They found that the recursive filtering method produced flows comparable to SWAT during the low-flow period but overestimated baseflow in the high-flow method. Therefore, for the catchments examined here we chose a β range of values from 0.800-0.995, with one run including the suggested optimal value, as it will be seen/discussed later.
Infiltration. Frozen ground is common to northern regions and infiltration of meltwater into frozen ground is a complex process that involves exchange of heat and mass flow with phase change [22]. To estimate DRG (see Equation (4)), it is required to estimate the amount of water infiltrating into frozen ground. According to Nicholaichuk and Gray [23], infiltration of snowmelt into frozen ground involves three possible water transfer and storage conditions depending on the surface entry conditions for the study area: restricted, unlimited, and limited. Restricted and unlimited water flows result in most all, or no snowmelt being infiltrated, respectively. In contrast, limited water flows result in some snowmelt being infiltrated and is influenced primarily by the soil physical properties, such as the available water storage capacity. Gray et al. [22] describe an algorithm that calculates areal infiltration into frozen ground on a watershed or hydrological response unit that combines the above three infiltration conditions. For each landscape unit in the watershed, the fraction classified as unlimited, restricted and limited infiltration varies dynamically according to the water holding capacity of frozen soil and cumulative infiltration is computed accordingly. The equation used to quantify the Infiltration for frozen, unsaturated soils of limited infiltrability is described as: where INF if the cumulative infiltration (mm), S 0 is the soil moisture content at the soil surface, C is the land cover type coefficient, t 0 is the infiltration opportunity time (h), S I and T I are the average soil saturation (mm 3 /mm 3 ) and temperature (K), respectively, of the first 40-cm soil layer at the start of infiltration. S 0 is usually taken as 1 due to low infiltration rate for frozen soil in areas where snow ablates rapidly. For prairie land covers C is 2.10, while for the forest covers it is 1.14.
The variable t 0 is the time it takes for the snow cover to completely melt assuming that the melt is continuous, there is a small storage space, and little evaporation [22]. The following equation was used to estimate the t 0 : Here, we assumed that limited infiltration conditions prevail during the snowmelt period and, accordingly, Equation (8) was used to calculate total infiltration into frozen soil. Note that depending on the value taken by S I , which is the soil saturation at the onset of the snowmelt, solving Equation (8) would result in restricted infiltrability if the soil is completely saturated (S I = 1, therefore INF = 0) and unlimited infiltrability would result as S I = 0. Further information on the infiltration into frozen ground algorithm is given in Gray et al. [22] and in Zhao and Gray [24]. Figure 2 proposes a linear relationship between CF and SWE. As mentioned above, CF corrects for the bias in SWE Glob max estimates based on the direct runoff that occurs during the spring melt season. The same correction factor cannot be applied to correct SWE Glob observations during the winter season, which presents SWE information when SWE is below SWE max . Previous research e.g., [25][26][27] has shown that PM starts underestimating SWE under deep snow conditions. It is therefore conceivable that the CF factor should vary between 1 (no correction) when SWE is above a certain threshold value, called SWE thres , and CF max , the correction factor corresponding to SWE max . some snowmelt being infiltrated and is influenced primarily by the soil physical properties, such as the available water storage capacity. Gray et al. [22] describe an algorithm that calculates areal infiltration into frozen ground on a watershed or hydrological response unit that combines the above three infiltration conditions. For each landscape unit in the watershed, the fraction classified as unlimited, restricted and limited infiltration varies dynamically according to the water holding capacity of frozen soil and cumulative infiltration is computed accordingly.
The equation used to quantify the Infiltration for frozen, unsaturated soils of limited infiltrability is described as: where INF if the cumulative infiltration (mm), So is the soil moisture content at the soil surface, C is the land cover type coefficient, t0 is the infiltration opportunity time (h), SI and TI are the average soil saturation (mm 3 /mm 3 ) and temperature (K), respectively, of the first 40-cm soil layer at the start of infiltration. So is usually taken as 1 due to low infiltration rate for frozen soil in areas where snow ablates rapidly. For prairie land covers C is 2.10, while for the forest covers it is 1.14.
The variable t0 is the time it takes for the snow cover to completely melt assuming that the melt is continuous, there is a small storage space, and little evaporation [22]. The following equation was used to estimate the t0: Here, we assumed that limited infiltration conditions prevail during the snowmelt period and, accordingly, Equation (8) was used to calculate total infiltration into frozen soil. Note that depending on the value taken by SI, which is the soil saturation at the onset of the snowmelt, solving Equation (8) would result in restricted infiltrability if the soil is completely saturated (SI = 1, therefore INF = 0) and unlimited infiltrability would result as SI = 0. Further information on the infiltration into frozen ground algorithm is given in Gray et al., [22] and in Zhao and Gray [24]. Figure 2 proposes a linear relationship between CF and SWE. As mentioned above, CF corrects for the bias in SWE max Glob estimates based on the direct runoff that occurs during the spring melt season. The same correction factor cannot be applied to correct SWE Glob observations during the winter season, which presents SWE information when SWE is below SWEmax. Previous research e.g., [25][26][27] has shown that PM starts underestimating SWE under deep snow conditions. It is therefore conceivable that the CF factor should vary between 1 (no correction) when SWE is above a certain threshold value, called SWEthres, and CFmax, the correction factor corresponding to SWEmax.  Whether or not this linear relationship holds true would require a detailed analysis of GlobSnow derived SWE against ground measurements. Nevertheless, for simplicity a linear relationship was applied here and provided reasonable estimates of corrected SWE values when compared against observed ground measurements (see result section for details). Here, we calculated CF only during the build-up of the snowpack. Initially, we defined SWE thresh as 150 mm following Luojus et al. [27] However, we observed that some years GlobSnow reported values of SWEmax below 150 mm and our associated CF values were smaller than 1 suggesting that SWE thresh should be smaller than 150 mm. This is in line with Larue et al. [7] in their analysis of GlobSnow data in Eastern Canada, who observed that PM SWE starts to be underestimated at depths above 100 mm. Accordingly, we defined SWE thresh value as 100 mm.

Study Area
The study was conducted on eight watersheds located in south central Quebec, Canada, between 50 to 54 • N latitude and 67 to 70 • W longitude, see Figure 3. All watersheds have natural flows, except for the Manic 5 watershed, which has regulated flows. The watersheds' areas vary from 795 km 2 (02PG022-Ouelle River Watershed) to 24,698 km 2 (Manic 5 watershed), see Table 1. All watersheds are mostly forested; for example, forests in 02PA007 (Batiscan River) and Manic 5 watersheds occupy respectively 87 and 83% of the total surface area. The remaining landcover is characterized as either open water (including wetlands), agricultural, or urban areas; for example, 14% of Manic 5 watershed consists of open water and wetlands, while only 7% of 02PA007 watershed is open water and wetlands. The Manic 5, 02RH035 (Aux Écorces River), 02NE011 (Croche River), and 02PA007 watersheds are in the Canadian Shield, which is characterized by small lakes and rolling hills as the result of glacial erosion from the last ice age. 02PL005 (Upper Bécancour River), 01010000 (St-John River), 01AD003 (St-Francis), and 02PG022 watersheds are located in the Appalachian Region of Eastern Canada, which is characterized by rolling hills and mountains. These watersheds are characterized by heavy winter snowfall and snowmelt runoff in the spring. Average maximum snow water equivalent (SWE) varies from 217 mm for 02PL005 watershed up to 298 mm for 02RH035 watershed, and generally increases from south to north because of lower air temperature with increasing latitude. Annual peak flow predominantly occurs during the spring thaw. The Manic 5 watershed, the largest in size and with heavy snowpack, has mean annual discharge is 529 m 3 /s [28] and annual peak flows exceeding on average 2000 m 3 /s, occurring in May. SWE can reach values of 300 mm or more [29]. Flow from Manic 5 drains into a 2000 km 2 hydropower reservoir, with an installed capacity of 2.6 GW [30].
Hydrometeorological information about these watersheds can be found in Table 1.  Hydrometeorological information about these watersheds can be found in Table 1.

GlobSnow
The GlobSnow V3 SWE dataset, as described in Luojus et al. [11], was selected for analysis. The product combines passive microwave radiometer data from Nimbus-7 Scanning Multichannel Microwave Radiometer (SMMR), the Special Sensor Microwave Imager/Imager (SSM/I) and the Special Sensor Microwave Imager/Sounder (SSMIS) Defense Meteorological Satellite Program (DMSP), along with ground based synoptic snow depth observations using Bayesian data assimilation in the HUT snow emission model [11]. The algorithm first solves the HUT model using weather station measurements of snow depth to obtain effective snow grain size by fitting simulated to observed satellite brightness temperature measurements at the locations of weather stations. This information is then interpolated to a 25 km × 25 km grid cell and used in an assimilation procedure for the estimation of SWE [10,27]. Takala et al. in their study [32] applied a two-dimensional convolution window of 0.25 • × 0.25 • (25 km × 25 km) to change the calculation grid to 0.05 • (approximately 5 km) projected in EASE-grid. The 25-km gridded GlobSnow product is available from 1979 to present. The 5-km product is available from 2006 to present and provides daily SWE 12 h after global satellite data has been acquired [33]. According to Luojus et al. [11], the 25-km product has an overall root mean square error (RMSE) for shallow to moderate snowpack (SWE < 150 mm) of 32.7 mm for the period covering 1980-2016 with no apparent bias but uncertainty increases when SWE is above this threshold.
We retrieved the SWE information for the study watersheds from the GlobSnow database. Figure 4 illustrates the temporal evolution of SWE for three watersheds, namely 02RH035 (Rivière des Écorces), 01AD003 (St-Francis River), and Manic 5. The temporal coverage, average SWE max from GlobSnow and from ground observations for the watersheds illustrated in Figure 4 along with other study sites are presented in Table 2, showing that GlobSnow underestimates SWE compared to in situ measurements.
interpolated to a 25 km × 25 km grid cell and used in an assimilation procedure for the estimation of SWE [10,27]. Takala et al. in their study [31] applied a two-dimensional convolution window of 0.25° × 0.25° (25 km × 25 km) to change the calculation grid to 0.05° (approximately 5 km) projected in EASE-grid. The 25-km gridded GlobSnow product is available from 1979 to present. The 5-km product is available from 2006 to present and provides daily SWE 12 h after global satellite data has been acquired [32]. According to Luojus et al. [11], the 25-km product has an overall root mean square error (RMSE) for shallow to moderate snowpack (SWE < 150 mm) of 32.7 mm for the period covering 1980-2016 with no apparent bias but uncertainty increases when SWE is above this threshold.
We retrieved the SWE information for the study watersheds from the GlobSnow database. Figure 4 illustrates the temporal evolution of SWE for three watersheds, namely 02RH035 (Rivière des Écorces), 01AD003 (St-Francis River), and Manic 5. The temporal coverage, average SWEmax from GlobSnow and from ground observations for the watersheds illustrated in Figure 4 along with other study sites are presented in Table 2, showing that GlobSnow underestimates SWE compared to in situ measurements.     The fifth-generation European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric Reanalysis of global climate (ERA5) database provides hourly data estimates at a resolution of approximately 31 km for the global atmosphere, land surface, and ocean waves from 1950 to present. We used ERA5 to obtain the average soil saturation (S I ) and temperature (T I ) required as input for the Gray's infiltration model. Compared to its predecessor ERA-interim, ERA5 has been greatly improved [34]. Improvements of the land surface model of ERA5 include a point wise simplified extended Kalman filter for three soil moisture layers in the top 1 m of soil, and a one-dimensional optimal interpolation for soil temperature [35]. Additional improvements of ERA5 are explained in detail by Hersbach et al. [36].

Streamflow Data
We used streamflow information from the Hydrometeorological Sandbox-École de technologie supérieure (HYSETS), a comprehensive dataset containing daily hydrometeorological data for over 14,400 watersheds across north America [31]. Streamflow information includes station location, flow regime (i.e., natural or regulated), daily flow data, as well as maximum and minimum daily air temperature and precipitation weather gauges, along with the SWE of ERA5-land. The database includes data covering the period 1950-2018 depending on the type and source of data. The streamflow data in HYSETS were retrieved from national water agencies repositories. For our study watersheds, streamflow is collected by Water Survey Canada (WSC) and stored in the National Water Data Archive (HYDAT). Further details on the HYSETS configuration are given in Arsenault et al. [31].

In Situ Data
The in situ SWE data, which from hence forth will be used interchangeably with the term 'reference SWE data', comes from Hydro-Québec's (HQ) snow course measurements for Manic 5 watershed and from Québec's Ministère de l'Environnement et de la Lutte contre les changements climatiques (MELCC) for the remaining watersheds investigated in this study. The SWE measurements are the average value of 10 samples equally spaced along a 300-m snow course [29]. These measurements are usually taken biweekly once the snowpack is developed until the snowmelt begins. The location of the measurement sites is depicted in Figure 3.

Results
Correction factors (CF) values were computed for all watersheds and years under study according to Equation (5). This was done through estimation of the direct runoff DRH, calculated with Equation (3) by removal of the estimated baseflow from total runoff estimated from gauged stations. Direct runoff was also estimated from SWE max data from GlobSnow, that is DRG, and considering precipitation and watershed infiltration occurring during the snowmelt period. Results are shown in Figure 5 and Table 3.

Results
Correction factors (CF) values were computed for all watersheds and years under study according to Equation (5). This was done through estimation of the direct runoff DRH, calculated with Equation (3) by removal of the estimated baseflow from total runoff estimated from gauged stations. Direct runoff was also estimated from SWEmax data from GlobSnow, that is DRG, and considering precipitation and watershed infiltration occurring during the snowmelt period. Results are shown in Figure 5 and Table 3.  As expected, CF values of a given watershed vary from year-to-year in response to varying snowpack characteristics (e.g., stratigraphy, grain size), which are much simpli-  As expected, CF values of a given watershed vary from year-to-year in response to varying snowpack characteristics (e.g., stratigraphy, grain size), which are much simplified in the radiative transfer model (HUT) used to retrieve GlobSnow's SWE estimates. Any important deviation from this snowpack representation will affect SWE estimation; for example, the presence of ice lenses in the snowpack may have a significant effect on microwave scattering and emission properties of the pack. This in turn may be wrongly interpreted as a snowpack with smaller SWE than reality, producing a larger CF value (CF > 1) (Equation (5)). Therefore, year-to-year variations in CF values can be seen as a consequence of varying hydrometeorological conditions leading to a heterogeneous snowpack, which are not adequately simulated by the HUT model used to estimate SWE. Saturation of the microwave signal as the snowpack gets thicker will also increase deviation of the estimated SWE compared to the 'real' SWE and, therefore, it is expected that large CF values will be obtained for years with thicker snowpack.
The average CF factor of the watersheds, as show in Table 3, has values ranging from 1.29 to 1.91. These results indicate that the range is rather small. The standard deviation for the study ranged from 0.29 to 0.53 with higher deviations coming from watersheds furthest south of the St. Lawrence, which may be an indication of difficulty in retrieving SWE data from the satellite microwave signal. A possible cause for such difficulty is the presence of complex snowpack stratigraphy, e.g., ice lenses, resulting from mid-winter rainon-snow and freeze-thaw events, which are common in these regions. This is in contrast to watersheds located further north, where such meteorological events are less susceptible to occurs and therefore where snowpack is generally more homogeneous. Consequently, standard deviation of the retrieved CF factors is generally smaller. Figure 6 shows examples of scatterplots of available in situ maximum SWE data against the original and corrected GlobSnow products with the recommended baseflow separation factor β = 0.925 for the 02RH035, 01AD003, and Manic5 watersheds.
The limited in situ measurements affected the analysis; for example, the Manic5 watershed only had 11 years of in situ measurements over the 18 possible years. The figure confirms that the original GlobSnow product underestimates SWE max ; for example, the observed SWE max values of the 02RH035 watershed vary between 170 and 450 mm but the corresponding GlobSnow values range from 120 to 270 mm. After correction, using the recommended β of 0.925, the range increases from 130 to 355 mm. Similar behavior are observed for all watersheds under study. Despite the bias being reduced, an underestimation remains.
Predominantly, the WSC approach increased the estimation of SWE max , as summarized in Table 4. The WSC approach increased the range between minimal and maximal SWE max values, and reduced the bias when compared to the referenced SWE max . However, the degree of linear correlation between estimated and reference SWE max values did not increase for all watersheds, nor did the RMSE systematically decrease (for example, see 02PL005). An increase of RMSE means that the WSC approach increased the dispersion of the corrected SWE max values with relation to observed SWE max . This is related to uncertainties along the modelling chain leading to the corrected SWE max product. For example, a recommended value of the baseflow separation factor β is 0.925 for all watersheds, but studies have shown that β can differ depending on the watershed's characteristics [18]. Higher/lower β values will decrease/increase baseflow and, therefore, directly affect the magnitude of direct runoff (DRH) and the magnitude of the corresponding correction factor (CF). The issue of model uncertainty is addressed later in this paper.
susceptible to occurs and therefore where snowpack is generally more homogeneous. Consequently, standard deviation of the retrieved CF factors is generally smaller. Figure 6 shows examples of scatterplots of available in situ maximum SWE data against the original and corrected GlobSnow products with the recommended baseflow separation factor β = 0.925 for the 02RH035, 01AD003, and Manic5 watersheds. The limited in situ measurements affected the analysis; for example, the Manic5 watershed only had 11 years of in situ measurements over the 18 possible years. The figure confirms that the original GlobSnow product underestimates SWEmax; for example, the observed SWEmax values of the 02RH035 watershed vary between 170 and 450 mm but the corresponding GlobSnow values range from 120 to 270 mm. After correction, using the recommended β of 0.925, the range increases from 130 to 355 mm. Similar behavior are observed for all watersheds under study. Despite the bias being reduced, an underestimation remains.
Predominantly, the WSC approach increased the estimation of SWEmax, as summarized in Table 4. The WSC approach increased the range between minimal and maximal SWEmax values, and reduced the bias when compared to the referenced SWEmax. However, the degree of linear correlation between estimated and reference SWEmax values did not increase for all watersheds, nor did the RMSE systematically decrease (for example, see 02PL005). An increase of RMSE means that the WSC approach increased the dispersion of the corrected SWEmax values with relation to observed SWEmax. This is related to uncertainties along the modelling chain leading to the corrected SWEmax product. For example, a recommended value of the baseflow separation factor β is 0.925 for all watersheds, but  The time series of the corrected GlobSnow product and the original GlobSnow against in situ SWE is shown in Figure 7. Once the snowpack is developed, the corrected SWE estimation becomes greater than the original GlobSnow estimation and approaches the in situ SWE data. However, the corrected GlobSnow still underestimates SWE for some years and overestimates for others. Overestimation occurs in years that have gaps in the in situ SWE data around the time when GlobSnow estimates the snowpack to be at its peak. The time series of the corrected GlobSnow product and the original GlobSnow against in situ SWE is shown in Figure 7. Once the snowpack is developed, the corrected SWE estimation becomes greater than the original GlobSnow estimation and approaches the in situ SWE data. However, the corrected GlobSnow still underestimates SWE for some years and overestimates for others. Overestimation occurs in years that have gaps in the in situ SWE data around the time when GlobSnow estimates the snowpack to be at its peak.

Discussion
In this study, we proposed a hydrograph-based approach to correct the underestimation of SWE derived from remotely sensed passive microwave data (PM), specifically, the GlobSnow product. The chief advantage of our Watershed Scale Correction (WSC) approach is that it uses data that commonly already exists (i.e., flow at gauged stations and remotely sensed data) and that it does not require labour intensive and costly field measurements of SWE. This contrasts with more 'classical' approaches that consist of

Discussion
In this study, we proposed a hydrograph-based approach to correct the underestimation of SWE derived from remotely sensed passive microwave data (PM), specifically, the GlobSnow product. The chief advantage of our Watershed Scale Correction (WSC) approach is that it uses data that commonly already exists (i.e., flow at gauged stations and remotely sensed data) and that it does not require labour intensive and costly field measurements of SWE. This contrasts with more 'classical' approaches that consist of merging point SWE measurements to SWE maps derived from PM. Moreover, the WSC approach preserves GlobSnow's physically based approach, but additionally incorporates a scaling factor derived from hydrological processes (i.e., direct runoff). While purely statistical approaches, such as spatial interpolation (e.g., simple approaches such as Thiessen polygons to more sophisticated methods such as kriging) between point SWE measurements have no physical basis. Additionally, the WSC approach allows for the correction of daily SWE products, such as GlobSnow, at the watershed scale, which SWE surveys cannot achieve. Therefore, the WSC approach reduces the need to have frequent field surveys to have greater temporal information on the state of watersheds. Therefore, the WSC approach adds value in the continued interest of presenting watershed state information on SWE during the snowpack accumulation period of winter. This in turn is valuable information for operational water resources managers to mitigate the threats of spring flood and also for optimal reservoir management for hydroelectric production.
Because PM typically underestimates SWE for deeper snowpacks [7], a number of physically based approaches have been proposed to correct the underestimation and have been successful. However, such approaches are complex to implement and are computer intensive when applied at the watershed scale. For example, LaRue et al. [14] implemented a chain of models to assimilate PM observations and obtained substantial bias reduction compared to original SWE simulations [14]. In their approach, they coupled the physically based Crocus snow model, driven by a global-scale atmospheric model, to the Dense Media Radiative Transfer-Multilayers (DMRT-ML) model to calibrate the snow stickiness parameter of DMRT-ML. They then assimilated AMSR-E T B observations in the calibrated DMRT-ML model to improve SWE estimations for 12 weather stations in northeastern Canada where continuous SWE measurements were available. Our WSC approach is comparatively simpler and can be applied to correct SWE at the watershed scale while retaining the physics of the algorithm used to convert T B into SWE.
Correction factor-based approaches were proposed in the literature and have been shown to provide more accurate SWE estimations. However, these approaches were developed to correct snowfall and not SWE directly. For example, Kang et al. [37] adjusted snowfall using a correction factor that was obtained by minimizing the difference between simulated SWE from a distributed hydrological model and observed SWE from PM. However, their approach is limited to shallow snowpacks because they assumed the reference SWE estimated from PM is unbiased, which does not hold true as the snowpack gets thicker. To our knowledge, no other study has presented an approach based on correcting SWE based on observed flow. Our WSC approach is relatively simple, requiring less complex modeling, and, therefore, is more user friendly.

Baseflow Separation as a Source of Uncertainty
As in any modelling approaches, the WSC has uncertainties in the modelling chain that affect results represented by the range of CF values obtained (see Figure 5 and Table 3). One uncertainty pertains to the HUT model's simplification of the snowpacks complex structure (e.g., stratigraphy, crystallography) when converting the PM signal into SWE in the GlobSnow product (see results section for more details).
Another potentially important source of uncertainty relates to the baseflow separation technique used for estimating direct runoff. Nathan and McMahon [16] recommend that the baseflow filter coefficient (β) of the recursive filter approach should assigned the value of 0.925. However, Ref. [15] stipulate that the Hydrun toolbox used in our study (which is based on the digital recursive method) was not intended for cold regions winter analysis because it was originally developed by associating a runoff event to a rainfall input event. Therefore, it is likely that the volume of direct runoff estimated using this method will be misrepresented for watersheds where slower melt events occur. This limitation likely contributed to errors, however, future versions of the toolbox is expected to include more detailed characterizations of snowmelt processes [15].
Using the current iteration of the toolbox, this source of uncertainty can be reduced by optimizing the β factor for each watershed. To do so, indirect approaches are required because baseflow is seldom directly measured. Gan and Zuo [17] proposed to use a hydrological model to fit the β factor to match simulated baseflow; however, the resulting optimized β factor will be dependent of the structure of the hydrological model used to perform baseflow matching. In contrast, we optimized the β factor by minimizing the distance between GlobSnow's SWE max and observed SWE max and, thus, our method is not model dependent. However, the drawback is that in situ SWE max measurements are required, which negates the WSC approach's prime advantage of not requiring in situ SWE measurements.
The optimization process is summarized in Table 5 that details the % bias of the corrected GlobSnow SWE max relative the reference SWE max for different baseflow β factors. The % bias for the recommended β value of 0.925 is shown in this table. With this β value the average bias was reduced to 18%, notably better than the Globsnow products 33.5% existing bias. While bold values are % bias which were minimized by adjusting β. Results indicate that using a single value for β fails to optimally correct for the underestimation of GlobSnow SWE max values. Instead, β values vary from 0.800 to 0.995 for the 8 watersheds under analysis. Interestingly, the average β is 0.914, which is close to the 0.925 recommended value of [16]. Adjusted β values reduced the bias in the WSC approach as illustrated in Figure 6. Moreover, RMSE was reduced, and the degree of linear correlation overall increased, as shown in Table 6.
Optimizing β values therefore contributed to reducing overall uncertainty. Table 5. Performance of the WSC approach to replicate the maximum measured SWE for the studied watersheds as a function of baseflow separation factor β.

Watershed
In We demonstrate that with watershed specific optimized β values, one can improve SWE max estimation from GlobSnow. However, this is at the expense of requiring in situ SWE max measurements. One way to avoid requiring in situ SWE, would be to develop a relationship between β and watershed physiographic and/or hydroclimatic characteristics and apply such relationship to watersheds where SWE is not measured.
This was attempted here, although a small number of watersheds were analysed which limits the robustness of our findings. Therefore, more watersheds need to be incorporated into such analysis to confirm the reported relationships below.
As a first attempt, we investigated the relationship between β and in situ SWE max , see Figure 8. We found a strong linear relationship, with a R 2 value of 0.798; deeper snowpacks, which have larger volumes of spring melt runoff, had larger β values. Recall that larger β value is associated with a lower contribution of baseflow to total runoff, while a lower β value means that the baseflow contribution is more significant to total runoff. This behaviour aligns with the nature of infiltration in northern watersheds, which is restricted and, sometimes almost completely limited, by a layer of frozen ground. Therefore, saturation overland flows occur as the soil above the frozen soil layer becomes saturated and any additional snowmelt (and rainfall) causes runoff.
Water 2022, 14, x FOR PEER REVIEW 18 of 24 As a first attempt, we investigated the relationship between β and in situ SWEmax, see Figure 8. We found a strong linear relationship, with a R 2 value of 0.798; deeper snowpacks, which have larger volumes of spring melt runoff, had larger β values. Recall that larger β value is associated with a lower contribution of baseflow to total runoff, while a lower β value means that the baseflow contribution is more significant to total runoff. This behaviour aligns with the nature of infiltration in northern watersheds, which is restricted and, sometimes almost completely limited, by a layer of frozen ground. Therefore, saturation overland flows occur as the soil above the frozen soil layer becomes saturated and any additional snowmelt (and rainfall) causes runoff. Our second attempt was to unravel possible relationships between  values and watershed physiographic characteristics (i.e., geology/pedology) as this influences baseflow. Indeed, baseflow-which is not to confused with groundwater flow-is the proportion of river runoff from stored sources such as soils, surficial deposits and permeable rocks. We observed that the watersheds located in the Appalachian region (0PL005, 02PG022, 01AD003 and 01010000) had lower  values, ranging from 0.800 to 0.900, whereas the watersheds located in the Canadian Shield (0RH035, 02PA007 and 02NE011) had higher value of 0.995. These two regions are characterized by different geological and pedological characteristics that likely effect baseflow. The Canadian Shield is generally associated with thin soils and the occurrence of granitic, gneiss, and volcanic bedrock where the capacity to store water is small. For example, at study sites characterized as Canadian Shield in Ontario, [38,39] found that water storage capacity to be only 6 mm and 8 mm, respectively. In contrast, the Appalachian region is mostly characterized by sedimentary bedrock covered with fine-to-deep glacial deposits where flow through fractured sediment rock and superficial deposits are more likely to occur. We argue that these different geologic and pedologic regimes of these regions may explain differences observed in the  values of watersheds, where lower/higher  values are associated with higher/lower baseflow contribution to total runoff. Additional analyses using more watersheds are needed to validate this hypothesis. Our second attempt was to unravel possible relationships between β values and watershed physiographic characteristics (i.e., geology/pedology) as this influences baseflow. Indeed, baseflow-which is not to confused with groundwater flow-is the proportion of river runoff from stored sources such as soils, surficial deposits and permeable rocks. We observed that the watersheds located in the Appalachian region (0PL005, 02PG022, 01AD003 and 01010000) had lower β values, ranging from 0.800 to 0.900, whereas the watersheds located in the Canadian Shield (0RH035, 02PA007 and 02NE011) had higher value of 0.995. These two regions are characterized by different geological and pedological characteristics that likely effect baseflow. The Canadian Shield is generally associated with thin soils and the occurrence of granitic, gneiss, and volcanic bedrock where the capacity to store water is small. For example, at study sites characterized as Canadian Shield in Ontario, Refs. [38,39] found that water storage capacity to be only 6 mm and 8 mm, respectively. In contrast, the Appalachian region is mostly characterized by sedimentary bedrock covered with fine-to-deep glacial deposits where flow through fractured sediment rock and superficial deposits are more likely to occur. We argue that these different geologic and pedologic regimes of these regions may explain differences observed in the β values of watersheds, where lower/higher β values are associated with higher/lower baseflow contribution to total runoff. Additional analyses using more watersheds are needed to validate this hypothesis.

Other Sources of Uncertainty
In addition to snowpack structure and baseflow separation techniques, other sources of uncertainties in the WSC approach (see Equation (5)) include precipitation during the melt period (P), modelling infiltration into frozen ground (I), and total runoff as measured at the gaging station (RO total ).
The accuracy of total precipitation (P) during the melt period is source dependent. Here we used the ERA5 daily precipitation that, according to Xiong et al. [40], shows good accuracy especially for the winter in our study region. Although, Crossett et al. [41] indicate that ERA5 underestimates precipitation along the Atlantic coast but overestimates precipitation inland compared to in situ observational weather stations. Despite that, ERA5 and similar products are particularly well suited for regions where weather stations are scarce and for larger watersheds, such as Manic5. However, for smaller watersheds and where weather station density is high, such as the study sites located in the Appalachian region, it may be preferable to use in situ observations. In this study, P was estimated from the ERA5 global reanalysis and therefore contains its uncertainties. This may affect the accuracy of the CF derived, impacting the corresponding corrected GlobSnow product.
Frozen ground infiltration was estimated using Gray's model (Equation (8)), which requires average near-surface (40 cm) soil saturation (S I ) and temperature (T I ) that were retrieved from the ERA5 reanalysis product. The ERA5 reanalysis product uses soil texture information to simulate S I and T I , where soil texture is categorized using the seven soil types derived from the root zone data of the FAO/UNESCO Digital Soil Map of the World [42]. However, this coarse representation of the true soil types inevitably introduced errors in S I and T I, and incorporating these products into the limited infiltration equation may have enhanced the bias associated with the algorithm and led to additional errors. Future works will benefit from the addition of in situ soil data that more accurately represents infiltration.
The opportunity time (t 0 ) is another variable affecting estimation of the infiltration into frozen ground, which is calculated according to Equation (9). For simplicity, we determined t 0 from the ERA5 SWE max and snow cover extent (SCE) product. Although the time at which SWE max occurs is rather straightforward to obtain, the time at which snow cover completely disappears is more difficult to retrieve because of the relatively low spatial resolution of the product. One way to improve future estimations of t 0 would be to use high spatial resolution snow cover data sets when available.
Given that t 0 corresponds to the duration of the snowmelt period, which was used to calculate the amount of precipitation (P), any error in estimating t 0 will also have repercussions on the estimation of P. Ideally, one should use the time during which surface runoff is occurring to calculate P across this period. This could be retrieved from a baseflow separation analysis; however, in our work, given the uncertainty at estimating the β filter parameter, we decided to approximate the runoff duration time as the duration of the snowmelt period. For small to medium size watersheds, this will have a minimal impact as the time of concentration is in the order of hours to a few days. However, in larger watersheds such as Manic5, the time of concentration will be several days and, therefore, will potentially impact the estimation of total precipitation during the opportunity time.
A rating curve is used for converting surface water elevation into river discharge. The curve is a relationship between water levels and discharges at a specific cross-section of the stream or the river and is established using gauged discharges and water levels. Rating curves can lead to significant errors. These errors can arise from a suite of operational factors (e.g., the number of sampling points and location of gauged section) and the uncertainty in discharge measurements might be as high as 20% of the observed value [43]. Furthermore, since discharge measurements are often impracticable during high floods, extrapolation errors are generally introduced. Additionally, rating curves may provide poor estimates of river discharge under ice conditions, as the river cross-section is constricted and produces higher river stage compared to ice-free periods. Any increase/decrease of river flow will directly impact CF values, see Equation (5). For example, assuming in a hypothetical watershed a ±10% uncertainty on river discharge, where P = 20 mm, I = 0, SWE max = 250 mm, total runoff = 350 mm and baseflow = 30 mm would result in CF ranging between 1.20 and 1.52.

The SWC Approach as Part of an Ensemble Flow Forecast System
Although the above uncertainties result in a range of watershed-scale CF values, we argue that this is an appealing feature when applying the WSC approach to probabilistic flow forecasts. An ensemble of corrected GlobSnow SWE values can be used to set initial conditions of a hydrological model for generating ensemble flow forecasts. As mentioned above, a 5-km GlobSnow SWE product is available 12 h after global satellite data has been acquired (data available in the Copernicus site). Therefore, the raw SWE product can be used to generate near real time corrected products. This is particularly useful to agencies concerned by spring flow forecasts, e.g., hydropower companies and river forecast centers, as forecasting flow scenarios becomes possible with estimates of corrected SWE at near real time, along with forecasted temperature and precipitation information.
While the WSC approach strictly applies to correct SWE max estimates from PM data (here GlobSnow), in the context of operational flood forecasting, one needs to correct SWE at any time a flow forecast is required. This requires having CF varying over time during snowpack build-up. Here, we proposed a simple linear relationship between CF and SWE from GlobSnow, see Figure 2, with a lower bound of CF = 1 and the upper value calculated according to equation 5. Whether or not a linear relationship between CF and SWE truly reflects the expected evolution of CF across the snowpack buildup remains to be validated by comparing resulting SWE values against in situ SWE(t) measurements.
In a probabilistic forecast world, the WSC approach produces an ensemble of SWE max values. Combining the ensemble with the CF(SWE) relationship (Figure 2), results in an ensemble of corrected SWE values each time a forecast is issued. Therefore, the corresponding probability distribution of corrected SWE at a given time only depends the corrected SWE max distribution and on CF-SWE relationship. Beneficially, no exogeneous information (e.g., observed SWE at the time a forecast is issued) is involved in defining the distribution. With better knowledge on the trajectory of SWE as the snowpack gradually accumulates (e.g., whether a thick or a thin snowpack is anticipated at the onset of the spring melt), one may update the probability distribution of the corrected SWE max . In other words, having the prior distribution of SWE max , one would obtain the posterior distribution of SWE max with knowledge of SWE at time t using a statistical inference method such as Bayesian inference: where P(H|E) is the posterior probability of H given E, P(H) is the prior probability of H and P(E|H) is the probability of observing E given H. In this equation, H = SWE max and E = SWE at a given time t. P(H), P(E), and P(E|H) can all be obtained from historical observations of SWE, including SWE max . Alternatively, the most likely SWE max value in a given year could be computed/updated each time new exogenous information affecting SWE max becomes available. This would require leveraging historical data to develop a statistical model between SWE max and exogenous variables at time t, such as the current SWE value. Using this approach, the most likely SWE max value computed at the time a forecast is issued would be used to correct the current SWE value using the CF-SWE relationship (Figure 2).

Validation of the WSC Approach: An Issue of Scale Mismatch
In this study, we describe a novel approach for correcting a remote sensing derived SWE product that avoids using in situ SWE measurements. However, as with the de-velopment of any quantitative remote sensing algorithm, we are limited in our ability to effectively validate the approach due to the severe spatial scale discrepancy between the data source used as the benchmark (i.e., in situ SWE measurements), and model product (i.e., corrected GlobSnow SWE). As pointed out by Wu [44] (p. 1769), 'Scale effects constrain the accuracy of retrieval and limit the development of remote sensing applications'. Therefore, it is important to match the reference SWE and the corrected SWE product to a common scale, here the watershed scale, for proper validation. The simplest upscaling approach is the arithmetic mean (a member of the Area-Weighted Methods family), which was used in this work.
In situ SWE sites are carefully selected to be representative surrounding areas by avoiding local effects on SWE values; for example, sites where wind effects are significant (resulting in either underestimating or overestimating the 'true' SWE) are avoided. Despite such cautions, accurate local SWE measurements may not be representative of the watershed scale SWE if the number of measurement sites is small, as is the case for this study.
Benchmark SWE values were obtained from only one or two in situ SWE measurements for all watersheds under study, except for the Manic 5 watershed, where 5 stations were used (locations of in situ measurements illustrated in Figure 2). Therefore, the low number of in situ stations used to obtain the reference SWE may explain situations where poor fits between the reference and the corrected SWE were obtained, in addition to the various sources of uncertainty discussed above. For example, watershed 01AD003 in 2010 was found to have a poor fit, where the corrected GlobSnow SWE max differed from the reference SWE max by 140 mm, resulting in a 60% deviation from the reference (note that this was obtained using β = 0.925). An improved data archive or SWE max would not only allow for better validation of the WSC approach but would also help in better defining relationships between the optimal β factor and SWE max discussed above. Datasets such as SNODAS [45] and CanSWE [46] may be considered for setting the benchmark.

Conclusions
In this study, we demonstrate the usefulness of integrating additional physical hydrological processes, specifically direct runoff, in the bias correction of GlobSnow, a global SWE database derived from space-borne passive microwave sensors. Unlike other physically based approaches that are more computer intensive and complex to implement, the WSC approach can correct SWE at the watershed scale while preserving the physics of the temperature brightness algorithm used in SWE conversion.
As with all modeling approaches, there are issues relating to uncertainties with the WSC approach; of those, the most critical to the output was the baseflow separation technique used in estimating direct runoff. While the recommended value of the baseflow filter coefficient (β = 0.925) for the recursive filter approach, overall reduced bias for all studied watersheds, we found that optimizing the β value for watersheds better corrected for the underestimation of GlobSnow SWE max values reducing overall uncertainty. We found that the deeper the snowpack the larger the corresponding β value, which is in line with the limited nature of infiltration in northern watersheds. Incorporating in situ soil data will reduce uncertainty associated with using the Gray's model to determine infiltrability. Additional uncertainty comes from the rating curve as it may provide poor estimates of river discharge under icy conditions as CF is sensitive to any increase or decrease in the river flow.
Despite these uncertainties, the WSC approach can be used for probabilistic flow forecasts. GlobSnow produces a near real time 5-km SWE product and an ensemble of corrected GlobSnow SWE produced using the WSC approach could be used to initialize conditions of a hydrological model. This application may be particularly useful to water management organizations concerned about spring flows.
To conclude, our WSC approach performed acceptably, increasing the accuracy of SWE max estimates when compared to in situ measurements for the studied watersheds. However, due to a scale discrepancy between the data used as benchmark and model product, and the representability of the local measurements to be used as a representative for the watershed, there were situations of poor fit. An improved benchmark would allow for better validating the WSC approach. Future studies will be dedicated to improving the benchmark and other assumptions made in evaluation of the approach.
Author Contributions: C.W. and R.L. contributed to the conception and design. C.W. performed Data analysis; C.W. prepared the original manuscript draft and did edits; R.L. reviewed and edited the manuscript and supervised the study. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded as part of the Canadian Natural Science and Engineering Research Council (NSERC) Industrial Research Chair on the Application of Hydrometeorological Data from Satellite Images to Improve Hydrological Forecasting. Financial support came from the following organizations: Hydro-Quebec, Brookfield and the City of Sherbrooke.
Data Availability Statement: Observed SWE data was obtained through private communications and is not available for sharing.

Conflicts of Interest:
The authors declare no conflict of interest.