Interrelationships between MODIS/Terra remotely sensed snow cover and the hydrometeorology of the Quesnel River Basin, British Columbia, Canada

A spatial filter (SF) method is adopted to reduce the cloud coverage from the Moderate Resolution Imaging Spectroradiometer (MODIS) 8-day snow products (MOD10A2) between 2000–2007 in the Quesnel River Basin (QRB) of British Columbia, Canada. A threshold of k = 2 cm of snow depth measurements at four in-situ observation stations in the QRB are used to evaluate the accuracy of MODIS snow products MOD10A1, MOD10A2, and SF. Using the MOD10A2 and the SF, the relationships between snow ablation, snow cover extent (SCE), snow cover fraction (SCF), streamflow and climate variability are assessed. Based on our results we are able to draw several interesting conclusions. Firstly, the SF method reduces the average cloud coverage in the QRB from 15% for MOD10A2 to 9%. Secondly, the SF increases the overall accuracy (OA) based on the thresholdk = 2 cm by about 2% compared to MOD10A2 and by about 10% compared to MOD10A1 at higher elevations. The OA for the four in-situ stations decreases with elevation with 93.1%, 87.9%, 84.0%, and 76.5% at 777 m, 1265 m, 1460 m, and 1670 m, respectively. Thirdly, an aggregated 1C rise in average air temperature during spring leads to a 10-day advance in reaching 50% SCF (SCF 50%) in the QRB. The correlation coefficient between normalized SCE of the SF and normalized streamflow is −0.84 (p<0.001) for snow ablation seasons. There is a 32-day time lag for snow ablation to impact the streamflow the strongest at the basin outlet. The linear correlation coefficient between SCF 50% and 50% normalized accumulated runoff (R 50%) attains 0.82 Correspondence to: J. Tong (jtong@unbc.ca) (p<0.01). This clearly demonstrates the strong links that exist between the SCF depletion and the hydrology of this sub-boreal, mountainous watershed.


