Use of satellite-derived data for characterization of snow cover and simulation of snowmelt runoff through a distributed physically based model of runoff generation

A technique of using satellite-derived data for constructing continuous snow characteristics fields for distributed snowmelt runoff simulation is presented. The satellite-derived data and the available ground-based meteorological measurements are incorporated in a physically based snowpack model. The snowpack model describes temporal changes of the snow depth, density and water equivalent (SWE), accounting for snow melt, sublimation, refreezing melt water and snow metamorphism processes with a special focus on forest cover effects. The remote sensing data used in the model consist of products include the daily maps of snow covered area (SCA) and SWE derived from observations of MODIS and AMSR-E instruments onboard Terra and Aqua satellites as well as available maps of land surface temperature, surface albedo, land cover classes and tree cover fraction. The model was first calibrated against available ground-based snow measurements and then applied to calculate the spatial distribution of snow characteristics using satellite data and interpolated ground-based meteorological data. The satellite-derived SWE data were used for assigning initial conditions and the SCA data were used for control of snow cover simulation. The simulated spatial distributions of snow characteristics were incorporated in a distributed physically based model of runoff generation to calculate snowmelt runoff hydrographs. The presented technique was applied to a study area of approximately 200 000 km 2 including the Vyatka River basin with catchment area of 124 000 km 2. The Correspondence to: A. Gelfan (hydrowpi@aqua.laser.ru) correspondence of simulated and observed hydrographs in the Vyatka River are considered as an indicator of the accuracy of constructed fields of snow characteristics and as a measure of effectiveness of utilizing satellite-derived SWE data for runoff simulation.


Introduction
The spatial variations of snow characteristics play a significant role in the hydrological cycle of river basins and snowmelt runoff generation.Information on these variations can be used in distributed physically based models of runoff generation for simulation of the spatial peculiarities of hydrological processes and improving prediction and forecasting of snowmelt floods.However, the existing network of ground-based stations most often does not provide measurements of snow cover at the required spatial and temporal resolution.Spatial maps of snow cover derived from satellites provide a promising opportunity to enhance the assessment and monitoring of the spatial and temporal variability of snow characteristics, particularly in areas with a sparse network of meteorological stations.Reliability of these products and their spatial resolutions have noticeably improved during the last years.To a large extent this improvement was due to the launch of several new satellite observing systems with advanced capabilities such as the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard Terra and Aqua satellites and the Advanced Microwave Sounding Radiometer for Earth Observing System (EOS) satellites (AMSR-E) onboard Aqua (König et al., 2001;Hall et al., Published by Copernicus Publications on behalf of the European Geosciences Union.L. S. Kuchment et al.: Characterization of snow cover and simulation of snowmelt runoff 2002).However, the potential of satellite snow data is limited by a number of environmental factors (cloudiness, land cover type, terrain peculiarities, etc.) as well as by insufficient, in many cases, accuracy of satellite measurements.Clouds create discontinuity in the spatial distribution and in time series of snow data.Dense forest vegetation complicates snow cover identification and mapping due to snow interception and ability to mask snow cover on the forest floor.The accuracy of satellite measurements of snow water equivalent (SWA) significantly depends on the actual properties of the snow, especially the amount of liquid water.
A possible way to improve characterization of the snow spatial distribution and temporal variability consists in coupling satellite snow cover products with ground-based meteorological measurements and snow pack models.Most of the studies are focused at using satellite measurements of snow covered area (SCA) (e.g.Rodell and Houser, 2004, Andreadis and Lettenmaier, 2006, Dressler et al., 2006, Kolberg et al., 2006).US National Operational Hydrologic Remote Sensing Center (NOHRSC) has developed a Snow Data Assimilating System (SNODAS) where the ground-based, airborne and satellite snow observations were assimilated into the snow model to obtain the snow cover characteristics at 1 km spatial resolution and hourly temporal resolution (Carroll et al., 2006).A joint US Air Force/NASA blended, global snow product named "ANSA" utilizing both MODIS, AMSR-E and QuikSCAT sensor data has been developed and presented in (Foster et al., 2007).
In this paper we present a technique for constructing space-time continuous fields of snow cover characteristics (SWE, snow depth, snowmelt, etc.) on the basis of a physically based model of snow pack and with the use of satellite measurements of SWE.Satellite SCA products are used for control of the calculated snow cover fields and for calibration of the snow model.Satellite products also include maps of land surface temperature and albedo, which are utilized by the model when these maps are available.The model is first calibrated against available snow measurements at meteorological stations.Next, satellite-derived SWE maps specifically corrected for forested areas are utilized as the initial conditions and ground-based meteorological data as boundary conditions to simulate the spatial distribution of snow characteristics.Finally, the simulated spatial distributions of snow cover characteristics are incorporated in a distributed physically based model of runoff generation to calculate snowmelt runoff hydrographs.The correspondence of simulated and observed hydrographs may be considered as an indicator of the accuracy of the constructed fields of snow characteristics and at the same time, as a measure of effectiveness of utilizing satellite-derived SWE data as the initial conditions for runoff simulation.It is obvious a priory that because of possible large errors in satellite SWE maps such utilizing can be effective if the ground-based meteorological network is sparse and does not represent properly spatial heterogeneities of snow characteristics.
The current technique was applied to a study area of about 200 000 km 2 located in the European part of Russia with 56 • N-60 • N and 48 • E-54 • E spatial coordinates.The study area incorporates the Vyatka River basin with a catchment area of approximately 124 000 km 2 (Fig. 1), for which this technique was used for simulation of snowmelt runoff hydrographs.
The study area has flat terrain and mixed vegetation cover.In its northern part more than 80% of the area is covered by coniferous and mixed forests.The southern part is mostly agricultural land with less than 10-15% forest cover (Fig. 2).The snow season lasts for about 5 months with seldom thaws.Maximum snow water equivalent accumulates in the end of March and ranges from 100-120 mm in the south of the region to about 200-250 mm in the north.Snow melting starts, on average, on the 27 March and ends in the beginning of May.

