Simulation of snow distribution and melt under cloudy conditions in an Alpine watershed

An energy balance method and remote-sensing data were used to simulate snow distribution and melt in an alpine watershed in northwestern China within a complete snow accumulation-melt period. The spatial energy budgets were simulated using meteorological observations and a digital elevation model of the watershed. A linear interpolation method was used to estimate the daily snow cover area under cloudy conditions, using Moderate Resolution Imaging Spectroradiometer (MODIS) data. Hourly snow distribution and melt, snow cover extent and daily discharge were included in the simulated results. The root mean square error between the measured snow-water equivalent samplings and the simulated results is 3.2 cm. The Nash and Sutcliffe efficiency statistic (NSE) between the measured and simulated discharges is 0.673, and the volume difference (Dv) is 3.9 %. Using the method introduced in this article, modelling spatial snow distribution and melt runoff will become relatively convenient.


Introduction
Snow is a very important component of the climate system because of its significant effect on surface energy and water balances; thus, its accurate representation is essential for a better understanding of climate effects on the hydrological cycle (Tarboton et al., 2000).Rather than using field measurements, remote sensing (RS) is now more widely used in modelling snow because of its unique advantages in obtaining large-scale ground information (Seidel and Martinec, Correspondence to: H.-Y. Li (lihongyi@lzb.ac.cn) 2004).The snow cover area (SCA) can be observed with higher spatial resolution at near-infrared and visible wavelengths (Hall and Riggs, 2006).The snow-water equivalent (SWE) distribution can be simulated by combining the energy balance method and the SCA data at each time step (Cline et al., 1998).The spatial and temporal variations in the snow are fundamentally driven by the energy and mass balance.If the energy inputs for snowmelt at each time step that snow was present could be assessed and it is considered that the snow duration (T ) is a sum of all time steps for which snow existed, then the SWE change ( SWE) in a complete snow accumulation-melt period could be regarded as a function of the snow duration (T ) and the energy input (E).

SWE = f (E, T )
(1) With a continuous observation of the SCA and accurate simulation of the spatial distribution of energy budgets, the snow duration (T ) and the energy input (E) for snowmelt at each grid could be obtained, and the snow processes could be theoretically calculated.However, the key to this method is obtaining snow existence durations under cloudy conditions.Cline et al. (1998) had simulated the spatial energy balance and obtained the total SWE distribution by this method, but only a few satellite images were used.Snowfall, cloudy conditions and presnowmelt season were ignored in the simulation.The assessment of snow existence under clouds is a known problem.For instance, the Moderate Resolution Imaging Spectroradiometer (MODIS) data that provide daily SCA information are used widely in mountainous SWE reconstruction, land surface models and hydrological models (Zaitchik and Rodell, 2009;Parajka and Bloeschl, 2008;Durand et al., 2008) because of the high temporal and spatial resolution and accurate identification of snow.However, the main disadvantage of the MODIS sensor is that it is unable to record observations in cloud-covered regions (Gafurov and Bardossy, 2009).Thus, some researchers had to discard scenes seriously affected by clouds (Su et al., 2008) or use the eight-day composite maximum snow cover tile product as a replacement (Tekeli et al., 2005).However, even in those images with a large area of cloud cover, many snow pixels were not affected.If those images are not used, then abundant data are wasted.It has been inconvenient to accurately represent snow accumulation and melt processes.In this paper, a linear interpolation method was developed to avoid the influence of clouds on SCA data and to calculate SWE change by utilising MODIS data.
The goal of this paper is to model continuous SWE and melt runoff by an energy balance method under cloudy conditions in an alpine watershed.The method consists of three steps: (1) estimating daily snow existence using MODIS data and meteorological observation; (2) modelling daily spatial energy balance processes and calculating the total SWE distribution using the estimated SCA, digital elevation model (DEM) and meteorological observation; and (3) calculating the daily SWE distribution based on calculated hourly snowmelt and snowfall observations.

Site description
The Binggou watershed, located upstream of the Heihe river basin in northwestern China, was chosen as one of the focal experimental areas within the WATER (Watershed Allied Telemetry Experimental Research) framework (Li et al., 2009).The altitude range of the Binggou watershed is from 3440 to 4400 m, and the surface area is approximately 30.27 km 2 .It is a seasonally snow-covered region.The mean depth of the seasonal snowpack is approximately 0.5 m, up to a maximum of 0.8-1.0 m.The snow redistribution is remarkable because of the interaction between blowing snow and complex terrain.There is more snowfall in spring and autumn than in winter.A snowfall period takes place from November to April, a rainy season occurs from May to August and snowfall occurs again after September.The mean air temperature at an altitude of 3450 m was approximately −2.5 Snow characteristic data were collected through a series of field measurements in the snow season of 2008.These datasets include basic snow properties, such as snow depth, density, grain size, temperature, emissivity and so on.The topographic information, such as elevation, slope and the aspect of the slope, were all obtained from the DEM of the Binggou watershed, with a 50 m resolution.

