Elevation based correction of snow coverage retrieved from satellite images to improve model calibration

Abstract. The most widely used method for snow dynamic simulation relies on temperature index approach, that makes snow melt and accumulation processes depend on air temperature related parameters. A recently used approach to calibrate these parameters is to compare model results with snow coverage retrieved from satellite images. In area with complex topography and heterogeneous land cover, snow coverage may be affected by the presence of shaded area or dense forest that make pixels to be falsely classified as uncovered. These circumstances may have, in turn, an influence on calibration of model parameters. In this paper we propose a simple procedure to correct snow coverage retrieved from satellite images. We show that using raw snow coverage to calibrate snow model may lead to parameter values out of the range accepted by literature, so that the timing of snow dynamics measured at two ground stations is not correctly simulated. Moreover, when the snow model is implemented into a continuous distributed hydrological model, we show that calibration against corrected snow coverage reduces the error in the simulation of river flow in an Alpine catchment.


Introduction
The snow dynamic and its spatial variability act as a "large mountain" reservoir due to the shift in time of the water volume determining the flow regime in the rivers network. Its role is well known to be very important on the river flow hydrograph for low as well as for high discharge regimes (Gurtz et al., 2002). For this reason and for its operative impact Correspondence to: C. Corbari (chiara.corbari@mail.polimi.it) on the management of water resources in mountain basins, a considerable amount of research effort has been done leading to several major improvements in the modelling technique (Verbunt et al., 2003;Huwald et al., 2006;Shamir and Georgakakos, 2006).
Despite the important role played by the snow and despite the capability reached by scientists and technicians to model its dynamic, one of the main uncertainties lies in the determination of the snow covered areas and of their water equivalent. This is due to the well known difficulties of monitoring snow parameters and of detecting the snow coverage on areas of complex topography or densely forested (Bocchiola and Rosso, 2007;Simpson et al., 1998;Baral and Gupta, 1997). Nowadays, the potentiality of images from remote sensing is well assessed by the scientific community (Baumgartner et al., 1994;Metsiimiiki et al., 1998;Parajka and Bloschl, 2008), but the conversion from qualitative data to quantitative ones is still very complex. This is due to the spatial and radiometric integration underlined in the pixel value of digital number given for any satellite image. This problem is evident in several fields of land observation.
Several methodologies have been considered to improve the use of satellite images ranging from the field of numerical modelling and ground measurements of snow parameters (Dressler et al., 2006;Bitner et al., 2000;Lee et al., 2005;Ranzi et al., 1999, Tachner et al., 1998 to the development of empirical relationships between flow discharge and snow accumulation and melt (Martinec, 1960).
In literature different techniques have been developed for operational snow cover monitoring (Romanov et al., 2003;Foppa et al., 2007). This work proposes a simple elevation based procedure to correct snow coverage retrieved from NOAA AVHRR (Advanced Very High Resolution Radiometer) satellite images (NOAA, 2000), from falsely uncovered pixels. Snow coverage is used to calibrate a distributed model, based on air temperature thresholds (Salandin et al., 2004), for the simulation of snow dynamic on alpine watersheds.
The proposed methodology is validated at local scale by the comparison of modelled and measured snow depth dynamic. However this approach gives only a qualitative verification of the timing of the snow dynamic because snow gauges do not provide information on the snow water equivalent. The methodology is also validated at basin scale where a good agreement between observed and simulated hydrographs at the river Toce basin outlet (Italy) is observed after model was calibrated against corrected images.

Site and data description
The subject area is the river Toce basin, a typical alpine basin with a total drainage area of about 1534 km 2 . The basin is located in the north of the Piedmont region in Italy (Fig. 1).
Climate conditions are typically humid, characterized by higher precipitations in autumn and spring. The annual average precipitation exceeds 2000 mm. Climatic characteristics, together with morphology and soil texture, frequently induced flood events in the past years. Most of the presented results are related to the Anza river basin, a tributary of the river Toce, with a total drainage area of about 261 km 2 .

Physiographic basin characterization
Available digital cartographic data include: the Digital Elevation Model (DEM) available in raster format at 100 m×100 m resolution; CORINE land cover maps (CEC, 1994;EEA, 2000) updated in the year 2000 available in vector format; pedologic characteristics for soils available in vector format. From the available basic thematic layers, basin parameters required for the application of the hydrological model have been derived at a spatial resolution of 500 m×500 m. These include: Curve Number (Soil Conservation Service, 1986), flow direction, slope and aspect, residual and saturated soil moisture, albedo, pore size distribution index, saturated hydraulic conductivity, wilting point, field capacity and soil depth.

Hydrologic and meteorological data
For this study, meteorological and hydrologic measured ground data were collected by the telemetric monitoring system of the Regione Piemonte. Rainfall, air temperature, incident short wave solar radiation and air relative humidity data are available from 1 January 2000 to 31 December 2004 at

Snow data: satellite images and ground measurements
The satellite images used to retrieve snow coverage are taken by AVHRR (Advanced Very High Resolution Radiometer) on board the polar satellites NOAA. The AVHRR (NOAA, 2000) is a six channel radiometer that detects energy in the visible and in the infrared spectral bands. These images are characterized by a nominal spatial resolution of 1.1 km at nadir while the resolution decreases to 3 km at the very edge of the scan. For the selected images, the study area is centred in the satellite scan to overcome this problem. Selected images are georeferenced using common software for remote sensing analysis.
The images are updated by a new passage of the satellite twice a day, one in daytime and the other after the sunset. However for this study only daytime images were considered.
For the four years studied, up to 23 images referring both to the accumulation period and, most frequently, to the melting period were collected. Since in the visible band reflectivity of cloud and snow is very similar (Henderson-Sellers, 1984), it is difficult to distinguish between them. Therefore, only images without clouds were used. Ground measurements of snow height are given by 1 snow gauge located in the river Anza basin at an elevation over 2000 m a.s.l. and by another snow gauge in the river Toce basin; these data are available with a time step of 30 min for the period 2000-2003.

Snow cover classification from satellite for complex topography
Snow covered pixels were classified with a supervised approach from the visible band (Choudhury, 1979). However some anomalies were encountered mostly on those pixels affected by shadow or densely forested coverage. In fact, it was noted that in some circumstances, pixels that should have been covered by snow were classified as snow free pixels. Inconsistencies in snow coverage retrieved from satellite can be detected using in situ information. Ground observations from snow station measurements ( Fig. 1) were compared against snow coverage of the corresponding pixel from the 23 satellite images used in this study. At Antrona station, on a total number of 19 days in which station measured a positive snow depth, pixel was classified as snow free 10 times, with an error of 52.63%. At Macugnaga station, on a total number of 9 days in which station measured a positive snow depth, pixel was classified as snow free 7 times, with an error of 77.7%. Influence of topography and forest on pixel classification is shown in Fig. 2. A subarea of the entire basin, characterized by pronounced peaks and deep valleys as well as forested coverage, is presented. Digital elevation model was resampled at the same spatial resolution of NOAA images (1100 m). A map of shadowed pixels was produced, related to the time when the satellite image was taken. A module taking into account the shadow effect induced by topography was developed. The algorithm calculates the angle, ψ between the point with maximum elevation in the direction of the solar beam, denoted by coordinates x m , y m , z m , and the examined cell, denoted by coordinated x 0 , y 0 , z 0 ( Fig. 3): If ψ is higher then the sun elevation, the cell is shadowed. An elevation based correction was applied to snow coverage retrieved from satellite image. According to this method all pixels above a reference elevation were considered as snow covered. Reference elevation is fixed as the mean elevation, taken from the digital elevation model, of pixels classified as snow covered.
The corrected image was compared to the original one and a strong correlation is shown between shadowed pixels and pixels falsely classified as not covered by snow (Fig. 2). In this area, forest seems to have less influence on snow cover classification as shown by forest cover retrieved from CORINE land cover map. The effectiveness of this simple method to correct snow coverage can be assessed taking into account ground observations from snow station measurements. Both at Antrona and Macugnaga station, the comparison with corrected snow coverage, reveals a complete agreement (100%). All anomalies have been removed by the correction procedure. In Table 1 the percentage of snow covered pixels from raw satellite images (second column) and after elevation correction (third column) is reported.

The hydrological model
Hydrological simulations were performed using the FEST-WB distributed water balance model (Montaldo et al., 2007;Rabuffetti et al., 2008;Ravazzani et al., 2008). FEST-WB computes the main processes of the hydrological cycle: evapotranspiration, infiltration, surface runoff, flow routing, subsurface flow and snow dynamics (Fig. 4). The computation domain is discretized with a mesh of regular square cells in each of which water fluxes are calculated at hourly time step.
The model needs spatially distributed meteorological forcings: precipitation, air temperature, air relative humidity, and solar net radiation, sum of shortwave and longwave components. The observed data at ground stations are interpolated to a regular grid using the inverse distance weighting technique. Spatial distribution of air temperature local measurements takes into account the reduction of temperature with altitude, with a constant lapse rate of −0.0065 • C m −1 . Thermal inversion phenomena are neglected. Shortwave net radiation involved in evaporative processes is calculated considering the effect of topography (Mancini et al., 2005). Longwave net radiation is evaluated as a function of air temperature (Goudriaan, 1977).
The snow module of FEST-WB includes snow melt and snow accumulation dynamics. The partitioning of total precipitation, P , in liquid, P l , and solid, P s , phase is a function of air temperature, T a (Tarboton et al., 1994): where α P is computed by:  where T low and T up are air temperatures below/above which all precipitation falls as snow/rain, respectively, to be found by means of calibration.
The snow melt simulation is based on the degree day concept (Martinec, 1960). The melt rate in m/s, M s , is proportional to the difference between air temperature and a predefined threshold temperature, T b : where C m (m • C −1 s −1 ) is an empirical coefficient depending on meteorological conditions and geographic location. The terrain covered by snow is supposed to be frozen and hence the melted water is prevented from infiltrating into the soil. Conversely, the liquid fraction of snow water equivalent, R s , sum of melted water and liquid precipitation, is supposed to flow cell by cell through the snow pack with a linear reservoir routing scheme (Ponce, 1989) with a celerity of 1.67 10 −3 m/s (Salandin et al., 2004). When R s reaches a cell not covered by snow, it is added to the liquid precipitation of that cell.
Runoff is computed for each elementary cell according to a modified SCS-CN method extended for continuous simulation  where the potential maximum retention, S, is updated cell by cell at the beginning of rainfall as a linear function of the degree of saturation, ε where S 1 is the maximum value of S when the soil is dry (AMC 1). The dynamic of the soil moisture θ, for a cell not covered by snow, is described by the water balance equation: where R is surface runoff flux, D is drainage flux, ET is evapotranspiration rate and Z is the soil depth. Soil moisture in cells covered by snow is assumed not varying with time.
The actual evapotranspiration, ET, is computed as a fraction of the potential rate tuned by a function that, in turn, depends on soil moisture content (Montaldo et al., 2003). Potential evapotranspiration is computed with a radiation-based equation (Priestley and Taylor, 1972).
The surface flow routing, computed on those cells not covered by snow, is based on the Muskingum-Cunge method in its non-linear form with the time variable celerity (Montaldo et al., 2007). Subsurface flow routing, similarly to the method implemented for the routing into the snow pack, is computed with a linear reservoir routing scheme (Ponce, 1989) with a celerity calculated as a function of the soil saturated conductivity.
The model was subjected to a rigorous process of calibration and validation. The available data sets for hydrological simulations were divided into two not overlapping periods.
The year 2000 was used as the calibration period. The remaining period was used for the verification of the model performance without any further adjustment of the parameters. Calibration of the infiltration related parameters was made by comparing simulated and observed discharge at the river Toce outlet. For more details the reader can refer to Rabuffetti et al., 2008. The calibration of the accumulation parameters T low and T up in Eq. (3) and of the snow melt parameters, T b and C m in Eq. (4), is described in the following subsections.

Calibration and validation of snow accumulation parameters from satellite images
The calibration of the snow accumulation temperature parameters, T low and T up , in Eq. (3) has been used in literature, from −1 • C to 3 • C (Tarboton et al., 1996), from −1 • C to 7 • C (Braun, 1991), from 0.5 • C to 1 • C (US Army Corps of Engineers, 1956) and Martinec and Rango (1986) referred to a maximum value of 5.5 • C.
The available images data set was divided into two not overlapping periods. The calibration period is from 1 August 2000 to 1 June 2002, while the validation period goes from 1 August 2002 to 31 December 2003. Calibration and validation were performed on the Anza river basin.
The calibration of T low and T up was based on "trial and error" approach, varying the value in the range −5 • C÷. 7 • C. To this purpose, simulated snow water equivalent maps are transformed in binary data: when snow depth is greater than zero, the pixel is marked as "snow" while, when height is zero, the pixel is marked as "no-snow". The calibration was based on two objective indices: the minimization of the root mean squared error (RMSE) calculated on the number of snow covered pixels and the maximization of the correct performance index (CPI) . Root mean squared error is defined as follows: where n is the total number of images, N p is the number of simulated snow covered pixels in a map, and N * p is the number of observed snow covered pixels in a satellite image. RMSE does not take into account the spatial distribution of the snow covered pixels. CPI overcomes this limit. In fact, computation of CPI is based on pixel to pixel analysis, where contingency tables (Mason and Graham, 1999) for comparing the observed satellite images and the simulated distributed data for different combination of T up and T low are built up. The contingency tables are a means to assess the ability of the model to correctly catch the spatial distribution of the snow or its temporal dynamic (Sect. 5  hit when there is a perfect agreement between the observed ground cover condition ("snow" -"no-snow") and the simulated one, the CPI is defined as the ratio between the number of hits and the total number of outcomes. We note (Fig. 5 and Table 2) that the values of the calibrating parameters that maximise CPI and minimize RMSE are −5 • C for T low and −2 • C for T up . However these values have to be discarded because they exceed the range accepted by literature.
The results of the calibration based on the corrected images are summarised in Fig. 6a and Table 2. The values of the best fit calibrating parameters are 0 • C for T low and 0 • C for T up , in the range accepted by literature.
In Fig. 6b results for the validation period are reported. As for calibration period the values of parameters that retain maximum value of CPI are 0 • C for T low and T up. .

Calibration and validation of snow melt parameters from satellite images
The calibration of the snow melt parameter C m in Eq.(4), is based on the comparison of simulated snow cover extent with the one retrieved from satellite images, as for the accumulation parameters, while the parameter T b is considered equal to 0 • C (Hock, 2003). In literature a large variability of the C m parameter values exists ranging from 2.7 mm d −1 • C −1 to 11.6 mm d −1 • C −1 , varying from site to site related to differences in the importance of the radiation components, air temperature and wind velocity (Hock, 2003;Martinec, 1960).
After fixing T low and T up equal to 0 • C, the calibration and the validation of C m were based on the "trial and error" approach. In Fig.7 the results are reported and the value of the calibrating parameter that maximise CPI is 4.32 mm d −1 • C −1 .

Simulated snow dynamic and ground measurements
To assess the influence that a suboptimal set of parameters calibrated against raw snow coverage may have on the simulation of snow dynamic, we compare simulation results to ground measured snow heights. We recall that the output of snow model is the snow water equivalent. The purpose of the model is, in fact, to take into account the effect of snow dynamics on the river flow at the basin scale. It is not its intention to reproduce the complex processes that describe evolution of the internal state of the snow cover (Essery et al., 1999;Diamond and Lowry, 1953) as it is required, for example, for avalanche forecasting (Brun et al., 1992). Therefore, snow density is not known and direct comparison of simulated snow water equivalent and measured snow heights is not possible. However, a possible approach is to compare the "snow"-"no snow" binary series to assess the time variability of the snow dynamics, namely the snow persistence, the beginning of accumulation and the end of depletion periods. A comprehensive index can be calculated from contingency tables that are build up including simulated snow data for different combination of T up and T low and observed snow depth, similarly to method explained in 4.1. The snow melt parameters, C m and T b , are considered constant and respectively equal to 4.32 mm d −1 • C −1 and 0 • C in this analysis. In Fig. 8 contingency tables for five combinations of T low and T up are reported for the studied stations. For the Macugnaga station the combination of temperatures that better represents ground cover condition is T low = T up = 0 • C (CPI = 86.2%), with only 202 discrepant days on a total of 1458, confirming the results obtained from calibration of distributed model. The other combinations of temperatures, T up = 3 • C and T low = −1 • C and T up = -2 • C and T low = −5 • C, denote CPI of, respectively, 83.1% and 85.4%. For T up = 0 • C and T low = −3 • C CPI is equal to 83.7%, while for T up = 7 • C and T low = 1 • C CPI is 78.5%. Also for the Antrona Lago station, T low = T up = 0 • C, is the combination of parameters with the greater value of CPI (89%).

Distributed model simulation results
The effect of snow model calibration on discharge simulation are shown in Fig. 9. The graphs refer to a flood event period occurred in November 2002 at the river Anza outlet. In these graphs five series of discharge simulated by means of the FEST-WB model with different combination of temperature thresholds for snow accumulation, are reported. The snow melt parameters, C m and T b , are considered constant and respectively equal to 4.32 mm d −1• C −1 and 0 • C.
Data are extracted from the continuous simulation for the period 2000-2003. We see that the hydrological model with T low = −5 • C and T up = −2 • C simulated the highest volume and peak discharge. On the contrary, the hydrological simulation with T low = 1 • C and T up = 7 • induced great snow accumulation and, as a result, the lowest peak discharge and cumulated volume.
In Fig. 10 three flood events at Candoglia are shown. Two autumn events (October 2000 andNovember 2002) and a summer event (June 2002) are considered and the observed data are compared with the simulated discharges with T up = −2 • C and T low = −5 • C and with T low = T up = 0 • C.
After the calibration with the uncorrected satellite images, FEST-WB seems to overestimate peak discharges. Model improvement is evident after the calibration with the corrected satellite images (Table 3) in terms of flood volume and peak discharge errors and of the Nash and Sutcliffe (1970) model efficiency defined as follows: where n t is the total number of time steps, Q o,j and Q m,j are the observed and modeled discharges at time step, j , respectively, and Q o is the mean of the observed discharges. The Nash and Sutcliffe efficiency is commonly used to assess the predictive power of hydrological models. It can range from −∞ to 1. Essentially, the closer the model efficiency is to 1, the more accurate the model is. Table 3. Numeral criteria for the intercomparison of observed hydrographs and FEST-WB simulation with T up = −2 • C and T low = −5 • C and with T low = T up = 0 • C for three flood events at Candoglia. η indicates the Nash and Sutcliffe (1970) model efficiency, ε q indicates the peak discharge percentage error and ε V indicates the cumulated volume percentage error.

Conclusions
The paper presents a simple operative procedure to correct NOAA-AVHRR satellite images respect to shadowed pixels induced from high crest in mountain areas and forest cover. The correction consists in assigning snow coverage to those pixels with an elevation greater than a reference value, calculated from observed map as the average snow coverage elevation.
Snow coverage maps retrieved from satellite images were used for the calibration of a snow model, based on air temperature, integrated in a distributed hydrological model. The calibration procedure that uses raw classified images, however, led to the determination of temperature thresholds in a range not accepted by literature. Instead, when using corrected snow coverage maps, the best result was most frequently achieved when both two temperature parameters were set to 0 • C. Moreover, simulation of discharge hydrograph at main river section reached best performance when the temperature thresholds calibrated with the use of corrected snow extent were used.
The snow model and the methodology for its calibration were validated at local scale, too. Comparison between simulated snow water equivalent and observed snow depth at two gauges showed that temporal dynamic of snow was better simulated when the two temperature parameters were set to 0 • C.
The presented procedure is suitable also for application to those river basins where snow depth measurements are not available. is also supported by Italian Ministry of University and Scientific Research (2006), project "Assimilation of remote sensing and ground data for the calibration of distributed hydrologic models and flash flood forecasting". The authors thank the ARPA Regione Piemonte, Italy, ARPA Regione Lombardia, Italy, Ufficio Federale dell'Ambiente UFAM Berna, Switzerland, Ufficio Federale di Meteorologia e Climatologia MeteoSvizzera, Switzerland for providing the data used in the case studied. The contribution of the reviewers to the improvement of this paper is gratefully acknowledged.
Edited by: A. Gelfan