Satellite and ground-based information used in the study
On the basis of measurements from the Advanced Microwave Scanning Radiometer (AMSR-E) of NASA's Earth Observing System (EOS) Aqua satellite, NASA issues since 2002 daily maps of SWE for the entire globe with 0.2 • grid pixel size of the latitude-longitude projection (AMSR-E/Aqua Daily L3 Global Snow Water Equivalent EASE-Grids -AE-DySno).These maps are computed using the empirical relationship between SWE and the difference in brightness temperature of the land surface at 18 and 36 GHz using appropriate fractions of forest cover and assuming vertical polarization.The potential error in SWE estimates using this satellite and retrieval algorithm is about 25% (Chang and Rango, 2000), yet for forested areas these errors can be significantly larger.Moreover, the accuracy of the SWE estimated from the radiometric satellite measurements noticeably decreases during melt period when snowpack is saturated by melted water (Engen et al., 2004).Significant errors can occur during thawing and when the land surface contains a thin ice crust.In addition, study by Dong et al. (2005) shows that the SWE retrievals are not sensitive to thin snow packs (SWE<10 mm).Results of validation of daily maps AE-DySno against ground snow courses for several test areas in the European part of Russia have shown that the satellitederived values of SWE before snowmelt can deviate as much as 200% from the actual ground SWE values (Nosenko et al., 2006).
On the basis of measurements from MODIS aboard NASA's Terra and Aqua satellites, NASA issues since 1999 daily maps of SCA (MOD10 L2 and MYD10 L2, respectively) with a resolution of 0.01 • in the latitude-longitude projection.Every pixel in the map is classified as snow, snow-free land surface or undetermined.The latter class includes pixels that were obstructed by clouds, did not get  enough daylight to be properly classified or did not pass any of data quality control tests.Earlier estimates of the accuracy of MODIS-based snow cover maps range from 90% to 98% depending on the season and the surface type (Hall and Riggs, 2007).Most studies report a substantial increase in snow detection error in forested areas (e.g.Simic et al., 2004).
The satellite information used in this study includes AE-DySno and MOD10 L2 maps as well as be-daily land surface temperature maps MOD11 L2 and albedo maps of albedo MOD43C1 derived at 5 km resolution from MODIS-Terra reflective channels data accumulated over 16-day periods.The accuracy of land surface temperature estimates is close to 1 • C.
Besides current MODIS and AMSR-E data, we also used static datasets such as type of land surface and the type and density of forest vegetation cover.All latter data products are based on observations from AVHRR aboard NOAA's satellite (Hansen et al., 2000).They have been produced at the University of Maryland Department of Geography and were acquired from http://glcf.umiacs.umd.edu/data.
The full list of the satellite products used in the present work is given in Table 1. Figure 3 shows examples of SWE, SCA, surface temperature and albedo maps.As can be seen from these maps, the satellite data reveal significant spatial heterogeneity of the considered variables.
Meteorological data were collected from 19 ground-based weather stations within the study area shown in Fig. 1.The data included 6-hourly observations of air temperature, air humidity, precipitation amount, atmospheric pressure, cloud cover and wind speed.Moreover, at the same stations, daily values of snow depth were made.