Assessment of snow existence under cloudy conditions
MODIS snowcover products, MODIS/Terra Snow Cover Daily L3 Global 500 m Grid (MOD10A1), were collected.Version 005 of MOD10A1 contains a new fractional snowcover product developed from a linear fit of binary Thematic Mapper snowcover to the normalized difference snow index of MODIS (Painter et al., 2009).Data fields in MOD10A1 data, SCA and snow cover fraction (SCF), were used to model snow distribution.The MODIS data were resampled to match the 50 m resolution of the DEM.The spatial distribution of clouds is largely different between daily MOD10A1 scenes, while SCAs and positions changed relatively slowly.If the duration of cloudiness is short enough within a grid of the study area, then accurately estimating snowcover fractional values for the grid on cloudy days by interpolation of a series of cloudless SCF values for the grid is possible.
A simple linear method was used to obtain the daily SCF in the study.Considering snowpack evolution at each grid of the MOD10A1 scene, SCA depletion is continuous in each grid, except when snowfall occurs.Assuming a cloudy period from 1 to n − 1 days and that 0 to n days are cloudless, the SCF value of day x [SCF(t x )] in this cloudy period can be obtained by a linear interpolation method, according to Equation (2): where SCF(t 0 ) is the SCF value of the day 0 and SCF(t n ) is the SCF value of day n.Because Eq. ( 2) cannot hold for any number of cloudy days, 4 days was set as the threshold.

Methods
A distributed physical snow model was designed to model daily snow distribution and melt.Meteorological data, MODIS snowcover data and the DEM were the major inputs for the model.Spatial meteorological elements were calculated using continuous hourly meteorological observations and the DEM of the watershed.Continuous SCF datasets were used to simulate daily snow ablation by a distributed energy balance method.Finally, the total SWE change over the entire snow season was obtained.The daily snow distribution was calculated based on the total SWE change, the daily snowfall records and the calculated daily snowmelt.

Spatial SWE calculation
At the grid, the SWE at time t (M grid ) is estimated from where P is the snowfall water equivalent rate (m s −1 ), S is the grid area (m 2 ), R grid is the melt-water output (m s −1 ) at the grid and E grid is the snowpack mass loss of sublimation (m s −1 ).Assuming sublimation and melt are homogeneous in a grid, where R is the melt runoff rate at the point (m s −1 ), E is the snow sublimation rate (m s −1 ) and SCF is the snow fraction proportion.

Energy balance method
The energy budget for the snowpack can be written as where Q melt is the energy budget for the snowmelt outputs (J) and Q i is the internal energy change of the snowpack (J).S net , L net , H , LE, G, and Q p are the net shortwave radiant energy input (J), the net longwave radiant input (J), the sensible heat(J), the latent heat of sublimation (J), the ground heat conduction (J) and the precipitation sensible and latent heat (J), respectively.Lf is the latent heat of fusion (J kg −1 ).Before melt runoff can be released from the snowpack, the average temperature of the snowpack must be raised to the melting point and its maximum liquid-water-retaining capacity must be reached.Further net energy inputs (Q melt ) produce melt-water outputs: where ρ w is the liquid water density (kg m −3 ) and B is the thermal quality of the snow (0.97).

SWE adjustment with MODIS snowcover data
The existence of snow is one of the necessary conditions for snowmelt to occur.In modelling the SWE, the snow depth was also modelled for the energy exchange computation at each time step.The modelled snow depth of a grid is possibly 0, but there is actually snow in the calculated SCA image because of external factors, such as blowing snow and internal errors of the model.In this case, the SWE was estimated from the SCF by an empirical method (Yang et al., 1997;Essery and Pomeroy, 2004): where M avg is the average SWE of the grid, M ini is the premelt SWE of the grid and Cs is the coefficient of variation.
The distribution of snow depth was measured in a few experimental regions in the Binggou watershed in the 2008 snow season.We noticed that land surface can be covered fully by a snowpack with an average depth of approximately 20 cm in these regions.M ini is then assumed to be 20 cm, and the Cs = 0.95, where 0.95 is an empirical value.
Because the snowpack is regarded as a single layer, an iterative algorithm that has been used in the UEB model (Tarboton, 1996) was used to model the snow surface temperature.First, the model was calibrated for the point scale at the DY station without using RS data.The calibrated parameters are shown in Table 1.After calibration, these parameters were used to simulate the snow processes of the watershed.

