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

Abstract

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 quantified 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. Infiltration 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., flow 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 filter coefficient β of 0.925 for baseflow separation. Results show the usefulness of integrating direct runoff for bias correction of existing GlobSnow product at the watershed scale. In addition, potential benefits are offered using the recursive filter approach for baseflow 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 flood-related damages.

1. 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 important ecologically in terms of water availability. This is so because groundwater replenishment depends on the amount of snowmelt and unusually low SWE_max will affect summer low flows [2].

SWE, including SWE_max, are traditionally estimated from manual point snow surveys. These surveys are limited in their spatial (i.e., individual observation points) and temporal coverage (i.e., typically only begins once the snowpack is considered well developed) and, therefore, do not always well represent watershed state [3]. For smaller and more homogeneous environments, the interpolation of spatially sparse point snow observations from either snow courses or weather stations may be adequate for estimating SWE and SWE_max at the watershed scale. However, for larger and more heterogeneous watersheds, alternative methods are needed in order to capture the high spatial (and temporal) variability of the snowpack, combined with terrestrial heterogeneity, and therefore, estimate SWE and SWE_max at the watershed scale.

Passive microwave (PM) sensors have been used for decades to improve the spatial coverage of SWE estimates. These sensors quantify the brightness temperature naturally emitted from Earth across the microwave wavelengths and can be used to estimate SWE because snow grains of a developed snowpack scatters microwave radiation as it is transmitted through the snowpack [4]. This method has three main limitations Firstly, the microwave brightness temperature decreases as the snowpack gets deeper [5]. This methodology is limited to dry snow conditions, as wet snow is opaque to microwaves and the quantified microwave radiation originates from a thin layer at the top of the snowpack [6]. Secondly, this methodology underestimates SWE under deep snow conditions (i.e., when SWE exceeds approximately 150 mm). The main reason is that higher brightness temperature values are recorded as snow absorbs radiation from the bottom layer and emits a signal stronger than that of dry snow covering soil, which is not properly accounted for [7]. The exact limit depends on snowpack stratigraphy and grain size. Thirdly, Space-borne microwave radiometers, like the Advanced Microwave Scanning Radiometer-EOS (AMSR-E) and the Special Sensor Microwave Imager/Sounder (SSMI/S) provide SWE estimates with full spatial coverage under dry snow conditions. However, these estimates suffer issues related to the coarse spatial resolution of these sensors (approximately 25 km²). Snowpack properties are highly variable at spatial scales smaller than the satellites’ footprint (<25 km²), driven by changes in vegetative and topographic controls [4,8], this fine-scale heterogeneity results in measurement errors in SWE estimates. Therefore, stand-alone algorithms for these PM observations can be suitable for mesoscale to continental scale applications but may not provide sufficient accuracy at the watershed scale for operational hydrology application [9], save for very large watersheds.