Snow pack model, its calibration and validation using ground-based data
The snow pack model (Gelfan et al., 2004;Kuchment and Gelfan, 2004) used in this study is a system of vertically averaged equations which includes description of temporal changes of the snow depth, contents of ice and liquid water taking into account snow melt, sublimation, refreezing melt water, and snow metamorphism: where H represents the snow depth; i and w signify the volumetric content of ice and liquid water, respectively; T s denotes the temperature of snowpack; S represents the snow melt rate; S i signifies the rate of freezing of liquid water in snow, E l denotes the rate of liquid water evaporation from snow; E s represents the rate of snow sublimation, Q a signifies the net heat flux at the snow surface; Q g denotes the ground heat flux; X s and X l represent the snowfall and rainfall rate at the snow surface, respectively; V signifies the snowpack compression rate; R denotes the snowmelt outflow from snowpack; c s represents the specific heat capacity of snow; ρ w , ρ i , and ρ 0 signify the densities of water, ice, and fresh-fallen snow, respectively; L denotes the latent heat of ice fusion.The snow melt rate, S, is computed from the energy balance of the snowpack at T s = 0 • C as: where Q sw represents the net short wave radiation; Q lw signifies the downward long wave radiation; Q ls denotes the upward long wave radiation from snow; Q T represents the sensible heat exchange; Q E signifies the latent heat exchange; Q P denotes the heat content of liquid precipitation.
Radiation fluxes are calculated differently depending on the properties of the vegetation cover.Components of Q a (in W m −2 ) for an open site are calculated as follows.The net short wave radiation and long wave radiation are expressed with the empirical relationships (Eqs.6-9) where the observed air temperature, air humidity, wind speed, precipitation, and cloudiness are used as input: where Q 0 =1000β represents the short-wave radiation flux under clear sky conditions for the day and the hour in question; β denotes the angle of short-wave radiation above the horizontal in radians, calculated as a function of the local latitude, the declination, and the hour angle; and α s signifies the snow albedo calculated from: in which α 0 represents an empirical coefficient; ρ s denotes the density of snowpack equal ρ s = ρ i i + ρ w w; N and N 0 signify the total and the lower level cloudiness (ratiometric), respectively; where σ represents the Stefan-Boltzmann constant (W m −2 K −4 ); e a denotes the air vapor pressure (mb), ε s signifies the effective emissivity of the snowpack taken equal to 0.99 in this study The turbulent fluxes of sensible and latent heat are calculated using: where r a represents the aerodynamic resistance; e s denotes the saturation air vapor pressure at the temperature of snow surface; ρ 0 signifies the air density; c a denotes the specific heat capacity of air; P a represents the atmospheric pressure; L s signifies the latent heat of sublimation of ice in the absence of liquid water or the latent heat of vaporisation when liquid water is present in snow; q T and q E represent the wind-less convection coefficients for the sensible and the latent heat fluxes, respectively.The coefficients q T and q E enable heat transfer to occur even when wind speed is zero.
L. S. Kuchment et al.: Characterization of snow cover and simulation of snowmelt runoff The aerodynamic resistance is calculated as follows: where U z represents wind speed at the height z; z 0 denotes the snow surface roughness assigned as 0.005 m; κ represents von Karman's constant, and Ri is the Richardson number whose value is estimated using: where g is the acceleration due to gravity.
The heat flux Q P caused by the liquid precipitation is expressed as: where c w represents the specific heat capacity of water.
The ground heat flux Q g is found from Eq. ( 28) describing the vertical heat transfer in soil.
At the forest floor, Eqs. ( 8) and ( 9) are modified to take into account the effect of canopy coverage and the type of vegetation on radiation fluxes.Net short wave radiation (Q * sw ) and downward long wave radiation (Q * lw ) fluxes on sub-canopy snow surface are calculated as: where C c represents the canopy coverage (ratiometric);k sw denotes the transmissivity through the canopy; Q lc specifies the long-wave radiation emitted by the canopy (upward and downward), calculated as ε c σ T 4 c , where T c represents the temperature of canopy ( • K) and assumed equal to air temperature; ε c denotes the emissivity of the canopy taken as 0.96.
Transmissivity is calculated using: where LAI denotes the leaf area index; Q ext represents the extinction efficiency estimated as: When calculating the turbulent heat exchange under forest canopy, it was assumed that air temperature and air humidity in forest are the same in the forest and in the open area adjoining to the forest.The wind speed in forest U * z is defined as: where k u denotes the coefficient of wind shield calculated as: It was assumed that the precipitation type changes from liquid to solid at an air temperature of 0 • C. Depending on the snowfall rate and the forest type part of the snowfall may be intercepted by the forest canopy.For a single snowfall onto a snow-free canopy, the amount of snow intercepted by the forest canopy is defined as: where X * s represents the snowfall rate above the canopy; I max = S p LAI 0.27 + 46ρ −1 0 represents the interception capacity, found as a function of the species snow loading coefficient, S p .The density of fresh-fallen snow ρ 0 (in kg m −3 ) depends on the air temperature.
The rate of sublimation of intercepted snow E i is derived from the energy balance of the intercepted snow as: where R i net represents the net radiation absorbed by the intercepted snow; Q i T denotes the sensible heat exchange for the intercepted snow.
The net radiation, R i net , absorbed by the intercepted snow is expressed as: 22) where α c represents the canopy albedo taken equal to 0.12 in this study.
The component Q i T is calculated from the ice bulb temperature and wind speed (Gelfan et al., 2004).
The rate of liquid water freezing in snowpack S i is calculated as: The snowpack compression rate V is found from (in cm s −1 ): where the density of snowpack, ρ s , is in g cm −3 ; v 1 , v 2 , and v 3 specify the coefficients equal 2.8×10 −6 cm 2 s −1 g −1 ; −0.08 • C −1 ; 21 cm 3 g −1 , respectively.The outflow of liquid water from snow is calculated as: where w max represents the holding capacity of snowpack related to its density ρ s as: Equations (1-4) were numerically integrated by the implicit finite-difference scheme.Numerical experiments were carried out to assess influence of time step of integration on the Hydrol.Earth Syst.Sci., 14, 339-350, 2010 www.hydrol-earth-syst-sci.net/14/339/2010/ simulation results.It was shown that decreasing of time step smaller than 6 h does not lead to increasing of simulation accuracy.As a result, time step of integration was 6 h, i.e. the same as the time-step of the input meteorological data.The snow pack model was calibrated and validated using available measurements of snow depth at 19 meteorological stations.To calibrate the model we used snow depth data obtained during the time period from 1 November 2001 to 30 May 2002.The values of α 0 (Eq.7), q T (Eq.10) and q E (Eq.11) were calibrated against snow depth and found equal to 1.01, 0.98 J m −2 s −1 and 0.11 J Pa −1 s −1 , respectively.Monte Carlo simulation was used to randomly generate one hundred parameter sets, from which optimal parameter set was identified.The parameters were assumed to be uniformly distributed within the pre-determined intervals.The snow depth measurements at the same meteorological stations during the time period from 1 November 2002 to 30 May 2005 were used to assess the model performance.The standard deviations of differences (errors) between calculated and observed snow depths were estimated for each of 19 meteorological stations.The overall mean standard deviation of errors was equal to 7 cm for the calibration period and increased to 9 cm for the validation period.As an illustration of the obtained results, Fig. 4 compares the results of snow depth simulations at several meteorological stations for the period of 1 November 2002 to 30 May 2003.
In order to test the validity of the model calibrated against the snow depth measurements in predicting SWE, we used the snow survey observations within the Vyatka River basin for 9 snow seasons: from 1971-1972 to 1979-1980.Both snow depth and SWE observations are available for that period.The model demonstrated satisfactory accuracy in predicting SWE: the overall mean error is 1.2 mm with mean standard deviation of about 10.6 mm.