Simulation of discharge
A simple conceptual method was used to estimate discharge according to the following equation, written for a time lag between the daily melt cycle and the resulting discharge cycle of 6 h (Martinec et al., 2008): where D is the average daily discharge (m 3 s −1 ), "melt" is the total snowmelt water equivalent (m), "rain" is the daily rainfall depth (m), n and n + 1 are the sequence of days during the snow season, C snow is the runoff coefficient of snow and C rain is the runoff coefficient of rain.A l is the conversion constant from m m 2 d −1 to m 3 s −1 .k is the recession coefficient indicating the decline of discharge in a period without snowmelt or rainfall (Martinec et al., 2008).The statistical behaviour of the discharge decline in periods without snowmelt and rainfall was analysed to obtain the recession coefficient k.The results can be formulated as The Nash and Sutcliffe efficiency statistic, NSE (R 2 ), and the volume difference, Dv (%), were used to assess the accuracy of the simulation results.
where D i is the daily discharge measured, D i is the daily discharge computed and D is the average measured discharge of the simulation period.Dv is the percentage difference between the total runoff (%) measured and simulated, V is the seasonal runoff volume measured, and V is the seasonal runoff volume computed.

Wind speed, air temperature, humidity, precipitation
In the model, hourly air temperatures were linearly interpolated between the two meteorological stations as a function of elevation.The wind speed values were modified according to topographic slope and curvature relationships, as Liston and Elder (2006) suggested.The observed relative humidity was converted to specific humidity, and then, the spatial relative humidity was extrapolated using the specific humidity and the previously extrapolated air temperature (Cline et al., 1998).Precipitation data were extrapolated to the altitude of each grid by an altitude gradient of approximately 1 %, following previous observations (Yang et al., 1992).The critical temperature assumed to determine the precipitation form (snow or rain) is 1 • C.

Solar radiation
The total incoming solar radiation observed was separated empirically into direct and diffuse components (Bristow et al., 1985).The two components were distributed to the whole watershed.
where τ d is the atmospheric diffused transmission coefficient, τ t is the atmospheric total transmission coefficient and B is the maximum clear-sky transmissivity of the atmosphere (0.76) (Flerchinger, 2000).τ d is assumed to be the same over the small watershed (30.27 km 2 ) on an elevation level before solar radiation reaches the ground.By Eq. ( 14), the solar radiation was divided into diffused transmission and direct transmission.Direct and diffused solar radiation all provide energy inputs for snowmelt.However, these two types of radiation distributed energy with different behaviours.Direct radiation was influenced by terrain slope and aspect, while diffused radiation was influenced mainly by sky-view factors.Complex alpine terrain modified the exchange of direct and diffused solar radiation.The incoming solar radiation on a slope can be computed according to the following equations (DeWalle and Rango, 2008): S diff slope = S diff cos 2 k s /2 + (1 − cos 2 k s 2) α t S envi ( 16) where S slope is the incoming shortwave radiation flux density on the slope, S dire is the incoming direct shortwave radiation at normal incidence, S diff slope is the incoming diffuse shortwave radiation flux density to the slope, Z is the angle between the solar beam and a perpendicular to the slope, S diff is the diffused shortwave radiation flux density on a horizontal surface, k s is the slope inclination angle and (cos 2 k s /2), (1−cos 2 k s /2) represent the sky-view factor and the fraction of the view from the slope occupied by the surrounding terrain, respectively.S envi is the total shortwave radiation flux density from the surrounding terrain.Albedo was calculated as the average of two reflectances in the visible and near infrared bands.The reflectance in each band is a function of the snow surface age and the solar illumination angle (Dickinson et al., 1993;Tarboton, 1996).When the snow depth is less than 10 cm, the albedo is taken as the composition of the ground albedo and the snow surface albedo (Dickinson et al., 1993).When solar radiation entered the snowpack, the net radiation flux at different depths was calculated with empirical equations (Anderson, 1976;Flerchinger and Cooley, 2000).