Methods that combine remote sensing with ground-based measurements, for example, the algorithm in described in Pulliainen [10], can help overcome limitations of either method alone and have further improved SWE estimates. Pullianinen’s algorithm combines passive microwave radiometer data with ground-based synoptic snow depth observations in the Helsinki University of Technology (HUT) snow emission model to estimate SWE [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 pixel resolution and used in an assimilation procedure for estimating SWE [10]. GlobSnow applies this approach using PM data from Nimbus-7 Scanning Multichannel Microwave Radiometer (SMMR), The Defense Meteorological Satellite Program (DMSP) SSM/I and DMSP SSMI/S to produce a 25-km SWE product over the Northern hemisphere, as described in Takala et al. [12] and Luojus et al. [11]. Further refinement to the original GlobSnow dataset uses land cover/land use data [13] to downscale SWE estimates to a 5-km product, which is available in near-real time.

With these advances, the GlobSnow SWE product is at a spatial resolution appropriate for hydrological applications at the watershed scale; however, limitations regarding the issue of SWE underestimation in deeper snowpack remains. Several methods have been 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.

2. Materials and Methods

2.1. 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:

DRH = DRG

(1)

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:

$${\mathrm{SWE}}_{\max }^{\mathrm{corr}}{=\mathrm{SWE}}_{\max }^{\mathrm{Glob}}\ast \mathrm{CF}$$

(2)

where, ${\mathrm{SWE}}_{\max }^{\mathrm{Glob}}$ is the SWE_max estimated from GlobSnow, ${\mathrm{SWE}}_{\max }^{\mathrm{corr}}$ 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:

$${\mathrm{DRH}=\mathrm{RO}}_{\mathrm{total}}-\mathrm{B}$$

(3)

$$\mathrm{DRG}=\left({\mathrm{SWE}}_{\max }^{\mathrm{corr}}-\mathrm{I}+\mathrm{P}\right)\times \mathrm{A}$$

(4)

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.

Combining Equations (1)–(4), one obtains the equation for CF:

$$\mathrm{CF}=\frac{\frac{{\mathrm{RO}}_{\mathrm{total}}-\mathrm{B}}{\mathrm{A}}-\mathrm{P}+\mathrm{I}}{{\mathrm{SWE}}_{\max }^{\mathrm{Glob}}}$$

(5)

Given that ${\mathrm{SWE}}_{\max }^{\mathrm{Glob}}$ 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:

• ${\mathrm{SWE}}_{\max }^{\mathrm{Glob}}$ 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.
• 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.
• Flow measurements also have errors related to the rating curve used to convert stream water level into streamflow.
• 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).

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]:

$${\mathrm{B}}_{\mathrm{t}}=\mathsf{\beta }·{\mathrm{B}}_{\mathrm{t}-1}+\frac{1-\mathsf{\beta }}{2}·{(\mathrm{Q}}_{\mathrm{t}}{+\mathrm{Q}}_{\mathrm{t}-1}),$$

(6)

$${\mathrm{q}}_{\mathrm{t}}{=\mathrm{Q}}_{\mathrm{t}}-{\mathrm{B}}_{\mathrm{t}},$$

(7)

where ${\mathrm{B}}_{\mathrm{t}}$ is the baseflow at time t, ${\mathrm{q}}_{\mathrm{t}}$ is the direct runoff, ${\mathrm{Q}}_{\mathrm{t}}$ represents the original streamflow and $\mathsf{\beta }$, is a filter coefficient which affects the magnitude of the baseflow values. For example, increasing $\mathsf{\beta }$ produce a lower baseflow. Researchers used prescribed values for $\mathsf{\beta }$ to obtain optimal results [17,18]. Gan and Zuo [17] set the $\mathsf{\beta }$ 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 $\mathsf{\beta }$ value should be 0.925. It is advised to thoroughly investigate the $\mathsf{\beta }$ 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 $\mathsf{\beta }$ 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:

$${\mathrm{INF}=\mathrm{C}\mathrm{S}}_0^{2.92} {\left(1-{\mathrm{S}}_{\mathrm{I}}\right)}^{1.64}(\raisebox{1ex}{$273.15-{\mathrm{T}}_{\mathrm{I}}$}\!\left/ \!\raisebox{-1ex}{$273.15$}\right.{)}^{-0.45}{\mathrm{t}}_0^{0.44},$$

(8)

where INF if the cumulative infiltration (mm), S₀ is the soil moisture content at the soil surface, C is the land cover type coefficient, t₀ is the infiltration opportunity time (h), S_I and T_I are the average soil saturation (mm³/mm³) and temperature (K), respectively, of the first 40-cm soil layer at the start of infiltration. S₀ 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₀ 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₀:

$${\mathrm{t}}_0{=\mathrm{Time}}_{\mathrm{SWE}\left(\mathrm{MAX}\right)}-{\mathrm{Time}}_{\mathrm{SWE}\left(\mathrm{ZERO}\right)},$$