Modeling spatial fields of snow characteristics
The calibrated and validated snow pack model was applied to simulate the spatial distribution of snow pack characteristics over the study area during the time period from 1 March, when the SWE is typically close to the maximum value and errors of its satellite measurements have to be minimal, to 30 June of 2003-2005.The model was run for each 0.01 • grid cell within the study area.To obtain the input data for the snow pack model, the daily ground-based meteorological observation data were interpolated to each pixel using the inverse distance squared method.
In the grid cells where the MODIS-based surface temperature and albedo data were not available, we calculated surface temperature and albedo using Eqs.( 4) and ( 7), respectively.
The initial spatial distribution of SWE for the open, i.e. forest-free, areas was assigned using AE-DySno SWE maps on 1 March.The SWE values for all forested pixels as well as the missed SWE values for open pixels were estimated by interpolation from open pixels where satellite-derived SWE values were available.The interpolated SWE was further multiplied by the factor k snow , representing the average ratio of the pre-melt SWE in the forest to the pre-melt SWE in the neighboring open area.The value of k snow was obtained from numerical simulation of snow accumulation in the forest with the snow cover model driven by meteorological observations during 14 snow seasons from 1970-1971 to 1979-1980 and from 2000-2001 to 2003-2004.The simulation results have shown that in the deciduous forest the value of k snow is close to unity.In coniferous forest this value changes with the canopy coverage.In the sparse coniferous forest, the pre-melt SWE exceeds the SWE in the open area by 9-11% because of smaller snow evaporation and snowmelt during a period of snow accumulation.At the same time, in the dense coniferous forest, interception and sublimation of intercepted snow during this period is larger and thus the pre-melt SWE L. S. Kuchment et al.: Characterization of snow cover and simulation of snowmelt runoff is 6-13% less to open ones.The numerical simulations have shown that the variations of the value of k snow from year to year for the same forest area are small and can be negligible.
Assuming that the spatial variability of snow density is much less than variability of SWE, the initial snow density on 1 March was accepted as constant for the whole area and defined from the measurements at the Kirov meteorological station.The initial volumetric moisture content of snow was assumed equal to zero.
Using the assumptions formulated above, we simulated daily maps of SCA and compared them with NASA MOD10 L2 SCA maps for the snowmelt season of the years 2002 and 2004.The comparison was carried out in order to refine the snow model parameters adjusted against the ground-based point measurements (see previous section) and use the refined parameters for further characterization of snow fields.The region was divided on to 19 Thissen polygons according to the location of the meteorological stations.Simulated, SCA calc , and satellite-derived, SCA sat , values were estimated for each polygon and for the dates when most of the polygons were free of cloudiness.The value of SCA was calculated as the number of open pixels covered by snow divided by the total number of open pixels within the polygon.When calculating SCA sat , only free of cloudiness pixels were taken into account.Two criteria were applied to summarize the goodness of fit of the simulated and satellite snow maps for each selected date: (1) the normalized mean error where SCA sat denotes the satellite SCA estimated for the whole polygon.It was appeared that minimum values of the both criteria are achieved under almost the same values of the snow model parameters (namely, α 0 =1.03, q T =0.98 J m −2 s −1 , and q E =0.12 J Pa −1 s −1 ) as the values obtained through the model calibration against the ground-based point measurements.
As an illustration of the obtained results, temporal changes of the criteria NME and NRMSE are shown in Fig. 5 for 3 of 19 polygons.
One can see from Fig. 5 that NME and NRMSE are close to zero in the beginning of spring, then the both criteria increase in the period of intensive melt and return to small values in the end of melt season.In general, the model allowed us to reproduce temporal changes of SCA for the open areas with satisfactory accuracy.
The accuracy of this reproduction is decreasing if we include forested pixels in addition to open ones when SCA calculating.Figure 6 shows the temporal change of SCA for two polygons which are located correspondingly in the south-eastern and north-eastern corners of the study area and have the coniferous forest percentage of 9% and 76% respectively.It can be seen from this Fig. 6 that, for the sparsely forested area the calculated values of SCA are close to the values defined from the satellite data while for the area with dense forest there is a significant difference between these values.
Daily maps of simulated snow characteristics were constructed for the time period from 1 March to 30 June of 2002-2005.Figure 7 presents the maps of simulated distributions of SWE for three dates in the second half of April 2003 corresponding to the period of intensive snowmelt.These maps apparently present the spatial picture of snowmelt dynamics.