Longwave radiation and turbulent heat exchange
Outgoing longwave radiation was computed from the Stefan-Boltzmann equation using the snow surface temperature (T ss , • C).Incoming longwave radiation was estimated based on the air temperature (T a , • C): (17) where ε s is the snow emissivity (0.99), ε ac is the air emissivity under cloudy conditions and σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W m −2 K −4 ).ε ac is estimated from the following (Bristow et al., 1986;Campbell, 1985;Flerchinger, 2000): where ε a is the clear-sky air emissivity, C is the cloud cover fraction, a ε and b ε are empirical coefficients (Idso and Jackson, 1969) and τ t is the atmospheric transmissivity calculated from the shortwave radiation measured.The terrain influence on incoming atmospheric longwave radiation was modified by the factors (cos 2 k s /2) and (1−cos 2 k s /2), similar to Eq. ( 16) (DeWalle and Rango, 2008).The bulk aerodynamic approach was used to compute the sensible and latent heat exchange (Hood et al., 1999).Rainfall energy was also computed in the model (US Army Corps of Engineers, 1998).Ground heat conduction was assumed to be a small, constant value for heat flow to the snowpack.

Daily snowcover distribution
By the linear interpolation method, daily SCF images without clouds were obtained.To prove the effectiveness of this method, the frequency of different cloudy durations in the snow season of 2008 was counted.
For the illustration of the accuracy of the linear interpolation for gap-filling the snowcover time series, we sampled those grids with continued snow existence.It was assumed that those grids were obscured by clouds during different days, and the SCF values of these grids were linearly interpolated.The interpolation results were compared with the real observed SCF data, and mean errors (absolute value) corresponding to different gap-filling days were computed (Fig. 2).This exercise indicated that the mean errors are stable and small (approximately 5 %) in those situations for which the number of gap-filling days is less than 5.However, the mean errors are not stable and acutely change when the number of gap-filling days is greater than 8. Fortunately, the most frequent duration found was 1 day; this duration value accounted for 53.9 % of the total samples.The maximum cloudy duration was 15 days, a length observed for only 0.8 % of the samples, and a duration less than 4 days accounted for 91.0 % of the samples.This suggests that more than half the cloudy durations were only 1 day and the rest were less than 4 days.The shorter the cloudy duration was, the more reliable the interpolation results were.As such, the interpolation results were theoretically close to the actual conditions.

SWE distribution and change
Spatial energy inputs are the most important factors for snowpack melting.The net energy inputs to the snowpack above the whole watershed were negative before snowmelt occurrence (Fig. 3).Shortwave radiation accounted for the major part of the total energy inputs.After 17 March, the net energy inputs became positive, with increasing solar radiation and air temperature.Corresponding to the changes in the energy inputs, there were three marked snowmelt processes.
The total spatial SWE change at the grids was calculated using the SCF data calculation and the energy balance method, the distribution of which has marked topographic features.The snow distribution demonstrated strong heterogeneity because of its redistribution by blowing snow and inhomogeneous precipitation.The SWE distribution had an increasing trend with higher altitudes.The daily snowmelt and discharge were likewise simulated.March, and so on).

Validation by spatial samplings
There were seven field sampling regions in the Binggou watershed (A, B, D, E, F, H, I) (Fig. 1).Snow depths were sampled in these regions, and the statistical results were used to validate the calculated spatial SWE.Besides snow depth, snow density was also measured on different vertical snow layers.We computed the SWE by measuring the snow depth and the density, then the modeled and measured SWEs were compared.The sampling scheme, including each sampling region, was divided into several subregions with a size of 30 m, and the snow depth and the density was randomly sampled at each subregion.In this study, the average SWE of the total samples in each grid with a size of 60 m was counted to validate the simulated snow depth using a 50 m resolution.Thirty-two SWEs were collected and averaged during March 2008 (11 March, 14 March, 15 March, 17 March, 19 March, 22 March, 24 March, 29 March, 30 March).
The simulated SWEs agreed well with the measurements in most cases (Fig. 4).The root mean square error (RMSE) was 3.2 cm.Obvious errors occurred, mainly in the H and I regions.These overestimated SWEs are marked in the dotted ellipse in Fig. 4. Regions I and H are located at the lower altitudes (3450 and 3528 m, respectively) by the side of the channels in the Binggou watershed.Snowpacks just begin to melt in the middle of March.For higher air temperature and increasing solar radiation, snow at lower altitudes melts faster than in other places; hence, the SCA changes considerably in the H and I regions in a single day.However, there was only one MODIS image every day.In this case, the snowmelt may be overestimated because one unchanged snowcover image was used to calculate hourly snowmelt all day.
As described above, spatial and temporal SWE distribution can be retrieved effectively using MODIS data and the energy balance method.