(9)

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 ${\mathrm{SWE}}_{\max }^{\mathrm{Glob}}$ estimates based on the direct runoff that occurs during the spring melt season. The same correction factor cannot be applied to correct ${\mathrm{SWE}}^{\mathrm{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.

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.

2.2. 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² (02PG022—Ouelle River Watershed) to 24,698 km² (Manic 5 watershed), see

Table 1. Physiographic and hydrometeorological data of the watersheds under study.

				Long Term Averaged Measured Values
Watershed ID	Station Name	Data Period	Area (km2)	SWEmax (mm)	Peak Flow (m3/s)
02RH035	aux Écorces River at highway-bridge 169	1985–2013	1110	298	166
01AD003	St. Francis River at outlet of Glasier lake	1985–2016	1359	231	200
Manic5	Manicouagan 5	2006–2016	24,698	269	2018
01010000	St. John River at Ninemile Bridge	1985–2016	3473	225	690
02NE011	Croche River downstream of Changy Brook	1985–2013	1570	274	190
02PA007	Batiscan River downstream of des Envies River	1985–2013	4480	280	551
02PL005	Bécancour River upstream Palmer River	1985–2013	919	217	214
02PG022	Ouelle River near St. Gabriel de Kamouraska	1986–2013	795	252	207

. 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³/s [28] and annual peak flows exceeding on average 2000 m³/s, occurring in May. SWE can reach values of 300 mm or more [29]. Flow from Manic 5 drains into a 2000 km² hydropower reservoir, with an installed capacity of 2.6 GW [30].

Hydrometeorological information about these watersheds can be found in Table 1.

2.3. Data

2.3.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. GlobSnow and in situ seasonal maximum SWE_max data.

Watershed ID	Data Period	Average SWEmax
		GlobSnow Product	In Situ Measurements
02RH035	1985–2013	171	298
01AD003	1985–2016	178	231
Manic5	2006–2016	185	269
01010000	1985–2016	162	225
02NE011	1997–2013	158	274
02PA007	1985–2013	162	279
02PL005	1985–2013	151	216
02PG022	1985–2013	181	252

, showing that GlobSnow underestimates SWE compared to in situ measurements.

2.3.2. ERA5

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].

2.3.3. 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].

2.3.4. 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.

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. Computed CF values for all watersheds under study.

Watershed ID	Area	GlobSnow Average SWEmax (mm)	CF
			Average	Std Dev	Min	Max
02RH035	1110	171	1.37	0.43	1.01	2.57
01AD003	1359	178	1.43	0.37	1.01	2.52
Manic5	24,698	185	1.43	0.41	1.01	2.51
01010000	3473	162	1.78	0.53	1.10	3.47
02NE011	1570	158	1.29	0.29	1.02	2.27
02PA007	4480	162	1.36	0.33	1.01	2.66
02PL005	919	151	1.91	0.52	1.20	4.08
02PG022	795	181	1.74	0.40	1.15	2.77

.

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 rain-on-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. Statistical summary of the original (Glob) and corrected (C-Glob) SWE_max. A negative bias means that GlobSnow and corrected GlobSnow underestimate SWE relative to the reference (in situ) SWE.

Watershed	Average SWEmax (mm)		Range (mm)		Bias (%)		Correlation Coefficient		RMSE (mm)
	Glob	C-Glob	Glob	C-Glob	Glob	C-Glob	Glob	C-Glob	Glob	C-Glob
02RH035	171	228	119–276	137–357	−42.4	−23.5	0.525	0.636	136.8	85
01AD003	178	241	108–273	137–410	−23.2	4.54	0.646	0.671	66.4	50.5
Manic5	185	263	112–240	177–321	−31.4	−2.2	0.227	0.691	103.3	38.6
01010000	162	264	104–277	114–417	−28.2	17.1	0.461	0.475	86.3	78.1
02NE011	158	196	93–277	93–314	−42.2	−28.3	0.756	0.639	123.5	93.1
02PA007	162	208	111–265	121–319	−42.1	−25.8	0.791	0.691	124.9	85.6
02PL005	151	275	103–232	142–558	−30.2	27.1	0.711	0.442	74.4	89.1
02PG022	181	292	127–281	132–441	−28.1	15.9	0.583	0.306	83.9	85.8