Using spatial snow characteristics for distributed modelling of runoff generation in the Vyatka River basin
The simulated fields of snow characteristics were used in the physically based distributed model of runoff generation in the Vyatka River basin (Kuchment et al., 2008).The  model is based on the finite-element schematization of the basin and describes the processes of snow cover formation and snowmelt, freezing and thawing of soil, vertical moisture transfer and evaporation, surface water detention, overland, subsurface and channel flow.In the Vyatka River basin 477 area and 84 channel finite elements were defined taking into account the basin topography, soil properties, and vegetation type distribution as well as river channel structure and stream gage allocation (Fig. 8).
To calculate the characteristics of snow cover, Eqs. ( 1)-( 26) were applied vertical water and heat transfer in the soil associated with soil freezing, thawing and infiltration of water is described with the following equations (Gelfan, 2006): where W , θ and I represent the total water content, liquid water content and ice content of soil, respectively (W = θ + ρ i ρ w I ); K = K(θ, I ) denotes the hydraulic conductivity of the frozen soil; T specifies the soil temperature; λ represents  The capillary potential ψ = ψ (θ,I ) and the hydraulic conductivity K = K(θ, I ) of the frozen soil are determined from the relations proposed in Gelfan (2006).Equations ( 27)-( 28) were numerically integrated by an implicit, four-point finite difference scheme.The corresponding difference equations were solved by the double-sweep method.At the lower boundary of the podzol soils typical for the Vyatka River basin, there is an impermeable layer at the depth of 1-2 m.The vertical water flux is assumed zero at this boundary.It is also assumed that the horizontal movement of water along the impermeable layer occurs if the soil moisture content exceeds the field capacity, θ f , of soil.As a result, the horizontal flux, R g , is calculated as: where z P is the soil layer in which θ > θ f .The detention of melt water by depressions at the catchment surface is calculated by the formula assuming exponential distribution of the storage capacity (Kuchment et al., 1986).
The rate of evaporation from an unfrozen, snow-free soil is calculated by the formula presented by Kuchment et al. (2000).
The kinematic wave equations are applied to describe overland and subsurface flow.To account for the subsurface flow, the following equations are used (Kuchment et al., 2000): where h is the subsurface flow depth; q is the subsurface flow discharge; i 0 is the slope of the layer with subsurface flow; K g is the horizontal hydraulic conductivity.
The kinematic wave equations are numerically integrated by the finite-element method coupled with the Galerkin method.
To calculate the water movement through the river channel elements, the advection-diffusion equation is applied.This equation was numerically integrated by the four-point implicit difference scheme.
To take into account subgrid effects, it was supposed that the saturated hydraulic conductivity is gamma-distributed within each finite element area.Mean value of this parameter within a finite element was assessed using the available experimental measurement data for the different types of soil in the Vyatka River basin.Coefficient of variation within a finite element was determined from the empirical formula (Kuchment et al., 1986) which relates this coefficient to the mean value of the saturated hydraulic conductivity.
Most of the model parameters were assigned on the basis of the available measurements of the basin characteristics including the basin topography and river channel data, soil, vegetation and snow constants as well as from empirical relationships that were derived and tested in (Kuchment et al., 2008) using mainly Russian experimental data and field observations.Six parameters influencing the processes of infiltration, soil moisture evaporation, detention in basin storage and flood routing were calibrated against the observed snowmelt flood hydrographs for the period from 1940 till 1959.Two parameters influencing snowmelt rate were calibrated against snow measurements.The validation was carried out by comparison of the observed and simulated hydrographs for the period from 1960 till 1980.The simulated and the observed hydrographs at the Vyatskie Polyany gage (the outlet gage of the Vyatka River basin) for the last ten years of the validation period are presented in Fig. 9.The standard deviation of errors of the simulated flood volumes and peak discharges are equal to correspondingly 1.1 km 3 and 486 m 3 s −1 .The Nash and Sutcliffe efficiency criterion for the flood volume and discharge simulations are 0.94 and 0.84, respectively.
The calibrated model of runoff generation in the Vyatka River basin was applied to simulate the snowmelt floods of 2003 and 2005 using the constructed fields of snow characteristics obtained by the technique presented in the previous section.The calculated and observed hydrographs for three gauges of the Vyatka River: Kirov, Kotelnich, Vyatskie Polyany are compared in Fig. 10.It can be seen from Fig. 10, that the simulated fields of snow characteristics have provided the satisfactory accuracy of hydrograph simulation.The flood volume and the peak discharge errors do not exceed 20% at the basin outlet.As a numerical experiment, we calculated snowmelt fields directly from the AE-DySno SWE maps and used these fields as inputs into the runoff generation model.It appears that the resulting hydrographs shown in Fig. 10 are significantly underestimated in comparison with the observed ones.This underestimation can be explained by large errors in satellitederived SWE in forested areas and for snowpack saturated with melt water.
The obtained results can be viewed as an indication that characterization of snow cover fields by the proposed technique can improve the representation of spatial distributions of SWE as compared to the respective fields obtained directly from satellite data.There is also an opportunity to further improve the accuracy of the proposed method using a more comprehensive calibration of the snowpack model with satellite measurements of SCA (if these data series are long enough) and the observed runoff hydrographs.