Introduction
The snow cover extent (SCE) and snow water equivalent (SWE -defined as the liquid water equivalent depth of snow) are important parameters for many land surface and hydrological models (Hall and Martinec, 1985;Jain and Lall, 2000;Déry et al., 2004Déry et al., , 2005;;Yang et al., 2003Yang et al., , 2007;;Liston, 1999;Kolberg and Gottschalk, 2006;Kolberg et al., 2006).This arises from the significant contribution of snow ablation to river discharge in high-latitude or alpine watersheds (Barnett et al., 2005;Schmugge et al., 2002;Dyer, 2008).However, during the latter half of the 20th century, snow cover has decreased significantly during spring over North America (Groisman et al., 2004;Déry and Brown, 2007).Due to increasing air temperatures, snow has begun to ablate about eight days earlier in northern Alaska compared to the mid-1960s (Stone et al., 2002.In addition, declines and earlier occurrences of maximum SWE have been observed in the Cordillera of western North America (Mote et al., 2005;Mote, 2006;Stewart et al., 2005;Déry et al., 2009).
The spatial and temporal variability of SCE and SWE are critical variables for many hydrometeorological studies.Yang et al. (2003) and Zhou et al. (2005) show that there exists a significant relationship between National Oceanic and Atmospheric Administration (NOAA) weekly SCE and Moderate Resolution Imaging Spectroradiometer (MODIS) remotely-sensed SCE changes and measured streamflow in the large Siberian watersheds and Upper Rio Grande River Basin, respectively.However, there are limited studies on the relationships between SCE and the hydrometeorology of the alpine sub-boreal forest owing to sparse in-situ observations and the limited accuracy of remote sensing products in this environment.Therefore, the Quesnel River Basin (QRB), a representative high-elevation watershed in the sub-boreal forest of British Columbia (BC), Canada, is selected as a test site to analyze the relationships between its SCE and hydrometeorology.
Tremendous technological achievements in the last decades make it possible for scientists to study the Earthatmosphere system at virtually every location on Earth with the rapidly increasing number of satellite platforms of remote sensors (König et al., 2001).Optical remote sensors such as MODIS on Terra can detect snow cover owing to the higher spectral reflectance of snow in the visible bands compared to the mid-infrared bands.The overall accuracy (OA -defined as the number of pixels correctly classified divided by the total number compared) of MODIS snow products is typically around 90% in some areas such as the Tibetan Plateau (Pu et al., 2007), Austria (Parajka and Blöschl, 2006), the Upper Rio Grande River Basin of Colorado (Klein and Barnett, 2003;Zhou et al., 2005), the Missouri and Columbia River Basins, USA (Maurer et al., 2003) and the Kuparuk River Basin of Alaska (Déry et al., 2005) in cloud-free conditions.However, cloud coverage significantly affects the OA of MODIS snow products owing to the inability of radiation in the visible and mid-infrared electromagnetic spectrum to penetrate clouds.To diminish the impact of clouds on the detection of snow by remote sensing, Parajka and Blöschl (2008) developed a spatio-temporal method applied to MODIS/Terra and MODIS/Aqua snow products.This approach successfully reduces the cloud coverage from an average of 63% in the daily snow products to about 4% with OA above 90% in Austria.However, since MODIS/Aqua was launched in July 2002, the spatio-temporal methods cannot be used to improve the snow products prior to this time.Furthermore, the OA of MODIS snow products in sub-boreal forests and mountainous areas such as the QRB remains unknown.
In this paper, an alternative simplified method entitled the spatial-filter method (SF) is proposed to reduce cloud coverage based on the MODIS/Terra 8-day snow products (MOD10A2) in the QRB.Then the SF is evaluated by comparing the different cloud coverage between MODIS/Terra daily snow products (MOD10A1), MOD10A2 and SF.In addition, the OA of these different snow products is evaluated according to the snow depth (SD) measurements at four insitu stations in the QRB.Finally, the relationships between snow ablation, SCE, streamflow and meteorology are explored based on the MOD10A2 and the SF datasets.

Study area
The QRB is located in the northern part of the Fraser River Basin (FRB), BC, Canada, and is centered near 52.5 • N and 121 • W. The QRB covers about 12 023 km 2 over a wide range of topography, meteorological conditions and land uses across the watershed.Elevation ranges from ∼ 500 m to 3000 m with an average value of 1375 m.The higher portion lies in the northeast, covered mainly by glaciers, alpine terrain and old growth forests.The southwest portion is the lower elevation area covered mainly by young forests, agricultural regions, and small settlements.The precipitation across the QRB ranges from 500-2500 mm yr −1 .Burford et al. (2009) provide full details on the land use, topography, ecological features, and hydroclimatological aspects of the QRB.The land cover distribution (forests, agriculture, surface water, alpine tundra and glaciers) in the QRB is highly representative of the other sub-basins of the FRB landcover types.In fact, the hypsometry of the QRB nearly matches that of central BC.Four in-situ observation stations operated by Environment Canada and the BC Ministry of Environment are located within the QRB.Streamflow for the Quesnel River (QR) is measured at a gauge situated near its outlet to the Fraser River (Fig. 1).

MODIS snow products
MODIS, with 36 discrete, narrow spectral bands from approximately 0.4 to 14.4 µm, is mounted on the Terra spacecraft, launched on 18 December 1999.Terra orbits across the equator at a local standard time of ∼10:30 a.m. on its descending node.The spatial resolution of MODIS bands range from 250 m to 1000 m with a spectral resolution of 0.01 µm to 0.05 µm for different bands.Bands 1 and 2 have 250 m spatial resolution; bands 3 to 7 have 500 m spatial resolution; and bands 8 to 36 have 1000 m spatial resolution.An identical MODIS instrument is mounted on the Aqua spacecraft, launched on 4 May 2002 to provide complimentary snow products.MODIS/Terra began to collect Earth system information on 24 February 2000 (http://modis.gsfc.nasa.gov/).
The SNOWMAP algorithm for MODIS snow products is described fully by other authors (e.g., Klein et al., 1998;Hall et al., 1995Hall et al., , 2002) ) so only a brief summary of its details is given here.Briefly, the SNOWMAP algorithm is based on the normalized difference snow index (NDSI), which is the ratio between at-satellite reflectance in MODIS band 4 (0.55 µm) and band 6 (1.6 µm).Normally, a threshold value of 0.4 in the NDSI combined with the normalized difference vegetation index (NDVI) and the measured reflectance of the surface is employed to classify snow pixels.The snow maps are available at different resolutions and projections such as 500 m daily and 8-day maximum data on a sinusoidal  Snow cover extent (SCE) and snow cover fraction (SCF) in this paper denote respectively the total area of the pixels classified as snow and the ratio of number of the pixels classified as snow to the total number of pixels in the entire QRB or some of its smaller subbasins.The MOD10A1 and MOD10A2 at a 500 m resolution are highly suitable for analyzing the SCE and SCF for a watershed of an area of about 10 000 km 2 .The MOD10A2 is an 8-day composite snow product from every 8-day period of the MOD10A1 products starting on the first day of every year.The MOD10A2 effectively forms a temporal filter of MOD10A1 data that minimizes the cloud coverage while maximizing the SCF.For MOD10A2, a pixel is classified as cloud only when the pixel is cloud-covered continuously all 8 days.However, if a pixel is observed as snow cover for any of the days, the pixel is classified as snow in MOD10A2.If a pixel is classified as no snow throughout all 8 days, it is then classified as no snow in MOD10A2.
Version 005 (V5) is the latest release of MODIS snow products that began to be reprocessed in mid-July 2006.The reprocessing includes all the MODIS snow products beginning from 24 February 2000 to the present, and was completed in September 2008.The V5 MOD10A1 and MOD10A2 data used here begin on 24 February 2000 and end on 31 December 2007.Tile h10v03 covers the entire study area of the QRB.The MODIS Reprojection Tool (MRT) is used to resize and reproject the tile data to BC Albers equal area projection with batch processing code.

In situ meteorological data
There are nine in-situ observation stations including automated snow pillows and weather stations in and around the QRB that are operated by Environment Canada and the British Columbia Ministry of Environment, Forest and Range (Fig. 1).Table 1 lists the details of the nine stations.All nine stations measure air temperature and precipitation, which are available online (http://www.climate.weatheroffice.ec.gc.ca/climateData/canada e.html or http:// a100.gov.bc.ca/pub/aspr/).All stations measure daily SD with the exception of Barkerville (1A03P) and Likely.Snow occurs only 5 to 10% of the time equating to no more than 15-20 8-day periods with snow over eight years at Quesnel A, Kersley A, and Williams Lake A. Therefore, a single error in classifying snow as no snow will lead to more than a 5% underestimation error (UE), which lacks statistical meaning.In other words, the sample sizes are insufficiently large to conduct a robust statistical analysis of the accuracy of the MODIS snow products at the removed stations.Hence, these stations with sporadic snow cover are omitted from the accuracy assessments.Therefore, only the SD data measured at Horsefly Lake/Gruhs Lake (HLGL), Barkerville (1090660), Boss Mountain Mine (BMM) and Yanks Peak East (YPE) are used to verify the accuracy of MODIS snow products (see Fig. 1).All air temperature and precipitation measurements at the nine stations are used to explore the relationships between SCE or SCF and hydrometeorological variables in the QRB.

Streamflow
The QR gauge station is one of at least 6 active gauges throughout the QRB; however, the other gauges only represent smaller drainage areas in the QRB rather than the integrated response of the entire watershed (Burford et al., 2009).
In addition, there are no dams or other major anthropogenic disturbances on the QR.Thus daily streamflow measurements at this gauge (I.D. 08KH006) from the Water Survey of Canada (http://www.wsc.ec.gc.ca/hydat/H2O/) over the period 2000 to 2007 are used in the analyses.However to explore the impacts of Quesnel Lake, a large feature of the QRB (Fig. 1), on the overall streamflow, the Horsefly River subbasin (∼2760 km 2 ) upstream of Quesnel Lake is also incorporated in some of the analyses.Missing streamflow (2000)(2001)(2002)(2003)(2004)

Spatial filter for MOD10A2
The SF is used only to reclassify the 8-day MOD10A2 data since they already form a temporal filter for the original MOD10A1 product.This process of reclassifying the MOD10A2 using the SF is more straightforward than the spatio-temporal method used by others (Parajka and Blöschl, 2008).The steps in the SF are outlined in Fig. 2. First, the snow maps are reclassified as snow, no snow, or cloud based on the 12 classified types of MOD10A2.Except for the snow and no snow pixels, all other pixels such as missing, erroneous, and sensor-saturated data are defined as clouds.Then, a cloud-covered pixel is replaced by the majority of noncloud pixels in the eight closest neighbourhood pixels.If the number of pixels of snow and no snow are equal, the center pixel is defined as snow.However, if all of the eight closest neighborhood pixels are cloud-covered, the center pixel is kept as cloud.

No
Cloud for all other 8 points

Evaluation of MOD10A1, MOD10A2 and SF products
To evaluate the OA of MOD10A2 and SF products, ground observations of snow distributions are needed.However, there are few in-situ snow depth observations in the subboreal forest owing to the remote and mountainous terrain.In the QRB, four in-situ SD observation stations (HLGL, Barkerville, BMM, and YPE) are considered as ground truth for the corresponding pixels in the remotely sensed snow products.
There are four possible outcomes comparing in-situ observations with MODIS snow products (Table 2).The OA, UE, and overestimation error (OE) of MODIS snow products are evaluated quantitatively based on comparisons with in-situ observations.When both in-situ observations and MODIS snow products simultaneously infer the same surface type, such as snow or no snow, MODIS snow products are deemed as correct.In this work, the OA, UE, and OE are as defined by Pu et al. (2007) and Parajka and Blöschl (2008).The OA of MODIS snow products is calculated by Eq. ( 1) through the whole period: where terms a, b, c and d are defined in Table 2.If the in-situ observations show snow (or no snow) when the corresponding MODIS data resolve no snow (or snow), this is called UE (or OE).The UE and OE of MODIS snow products are calculated by Eq. ( 2) and (3), respectively: All accuracy and error estimates are based on days when the in-situ observations are available and the corresponding MODIS snow products are not cloud.
Owing to the heterogeneity of snow cover, a threshold of SD (k) has to be established to classify the ground truth as snow or no snow.When the SD≥k cm, the land surface is classified as snow; otherwise, the land surface is defined as bare.The OA of the snow products based on the different thresholds is compared to determine which threshold is better for classifying the ground truth.Three different thresholds of k = 0 cm, 2 cm, and 5 cm are initially used to calculate the OA of the snow products.This shows that the highest OA is achieved with k = 2 cm (not shown).Therefore, a commonly used threshold of 2 cm (Stieglitz et al., 2003;Neumann et al., 2006) is adopted to classify the land as snow or no snow at in-situ stations.For this evaluation, the daily snow conditions at the ground truth sites are reclassified to match the periods of the 8-day MODIS snow products.During an 8-day period, the ground is defined as snow only when there are at least 4 days with snow on the ground (Zhou et al., 2005).
Since there are uncertainties in the locations of in-situ stations (resolution of about 0.01 o ) and the MODIS reprojection tool, three different patches (1×1, 3×3, and 5×5 pixel 2 ) are initially adopted to retrieve the remotely sensed snow products of the pixels corresponding to the in-situ stations.For the 3×3 pixel 2 or 5×5 pixel 2 , the central pixels are defined as the majority of the 9 or 25 closest neighbourhood pixels.Zhou et al. (2005) find that the retrieval method of 5×5 pixel 2 had the highest OA compared to in-situ ground observations by the Snowpack Telemetry measurements in the Upper Rio Grande River Basin, New Mexico, USA.The OA of MODIS snow products from the three different retrieval methods show that 5×5 pixel 2 has the highest OA in the QRB (not shown).Therefore, only the retrieval method for 5×5 pixel 2 is used to evaluate the MODIS snow products in the QRB.

Eight day in-situ data
To explore the interrelationship between SCE and river discharge, the daily streamflow data measured at the gauges are averaged into 8-day values beginning from 1 January every year.To compare the MOD10A2 SCE to the corresponding streamflow, the 8-day averaged and accumulated streamflow are calculated over the same periods as MOD10A2.To study the effects of climate factors such as precipitation and temperature on the SCE and streamflow, hydrological year precipitation and temperature are calculated based on the insitu observations.Since the snow accumulation begins in September nearly every year in the QRB, the "hydrological year" in this paper is defined from 1 September to 31 August of the following year (similarly to Wang et al., 2009) The date at which the SCF attains 50% during each snow ablation season is referred to as SCF 50% , whereas the dates at which the normalized accumulated runoff (ratio between accumulated runoff at a given time and total snow ablation season accumulated runoff) reaches 50% of its seasonal value during the snow ablation season will be denoted by R 50% .This paper focuses on the interrelationships between SCE, SCF and hydrometeorological conditions during snow ablation seasons, defined as the period from 1 March to 30 June with a total of 128 SCE values from 2000 to 2007.To quantify the relationships between SCF 50% and average air temperature in the QRB in a given period, we define the day of year of SCF 50% and R 50% beginning from the first day of every year.In addition, 8-day average air temperatures and 8-day accumulated precipitation through the hydrological years are calculated based on the nine in-situ stations to compare with the SCF 50% and R 50% .

Cloud coverage of different MODIS snow products
The SF reduces cloud coverage over the QRB (Fig. 3).coverage than MOD10A1 owing to the temporally filtered algorithm.The maximum cloud coverage difference between MOD10A1 and MOD10A2 is 66.2%, occurring on 2 February, while the minimum cloud coverage difference between MOD10A1 and MOD10A2 is 29.1%, occurring on 13 August.Although the cloud coverage of MOD10A2 is much lower than those of MOD10A1, there remains ≈15% cloud coverage throughout the whole year.Therefore, cloud coverage results in a spurious decrease of remotely sensed SCF since, for example, it decreases about 18.7% from 76.8% on 1 January to 58.1% on 17 January.However, the SF reduces some of the cloud coverage in MOD10A2.The maximum amount of cloud coverage reduced by SF is 19.7% on 13 March 2004 from 42.5% of MOD10A2 to 22.9% of SF.The average cloud coverage decreases from an annual average of about 15% of MOD10A2 to 9% using the SF.The SCF of the SF has an average about 5% more than those of MOD10A2 in winter owing to the elimination of cloud coverage (Fig. 4).

Accuracy of different MODIS snow products
The evaluation is based on in-situ observations at the four stations in the QRB.The threshold of k = 2 cm for the insitu SD is set to classify the ground as snow or no snow.The evaluation is based on cloud free days when the MODIS snow products corresponding to the in-situ stations that show no cloud is present and the days when the in-situ observations are available.
Table 3 lists the cloudy days for 5×5 pixel 2 of MOD10A1, MOD10A2, and SF snow products and total days when the in-situ observations are available.During the total 2864 days from 24 February 2000 to 31 December 2007, there are 180, 366, 492, and 418 missing daily SD for HLGL, Barkerville, BMM, and YPE, respectively.A total of 361 8-day in-situ SD during the period can be calculated according to the daily SD.The cloudy days represent the days when the remotely sensed snow maps show clouds for the pixels cor-responding to the in-situ observations during the days when the in-situ SD are available.For MOD10A1, the percent of cloudy days ranges from 64.2% to 71.9%.The fraction of cloudy days calculated from MOD10A2, which ranges from 10.5% to 17.7%, is considerably lower than that inferred from MOD10A1.The SF reduces the percentage of cloudy days from MOD10A2 to 12.1%, 8.0%, 4.7%, and 3.6% for YPE, BMM, Barkerville, and HLGL, respectively.
The OA, UE, and OE for 5×5 pixel 2 of MOD10A1, MOD10A2, and SF based on k = 2 cm are listed in Table 4.The OA for the four in-situ stations ranges from 67.4% to 90.1%, 75.9% to 90.1%, and 76.5% to 93.1% for MOD10A1, MOD10A2 and SF snow products, respectively.The OA for the four in-situ stations varies owing to the different topography.HLGL, located at 777 m, has the highest OA.The lowest OA occurs at the highest station YPE, at an elevation about 1670 m.BMM and Barkerville, which are located at 1460 m and 1265 m, respectively, have intermediate OAs.For the same in-situ station, the SF has about 3.0%, 1.2%, 1.6% and 0.6% higher OA than MOD10A2 and 3.0%, 1.9%, 10.5%, and 9.1% higher OA than MOD10A1 for HLGL, Barkerville, BMM, and YPE, respectively.
The error rate is mainly due to the UE, which implies that the snow on the ground is incorrectly classified as nosnow in the remotely sensed snow maps.The UE at YPE reaches about 32.1%, 23.4%, and 23.1% for MOD10A1, MOD10A2 and SF, respectively; however, the OE at YPE is relatively low at about 0.6%, 0.7%, and 0.5% for MOD10A1, MOD10A2 and SF.BMM exhibits a similar pattern as YPE with much higher UE than OE.However, the UE and OE at HLGL are nearly equal.This phenomenon suggests that the complex topography in the sub-boreal forest is a main factor leading to the UE of snow from MODIS (Tong et al., 2009).This may be due to the more variable atmospheric conditions and topographical shading that prevent the remote sensing instrument in obtaining clear signals from the earth's surface.A more detailed comparison is shown in Table 5 to assess the bias of disagreement for MOD10A2 and SF.For all in-situ stations, the accuracy of classifying land as land is higher than that of classifying snow as snow.The SF increases the accuracy of classifying snow as snow by about 1.8%, 5.8%, 4.8%, and 7.9% for YPE, BMM, Barkerville, and HLGL, respectively.The error of classifying land as snow (2.6-7.2% for MOD10A2, 2.7-4.0%for SF) is much lower than the error of classifying snow as land (11.4-24.6% for MOD10A2, 12.9-27.4% for SF).The error of classifying snow as cloud is higher than the error of classifying land as cloud for all in-situ stations, with differences from 2.5% to 10.4%.

Interrelationships between SCE and streamflow
In sub-boreal forests and mountainous areas such as the QRB, the SCE and SCF have a strong relationship with streamflow (Fig. 4b).
The 8-day annual averages of SCF and streamflow from 2000 to 2007 show significant opposite trends with correlation coefficients of about -0.746 (p<0.001) and −0.750 (p<0.001) for MOD10A2 and SF, respectively.However, during the summer (especially July and August), snow is not the main source of streamflow although there is a perennial snow cover and/or glaciers at higher elevations (>2500 m; Tong et al., 2009).The 8-day time series of SCF derived from MOD10A2 and SF for the QRB and normalized accumulated runoff of the QR during snow ablation seasons from 2000-2007 are shown in Fig. 5.
The SF shows higher SCF than MOD10A2 during the snow ablation seasons, which can reach peak values of 97%.During spring, a decrease in SCF induces higher runoff in the QRB.The correlation coefficients between SCE and lagged streamflow attain their most significant values at about −0.84 (p<0.001) and −0.79 (p<0.001) for a 32-day lag for the SF and MOD10A2, respectively (Fig. 6a).The correlation coefficients between SCE from SF and MOD10A2 and lagged streamflow are also calculated for the Horsefly River subbasin (Fig. 6b).The maximum lagged correlation coefficients for the Horsefly River subbasin is about −0.76 (p<0.001) at the 0-day lag date.Therefore, the 32-day lag in the QRB arises in part from the deep snowpacks in the basin's headwaters (see Table 7) and the presence of Quesnel Lake that forms a large reservoir of water, attenuating the spring freshet signal in the QR (Fig. 1).The scatter plots between (normalized) SCE and (normalized) streamflow are shown in Fig. 7.The correlation coefficient between the SCE of SF and streamflow are higher than that between SCE of MOD10A2 and streamflow (−0.80 vs. −0.76(p<0.001))owing to the elimination of cloud coverage by the SF.The correlation coefficients between normalized SCE and normalized streamflow are about 0.04 higher than those between SCE and streamflow.Regression analyses and corresponding coefficients of determination (R 2 ) for the snow seasons demonstrate that linear regressions fit better the data than the logarithmic regressions.6).To explore the reasons for the earlier snow ablation, the monthly average air temperature and accumulated precipitation for different hydrological years are compared (Fig. 8).For snow onset seasons (September to December), the average air temperature in 2003 is the second highest with 0.3 • C lower than the average air temperature in 2000 and 0.3 • C −1.3 • C higher than the average air temperatures in the other years.Furthermore, the accumulated precipitation in 2003 is about 140 mm-250 mm lower than the mean accumulated precipitation over 2000-2007, resulting in less accumulated snow compared to other years.precipitation are always less than −0.1, demonstrating there is no significant relationship between them.However, there exists a strong relationship between average air temperatures and the SCF 50% , with the most significant anticorrelation coefficient at −0.847 (p<0.01) on 8 June.The scatter plots between average air temperature from 24 February to 8 June and SCF 50% and between SCF 50% and R 50% are shown in Fig. 9b and c, respectively.Linear relationships between them exist, with correlation coefficients of -0.847 (p<0.01) and 0.815 (p<0.01),respectively.However, the pair of values in 2003 falls out of the linear relationship.After exclud-ing these outliers, the correlation coefficient between average air temperature and SCF 50% reaches -0.993 (p<0.001) on 8 June and the correlation coefficient between SCF 50% and R 50% reaches 0.913 (p<0.001).According to the linear regression, the slope between SCF 50% and average air temperature (d(SCF 50% )/dT ) is 10 days ( • C) −1 .This implies that an aggregated rise of 1 • C above the seasonal average air temperature leads to a 10-day advance in reaching 50% SCF in the QRB.As such, this provides a measure of predictability for both SCF 50% and R 50% during spring.

Concluding discussion
MODIS snow products have been used to analyze the snow distribution in different areas such as the Tibetan Plateau, the Upper Rio Grande River Basin and the North Slope of Alaska (Pu et al., 2007;Klein and Barnett, 2003;Zhou et al., 2005;Déry et al., 2005).However, cloud coverage always prevents MODIS sensors from fully detecting the underlying snow cover that induces snow mapping errors.This study applied a SF method to reduce the cloud coverage to improve the accuracy of snow maps from MODIS in the subboreal mountainous QRB.The SF method reduced the average cloud coverage in the QRB from 15% for MOD10A2 to 9% for SF products.
There are limited evaluations of MODIS snow products in sub-boreal mountainous forests.The measurements at four in-situ stations in the QRB were adopted to evaluate the accuracy of MODIS snow products.In all cases, the snow depth is not measured under a forest canopy, but rather on a grassy site.However, in the 5 km 2 area around the four in-situ stations, forests account for 80-95% of the landcover types at a 1 km resolution.Therefore, the statistical results should Hydrol.Earth Syst.Sci., 13,[1439][1440][1441][1442][1443][1444][1445][1446][1447][1448][1449][1450][1451][1452]2009 www.hydrol-earth-syst-sci.net/13/1439/2009/   represent the OA in the sub-boreal forest area.Comparisons based on the three thresholds (k = 0 cm, k = 2 cm, or k = 5 cm) for the in-situ SD, which were set to classify the ground as snow or no snow, indicated that the threshold of k = 2 cm performed better than the threshold of k = 0 cm and k = 5 cm.
The OA for the four in-situ stations decreased with elevation.For the SF, the OA were 93.1%, 87.9%, 84.0%, and 76.5% for HLGL, Barkerville, BMM and YPE, respectively.The SF not only reduced cloud coverage but also increased the OA by about 2% compared to MOD10A2.The OA of MOD10A1 was lower by about 10% at high elevations such as BMM and YPE compared to SF.The lower OA in higher elevations was mainly caused by the UE defining snow as no snow.The OA of MODIS snow maps in sub-boreal mountainous forests are comparable to that in other regions with different land cover types (Pu et al., 2007;Zhou et al., 2005).Snow cover within sub-boreal mountainous areas such as the QRB imposes a strong influence on streamflow.Snow ablation processes are strongly controlled by air temperature.The correlation coefficients between streamflow and SCE based on the SF were more significant than those from MOD10A2, owing to the effective reduction of cloud coverage by the SF during snow ablation seasons.The correlation coefficient between spring streamflow and SCE of the SF is about −0.84 (p<0.001),indicating that the snow cover forms an important resource of freshwater in the QRB.There was a 32-day lag period between snow ablation and its greatest impact on streamflow at the QR gauge.Many studies have examined the contribution of snow ablation to the spring freshet in high-latitude or alpine watersheds such as the Mackenzie River Basin of northern Canada (Dyer, 2008), Upper Rio Grande River Basin of Colorado (Zhou et al., 2005), and the large Siberian watersheds (Yang et al., 2003).However, the streamflow and SCE in the sub-boreal forest such as the QRB has higher correlation coefficients compared to other areas.The hydrological cycle is complex, involving processes such as precipitation, surface water runoff, surface infiltration and groundwater storage, and evaporation.Although the terrain of the QRB is not highly permeable, groundwater may have an impact on streamflow of the Quesnel River, particularly during the low flow season (winter).However, the majority of the annual discharge occurs during the spring freshet when the groundwater has much less impacts on discharge (Burford et al., 2009).This paper focuses only on the relationships between surface water storage such as snow and streamflow.Therefore, it is beyond the scope of this study to truly evaluate the goundwater contribution to the overall streamflow.An aggregated 1 • C rise in seasonal air temperatures led to a 10-day advance in reaching SCF 50% in the QRB.There was also a strong linear relationship between SCF 50% and R 50% , with a correlation coefficient of 0.815 (p<0.01).In general, climate change will have multiple, nonlinear impacts on snow accumulation and snowmelt runoff.Since the QRB is located in the sub-boreal forest where snowmelt forms the majority of the source of spring runoff, its timing is highly sensitive to air temperature.However, the results should be combined with other climate change scenarios for precipitation, evaporation, and earlier onset of photosynthesis in other areas.Gafurov and Bárdossy (2009) eliminate all of the cloud cover in MODIS snow products in the Kokcha Catchment, Afghanistan; however, the approaches are based on different assumptions of the temporal and spatial distribution of snow cover, some of which are not well understood in the complex topography of the sub-boreal forest.The results of this study demonstrate that the SF is a feasible, straightforward method that reduces the cloud coverage to improve the snow mapping from the original MOD10A2 product.The spatial filtered MODIS snow products in sub-boreal mountainous forests provide a better understanding and prediction of the snow cover distribution and the resulting streamflow in these regions.The methods and results developed here improve the understanding of snow distribution with elevation (Tong et al., 2009) and of the resulting hydrological response over complex terrain.This leads to an enhanced ability to model and predict land surface process changes caused by climate variability and change over complex topography.

Fig. 1 .
Fig. 1.Geographical map of the Quesnel River Basin (QRB) and surrounding region.Stars represent the stations within the QRB whose air temperature, precipitation and SD are used.Triangles represent the stations around the QRB whose air temperature and precipitation are used.The green circles with numbers 1, 2, 3, and 4 represent the locations for the gauges of Horsefly (I.D. 08KH031), Moffat (I.D. 08KH019), Horsefly above McKinley (I.D. 08KH010), and McKinley (I.D. 08KH020), respectively.The Horsefly River subbasin is the area in the lowest red line polygon.
for the Horsefly River gauge (I.D. 08KH031) are estimated by summing the measured streamflow at Moffat Creek (I.D. 08KH019), McKinley Creek (I.D. 08KH020), and the Horsefly River above McKinley Creek (I.D. 08KH010).The boundary of the Horsefly River subbasin and locations of these four hydrometric gauges are shown in Fig. 1.

Fig. 2 .Fig. 2 .
Fig. 2. Flow chart describing the spatial filter method employed to reduce cloud coverage in MOD10A2.

Fig. 3 .
Fig. 3. Comparison of snow maps of MOD10A1 (top), MOD10A2 (middle), and the SF (bottom) in the QRB within the same period.Snow maps of MOD10A2 (middle) and SF (bottom) on 17 January 2007 represent the 8-day composite products from MOD10A1 covering the period from 17 January 2007 to 24 January 2007.
50% in different years ranges from its earliest value near 25 March 2005 and 27 March 2003 to the latest value around 10 May 2002, amounting to a difference of 40 days in reaching this mark (Table

Fig. 5 .
Fig. 5. Normalized accumulated runoff of QR and 8-day maximum SCF from MOD10A2 and SF in the QRB during snow ablation seasons from 2000-2007.

Fig. 6 .
Fig. 6.Lagged correlation coefficients (p<0.001) between 8-day maximum SCE from SF and MOD10A2 (a) in the QRB and streamflow of QR and (b) in the Horsefly River subbasin and streamflow of Horsefly River during snow ablation seasons from 2000-2007.

Fig. 9 .
Fig. 9. (a) Correlation coefficients and (b) scatter plots between average air temperature within different periods and Julian day of SCF 50% and (c) between Julian day of SCF 50% and R 50% during snow ablation seasons from 2000-2007 in the QRB.
August 2000.The 8-day average air temperature and 8-day accumulated precipitation measured at the nine insitu stations are calculated for hydrological years 2000-2007.During snow ablation seasons from 24 February to 24 June with 16, 8-day values, the corresponding average air temperatures for given periods are calculated.The periods are aggregated from the first date 24 February to another date; for example the period from 24 February to 4 March is represented by 4 March; period from 24 February to 8 June is represented by 8 June, as discussed later.

Table 3 .
The number of cloudy days for MOD10A1, MOD10A2, and SF snow products for 5×5 pixel 2 and total number of days when the in-situ observations are available.The numbers in parentheses are percentage of cloudy days compared to the total days.

Table 5 .
Statistics of the intercomparison between MOD10A2, SF and in-situ observations from 2000-2007 for the retrieval method of 5×5 pixel 2 .The numbers in parentheses are the times that the stations observed as snow or land, respectively.

Table 7 .
The annual maximum snow depth measured at the four in-situ stations (unit: cm)