. 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.

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.

4. 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.

4.1. 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. Performance of the WSC approach to replicate the maximum measured SWE for the studied watersheds as a function of baseflow separation factor β.

Watershed	In Situ SWE Max	% Difference from In Situ SWEmax Measurements
		GlobSnow	Baseflow Separation Factor β
			0.800	0.825	0.850	0.900	0.925	0.975	0.995
02RH035	298.1	−42.4	−42.4	−38.0	−35.6	−29.2	−23.5	−8.1	−1.9
01AD003	230.9	−23.2	−19.1	−17.1	−14.3	−3.7	4.5	27.1	40
Manic5	269.2	−31.4	−24.3	−22.1	−16.8	−6.5	−2.2	13.4	21.6
02PA007	279.7	−42.1	−41.4	−40.4	−38.7	−31.7	−25.8	−8.4	−0.7
02PG022	251.8	−28.1	−9.2	−6.1	−2.0	8.8	15.9	34.0	42.4
02PL005	216.6	−30.2	−0.1	3.9	8.42	19.8	27.1	46.0	60.7
02NE011	274.0	−42.2	−42.2	−40.4	−39.3	−33.7	−28.3	−14.1	−5.4
1010000	225.4	−28.2	−10.1	−6.3	−1.8	9.5	17.1	35.8	46.1

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. Performance metrics of the WSC approach using the recommended β and the optimized β for the watersheds.

Watershed	Average SWEmax (mm)		Range (mm)		Bias (%)		Correlation Coefficient		RMSE (mm)
	β = 0.925	βopt	β = 0.925	βopt	β = 0.925	βopt	β = 0.925	βopt	β = 0.925	βopt
02RH035	228.1	292.5	137–357	173–484	−23.5	−1.9	0.636	0.653	85	55.2
01AD003	241.4	222.5	137–410	137–384	4.5	−3.7	0.671	0.697	50.5	44.2
Manic5	263.2	263.2	177–321	177–321	−2.2	−2.2	0.691	0.691	38.6	38.6
02PA007	207.5	281.5	121–319	121–443	−25.8	0.7	0.691	0.570	85.6	62.1
02PG022	291.9	246	132–441	130–365	15.9	−2.0	0.306	0.286	85.8	65.7
02PL005	275.2	216.5	143–558	139–443	27.1	0.1	0.442	0.360	89.1	58.0
02NE011	196.6	259	93–315	93–406	−28.3	−5.3	0.639	0.560	93.1	68.8
01010000	263	221	114–417	114–345	17.1	−1.8	0.475	0.485	78.1	59.3

. Optimizing β values therefore contributed to reducing overall uncertainty.

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² 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, 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.

4.2. 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₀) is another variable affecting estimation of the infiltration into frozen ground, which is calculated according to Equation (9). For simplicity, we determined t₀ 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₀ would be to use high spatial resolution snow cover data sets when available.

Given that t₀ corresponds to the duration of the snowmelt period, which was used to calculate the amount of precipitation (P), any error in estimating t₀ 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.

4.3. 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:

$$\mathrm{P}\left(\mathrm{H}|\mathrm{E}\right)=\frac{\mathrm{P}\left(\mathrm{E}|\mathrm{H}\right)\times \mathrm{P}\left(\mathrm{H}\right)}{\mathrm{P}\left(\mathrm{E}\right)}$$

(10)

where $\mathrm{P}\left(\mathrm{H}|\mathrm{E}\right)$ is the posterior probability of H given E, $\mathrm{P}\left(\mathrm{H}\right)$ is the prior probability of H and $\mathrm{P}\left(\mathrm{E}|\mathrm{H}\right)$ is the probability of observing E given H. In this equation, H = SWE_max and E = SWE at a given time t. $\mathrm{P}\left(\mathrm{H}\right)$, $\mathrm{P}\left(\mathrm{E}\right)$, and $\mathrm{P}\left(\mathrm{E}|\mathrm{H}\right)$ 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).

4.4. 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 development 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.

5. 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.