Conclusions
In this paper, we have presented a technique that uses satellite remote sensing data to construct fields of snow characteristics, which are continuous in time and space and to be used for distributed snowmelt runoff simulations.Within this technique, satellite land surface monitoring data and available standard ground-based meteorological measurements are incorporated in a physically based model of snow pack formation.The model provides adoption of available satellite data including snow water equivalent, fractional snow-covered area, albedo, snow surface temperature, types of the surface land and takes explicitly into account the heterogeneities of the river basin and spatial variations of its characteristics.As compared to snow products derived solely from satellite data the advantage of this approach to snow characterization consists in a detailed physical description of snow processes which includes the forest cover effects.This allows for obtaining additional information on snow cover that is absent in satellite measurements.
The proposed technique may be useful for runoff simulation in river basins where the ground-based meteorological network is sparse and does not represent properly spatial heterogeneities of snow characteristics.Because of large possible errors in satellite measurements of SWE during a snowmelt period, the satellite-derived SWE are used as the initial conditions.The correspondence of simulated and observed hydrographs may be considered as an indicator of the accuracy of the constructed fields of snow characteristics and at the same time, as a measure of effectiveness of utilizing satellitederived SWE data for runoff simulation.It is possible to assume that the suggested approach is especially designed to help physically based distributed hydrologic models realize their full potential in predicting spatial peculiarities of river runoff genesis.