Validation of snowmelt runoff by discharge measurements
By contrasting the measured and simulated discharges in the snow season, the total amounts and the temporal trend of calculated snowmelts were validated.Daily discharges in the snow season of 2008 were computed (Fig. 5).The simulated results agreed well with the measured discharge; the NSE was 0.673, and the volume difference Dv was 3.9 %.The refreezing phenomenon was not considered in computing the discharge.Abundant ice bodies, formed by refrozen snowmelt, were found in valleys and river channels in the Binggou watershed in the early spring of 2008.The thicknesses of these ice bodies varied from dozens of centimetres to a few metres.This phenomenon should have delayed the confluence of snowmelt runoff, the reason for which was quite possibly the inaccurate agreement of the simulated results with the measurement obtained during the earlier snowmelt season.

Discussion of other potential errors
There are some simplified parameterization schemes for spatial energy computation.These schemes consist mainly of simplifications of the terrain influence, ground heat changes and the spatial heterogeneity of precipitation.Terrain influence is an important factor in snow surface energy exchange in the alpine region.In this study, the effects of terrain on solar radiation were corrected by a simplified sky-view factor.By this method, the terrain features of solar radiation in an alpine watershed could be recognized.However, the influence of surrounding terrain shading on solar radiation was not considered in detail.Wind speed plays an important role in computing turbulent heat exchange.High-resolution wind field simulation is needed to predict snow transport and snowcover development over steep topography (Raderschall et al., 2008).Only a simplified interpolation method was used in this study.For more accurate wind-field simulation, another physical wind simulation model must be added.Soil heat flux density is usually regarded as a constant in modelling snow processes (Yang, 2008), a simplification that was also used in this study.In the snow season of 2008, the measured ground heat did not change very much in the Binggou watershed.Especially in the mid-winter, it was almost constant (approximately 10 W m −2 ).However, it fluctuated with increasing snowmelt during the later snowmelt season.When snow disappeared, the ground heat obviously increased.Although the magnitudes of ground heat are not large relative to other heat sources, such as shortwave radiation, a more accurate energy exchange computation must also be considered.Especially in the land-surface hydrology of cold regions, frozen soil plays an important role because the freezethaw cycle of soil significantly changes the soil hydraulic and thermal characteristics (Wang et al., 2010).The spatial heterogeneity of precipitation is relatively weak in a small watershed such as the Binggou; thus, only an elevation adjustment was used to correct the precipitation distribution.If a larger watershed was chosen to model the snow processes, the heterogeneity of precipitation should be considered in more detail.
Second, some elements in snow evolution processes were ignored in the computation, such as blowing snow sublimation and transport.Blowing snow sublimation plays an important role in snow hydrologic processes, and the redistribution of snowpack from interactions between blowing snow and terrain is prominent under high-wind speed conditions (Pomeroy and Li, 2000).Therefore, snow sublimation may be underestimated in the modelled results.However, because MODIS snowcover data were used to discriminate snowpack existence, even the redistributed snowpack can be recognized and included when calculating snowmelt.The errors from blowing snow transfer are relatively small in theory.

Conclusions
An energy balance method and MODIS data were used to model snow distribution and melt processes in an alpine watershed.As a result, the total SWE change over the entire snow season and the daily SWE distribution were obtained.Discharges were also simulated and then used to validate daily total snowmelt.Simulated snow-water equivalents were validated by field measurements.
To obtain continuous daily snowcover maps, a linear interpolation method was used.By this method, each MOD10A1 scene could be used and, thus, no data were wasted.The simulation results indicate that almost half of MODIS image waste could have been avoided in the Binggou watershed in the snow season of 2008.With further validation and improvement, this method can be used in large-scale snowcover mapping.By using the method introduced in this article, modelling spatial snow distribution and melt runoff and removing clouds effectively will become relatively convenient.
Using the methods described in this paper, the physical processes of snow ablation and the snowmelt distribution can be modelled accurately.SCA data from RS were used as an input for snowmelt modelling, and a reasonable combination of SCA data and the energy balance method was introduced in the paper.The results indicate that in the alpine region without abundant ground gauges, it is practical and relatively accurate to simulate spatial and temporal snow distribution using snowcover information from RS data and observations from several local weather stations.

Fig. 2 .
Fig. 2. Mean error between the interpolation and the observed snow cover fraction for increasing gap-length.

Fig. 3 .
Fig. 3. Different daily energy inputs to the snowpack of the total watershed in the 2008 snow season.

Fig. 4 .
Fig. 4. Validation of the simulated snow-water equivalent using measurements at sampling regions in the Binggou watershed (H(3/30) represents the measured snow depth in the H region on March, and so on).
• C in the 2008 snow season; the minimum temperature was −29.6 • C and the maximum was 19.9 • C.

Table 1 .
Adjusted parameters on the point scale.