Fig. 1 .
Fig. 1.Location of the study region (red points -meteorological gauges; blue points -hydrological gauges).The Vyatka River is highlighted in yellow.

Fig. 1 .
Fig. 1.Location of the study region (red points -meteorological gauges; blue points -hydrological gauges).The Vyatka River is highlighted in yellow.

Fig. 2 .Fig. 2 .
Fig. 2. Land cover classification within the study region Fig. 2. Land cover classification within the study region.

Fig. 3 Fig. 3 .
Fig. 3 The satellite-derived maps of SWE, SCA, land surface temperature and albedo for the study area

Fig. 4 .Fig. 4 .
Fig. 4. Seasonal change of the observed (points) and simulated (line) snow depths (in cm) at selected stations within the study area for the season from November 2002 to June 2003 Fig. 4. Seasonal change of the observed (points) and simulated (line) snow depths (in cm) at selected stations within the study area for the season from November 2002 to June 2003.

Fig. 5 Fig. 5 .
Fig. 5 Normalized mean error (NME) and normalized root mean square error (NRMSE) of simulated snow covered area in comparison with one obtained from NASA MOD10_L2 maps for three Thissen polygons surrounding meteorological stations Kilmez, Kirov, and Kirs (forested pixels are not taken into account ).01/03/2002-30/04/2002 Fig. 5. Normalized mean error (NME) and normalized root mean square error (NRMSE) of simulated snow covered area in comparison with one obtained from NASA MOD10 L2 maps for three Thissen polygons surrounding meteorological stations Kilmez, Kirov, and Kirs (forested pixels are not taken into account).1 March 2002-30 April 2002.

Fig. 8 .
Fig. 8. Finite-element schematization of the Vyatka River basin (bold lines represent the channel network; thin lines represent the boundaries of finite elements; green labels mark the forested elements).

Fig. 8 .
Fig. 8. Finite-element schematization of the Vyatka River basin (bold lines represent the channel network; thin lines represent the boundaries of finite elements; green labels mark the forested elements).

Fig. 9 .
Fig. 9. Comparison of the observed (solid line) and calculated (dashed line) hydrographs at the gauge Vyatskie Polyany of the Vyatka River (the last ten years of validation period).

Fig. 9 .
Fig. 9. Comparison of the observed (solid line) and calculated (dashed line) hydrographs at the gauge Vyatskie Polyany of the Vyatka River (the last ten years of validation period).

Fig. 10 .Fig. 10 .
Fig. 10.Comparison of the observed hydrographs (bold line) with the hydrographs calculated from the constructed SWE fields (thin line) and from AE-DySno SWE maps (dashed line).(The left column presents hydrographs of the 2003 flood, the right column presents hydrographs of the 2005 flood) Fig. 10.Comparison of the observed hydrographs (bold line) with the hydrographs calculated from the constructed SWE fields (thin line) and from AE-DySno SWE maps (dashed line).(The left column presents hydrographs of the 2003 flood, the right column presents hydrographs of the 2005 flood).

Table 1 .
List of satellite data products used in the study.