Snow cover analysis in Emilia-Romagna

In this paper we present an operational chain developed in the Emilia-Romagna region (Italy) to monitor snow cover and snow water equivalent (SWE) over the area managed by the Regional Catchment Technical Service. Remote sensing data from medium resolution sensors (MODIS and AVHRR/3) are used as input for snow cover detection algorithms. Hourly weather station measurements are used as input for a snow melt and accumulation model in order to estimate the snow water equivalent. The model reliability is mainly related to the network density of heated rain gauges detecting snow precipitation. The products are disseminated as web bulletins and they are used in monitoring chain and alert for the Civil Protection Agency.


Introduction
The present research demonstrates how to detect snow spatial extent and to estimates snow water equivalent (SWE) starting from MODIS half kilometer L1B images (MOD02HKM) and meteorological data acquired by automatic regional weather stations. MODIS data are processed using band ratios and a decision trees approach in ENVI ITT-VIS environment, while meteorological data are processed using the software PRAGA [Antolini and Tomei, 2006]. In case of lack of availability of MODIS L1B images, the procedure can be applied on AVHRR/3 data from the METOP constellation at 1.1 km spatial resolution. Also METEOSAT Second Generation (MSG) images can be classified for snow cover purpose. The ground resolution of 3 km is not optimal because the grid is too coarse, but the real time access to data and the frequency of acquisition (every 15 minutes instead of one per day) may be a good compromise in case of cloud cover present in the images acquired by other sensors. In fact it is easier to get a cloud free scene during the daylight. A typical example is the presence of fog that may be intense during the first hours of the morning and then decrease in the afternoon. Also RGB composite can be used to classify snow cover. The 721 composite distributed in the MODIS Rapid Response System can be quite interesting because they are published few hours after the satellite pass. Snow has a typically low spectral reflectance in the band of short-wave infrared at 1.6 μm [Hall et al., 1995]. In Figure 1 we show an example of snow's spectral signatures of snow collected from MODIS data acquired February, 2th 2006. The signature varies both in relation to exposure and to the state of the snow that can be: compacted, wet or fresh. Snow detection procedures from remote sensing images were first performed with Landsat Thematic Mapper data and were then developed for MODIS sensor [Hall et al., 1995]. The National Snow and Ice Data Center (NSIDC) are in charge of the distribution of MODIS snow cover dataset. Different products are available via on-line services and they are well described in NSIDC and by Hall et al. [2002]. Real time acquisition of snow coverage represent an important factor in snow analysis, especially when the data are used for landslides risk management, therefore the possibility of acquiring image data in a few hours after the snow sensing is crucial. For this reason in Emilia-Romagna we opted for the Level 1B data instead of an higher level of product as provided by NSIDC [Spisni et al., 2007]. Different models can estimate the SWE. Usually models need several kinds of inputs, such as meteorological data used to compute the output variables of the model. Snow data can also be used to estimate some parameters. Hydrological data are necessary in the validation stage and morphological input are needed for basin elevation zoning. Usually remote sensing data are used only for snow cover detection and model corroboration. Snowmelt models can generally be lumped into four categories: 1) finite-element based mass and energy balance models [Lehning et al., 1999]; 2) temperature index models [Maidment, 1993]; 3) hybrid temperature index models [Federer, 1995]; 4) "simple" mass and energy balance models [Tarboton and Luce, 1996].
An interesting model of the first category was developed by the Institute for Snow and Avalanche Research (SLF) and it is called SNOWPACK [Lehning et al., 1999, Lehning et al., 2002. Worldwide studies show that this model can be applied to different environments, such as: alpine, arctic, maritime and continental snow covered area. The model solves the mass and energy balance equations using a finite element numerical scheme. However, the drawback of this model is that it requires many input variables, regarding meteorological, snow, soil and vegetation data. To run the model it is therefore necessary to set up a complex monitoring network, which is not available in the Emilia-Romagna region.
The simple mass and energy balance models fit the 'ideal' models with less input variables. Unfortunately, these models still include a level of complexity that requires detailed computer modeling and long run times. The advantage of these models is that they physically represent the mass and energy transfer in a snow pack without requiring detailed, multi-layered, and in some cases, iterative numeric solutions. Typically, the level of complexity inherent in these models is determined by how they simulate the surface and snow pack temperatures. The model used in this work is the snow accumulation and melt (SAM) model [Brooks et al., 2007], which represents a simple approach to the SWE determination, similar to the point based Utah Energy Balance (UEB) model [Tarboton and Luce, 1996]. UEB model predicts changes in the snow water equivalent depth by performing iterative mass and energy balance calculations on a single layer that includes the snow pack and a thin (0.2 -0.4 m) soil layer. In the UEB model the snow surface temperature is obtained by balancing the energy flux at the snow surface with the pack heat conductance using an iterative solution. Anderson [1968] showed that it was possible to replicate snow surface temperature without an iterative solution using approximations. The SAM model, similar to the approach of Anderson, does not require any iterative solutions and predicts snow surface temperature using a simple damping depth approximation. Moreover, SAM does not need neither detailed hydrological input nor snow data from field campaign. The model is based on hourly meteorological data only, that are easily automatically collected. The output of the whole processing chain is used for the snow bulletin of the Emilia-Romagna region (http://www.arpa.emr.it/sim/?telerilevamento/innevamento) and for landslides and floods risk for Civil Protection's alerts. Figure 2 shows the workflow.

Study area
The study area is located in the northern part of Italy and its extent is about 25000 square km. The area of interest is defined as the whole river basins managed by the Regione Emilia-Romagna. The administrative boundary cover the whole Emilia-Romagna and small part of: Toscana, Marche, Veneto, Lombardia e Ligura (Fig. 3).
The area has two main rivers that flow to the Adriatic sea in the central/northern part of Emilia-Romagna coast: the Po and the Reno rivers. Then there is a group of small rivers located in the eastern part, which are called Bacini Romagnoli. They flow in the Adriatic sea in the central/southern part of Emilia-Romagna coast. The area has about 10000 km 2 of plain land and the remaining land is located in the Apennines mountain chain, which crosses from north to south, the whole Italian country. At low altitude the mountain zone has widely diffuse badlands and agricultural activities. At increasing altitudes, deciduous forest, such as oak and chestnut forest, and pasture prevail and they become the main land cover / land use.
In The weather stations used in this study belong to the ARPA-SIMC (Regional Environmental Protection Agency -Hydro-Meteo-Climate Service) real-time network, which at present consists totally of almost 500 stations. Hourly measurements of temperature, precipitation, relative humidity, wind speed and direction and short wave radiation are recorded by the stations. The data are automatically logged and transmitted to ARPA-SIMC via three systems: GPRS, general radio packet service and modem. Then they are stored in an ORACLE database and then controlled and validated every day by both an automatic and a human-supervised quality control.  [Wen, 2008] in case of non centered acquisition and georeferenced to UTM 32 WGS 84 projection coordinate system. All the procedures are done using MODIS Conversion Toolkit. The image is then resized over the study area.

Snow detection
The TOA reflectance image is then processed in two steps: first we calculate the normalized difference snow index (NDSI) [(b4-b6)/(b4+b6)], then a decision tree is run based on NDSI, bands 2 and 4. The NDSI is particularly sensitive to the presence of snow, because it exploits the low reflectance values of the band 6 compared to the band 4. The thresholds reported in literature have been calibrated over areas of higher snow intensity and coverage than our region [Riggs and Hall., 2004]. Therefore, thresholds have been corrected to better fit Emilia-Romagna types of snow. In Table 1 we report the parameters applied to the decision tree. The snow cover map may present an overestimation of snow localized along the shoreline or wetlands. Sometimes it can also be affected by multilayer clouds characterized by over freezing. In these cases, the classification is correct on screen to eliminate false snow cover pixels. All procedures are performed using the software ENVI (ITT-VIS).
In cases where MODIS L1B is not available, the classification depends on the quality of the other data source available, as described above, and may be vary from unsupervised 5 classes (land, water, cloud, snow and mixed pixels) K-Means algorithm, supervised Spectral Angle Mapper with the same classes of unsupervised classifier or object oriented analysis based on ENVI spatial feature extraction module (SFE / ITT-VIS). Object oriented operator is used only to generate the mask or vector layer containing snow cover detected due to the high contrast in RGB composite. Depending from image pixel size, the scale level applied is usually greater than 80 and merge level between 80-90. We use also the advance option of thresholds to better detect only the snow cover pixels. Typically we search for MODIS 721 composite (250 m or 500 m spatial resolution) in MODIS Rapid Response System (AERONET_Ispra Subset) downloadable after few hours from acquisition. The RGB composite shows the snow cover in cyan (because it reflects light mainly in the green and blue channel of the RGB composition) and using a feature extraction algorithm is quite easy to classify in comparison with other land cover. Some kind of clouds may present similar spectral response. If these images are cloud covered or not available, than we check AVHRR/3 Full Resolution Area Coverage (FRAC) L1B data from METOP at 1.1 km or MSG data at 3 km spatial resolution. Full disc MSG image is georeferenced using an ENVI (ITT-VIS) tool downloadable from the on line user contribution code library, then the classification is performed. MSG data are also available through our receiving station (http://www.arpa.emr.it/sim/?osservazioni_e_dati/satellite). The four classification schemes (MOD02, MODIS AERONET, AVHRR/3 and MSG) can guarantee to extract the snow cover extension in near real time, exploiting different remote sensing data distribution internet sites.

Model validation
The University of Bologna conducted a SAM model validation, comparing the snow distribution simulated by the model with the snow cover estimated from remote sensing images in about 90 days during the period 2001-2007. The test considered only the days without cloud cover. The comparison shows a good correlation between the two methods (R 2 = 0.74) and it is statistically significant (P <0.05), only 14.4% of the data show differences greater than 10%. Figure 4 represents the relationship between snow cover from satellite and from model in the 90 days analyzed. It is a linear regression with prediction interval set at 95%. In more detail, the model is not always exhaustive but in several instances because of lack of precise experimental input data. In some critical events (around 8%), the relationship between remote sensing and model data is out of the prediction interval. For example in December 2001 and March 2005 the model underestimates snow cover probably due to wrong precipitation data measured by non-heated rain gauges. On the contrary in November 2005 and January 2006 the model overestimates the snow cover, this is probably due to the overestimation of model simulation in the plain and to undetectable data in satellite observation. These are pixels with snow cover but not classified by algorithms, due to the low level of snow cover not detectable by the satellite sensor. It was found, as expected, that the maps produced are reliable mainly in mountain areas where the density of heated rain gauges is adequate. In plain snow cover estimated by remote sensing is more accurate than data estimated by the model, due to a lack of heated rain gauges. The data measured by non-heated rain gauges (prevalent in the study area) false the simulation: the non-heated rain gauges does not record precipitation during snowfall and on the contrary they record precipitation during the hottest hours of the following days, because the instruments measure the snow melt. One example of this anomaly is shown in Figure 5. It reports the data collected of temperature (°C) and precipitation (mm) by two weather stations: Sestola (985 m asl) with heated rain gauge and Sassostorno (971 m asl) with a non-heated rain gauge, in the period from 1 to 10 January 2009. The two stations are located with a distance of 6 km. The Figure 5 shows how in Sestola the precipitation events were recorded with temperature below 0°C, while in Sassostorno there is a shift of the precipitation events when temperature increase above 0°C. In January 2009, the Emilia-Romagna region was covered by 41 heated rain gauges working and 13 damaged. During the winter season heated rain gauges must be regularly checked in order to control that they are correctly working. This check is done by the user once a week, by means of PRAGA software. The test involves checking the coincidence of known snow events and data recorded by instruments. In case of shifted plots, the rain gauge is turned off from the model and its failure is reported to the maintenance staff. The best agreement between model and remote sensing data is observed in the pixels with higher values of SWE, whereas the pixels where SAM model estimates low SWE values (< 10 mm) are not much reliable. This was confirmed during a survey conducted in January, 29 th 2009, along the Lavino, Samoggia and Panaro valleys, which lay between 400 to 800 m asl. Figure 6 shows the comparison between field data and MODIS classification that confirms remote sensing snow detection was accurate, even with a shallow snow layer. The survey points recorded as no snow but classified as snow covered from satellite were partially covered by snow. For the same day SAM model overestimated the area covered by snow in the pixels with low SWE values. Since low SWE output values of the model are not much reliable, it could be useful to adjust them with remote sensing analysis, when this is reliable (no cloud cover days). In more detail, we impose a minimum threshold (i.e. 5 mm) for the definition of snow presence in the SAM model output. This causes a general underestimation of modeled snow cover with respect to the satellite data, that is balanced adding all snow covered pixels estimated by remote sensing and not by the model, setting them to the threshold value. Therefore, the value of SWE in these points has to be considered fairly accurate. The integration between satellite and model also allows increasing the data dissemination. Typically, bulletins are produced only on days with low cloud cover, whereas for the Civil Protection activities, the model is run daily. In Figure 7 is shown how the integration of of satellite and model data can simulate snow under clouds. The snow surfaces detected is then analyzed to extract the data based on topographic zoning at provincial level (Tab. 2). The zones were defined from SRTM [Void-filled seamless SRTM data V2, 2006] at 90 m spatial resolution. The procedure is performed using an IDL script (ITT-VIS). At catchments level estimations of snow cover, SWE (mm), volume of water estimated (m 3 ), mean SWE over the snow-covered area and over the whole area are produced. This 90 m SRTM, much more detailed than the 500 m, could be used in future studies for specific areas having a high potential for landslides.

Snow and landslide: application's example
Snow's melting is a major cause for trigger of landslides and in certain conditions can contribute to fluvial floods. The melting of snow brings an increase in the amount of soil moisture and saturation, with consequent rise of groundwater and increase the degree of saturation of the soil. The relationship between landslides and snow melt is illustrated by the monthly distribution of major landslides. Figure 8 shows the monthly distribution of the occurrence of landslides in the historical archive of the Emilia Romagna with the trend of average monthly precipitation.

Figure 8 -Annual landslide distribution in Emilia Romagna (Servizio Geologico Sismico e dei Suoli dell'Emilia Romagna) and mean monthly rainfall.
The monthly variations in the frequency of landslides are quite consistent with the monthly variations in the amount of precipitation. Monthly distribution of rainfall has two maxima (autumn and spring) and two minima (summer and winter) that follow the distribution of landslides. It should be noted however that even if the autumn rains begin in September and are relevant in the month of October, most of the instability is concentrated in the month of November. This phenomenon can be explained by the fact that rainwater takes a long time to soak dried clay soils and fractured brought by the summer months. This confirms how important is the degree of saturation of the soil at the time that meteorological event triggers landslides. The second maximum of landslide is identified in the spring. The frequency of landslides, that occur in March instead of April (which is the wettest month), can be explained considering that in March snow melt takes place. In addition to the thickness of the snowpack, for the trigger of landslides and the formation of fluvial floods it is very important the time when the melting occurs. Indeed, a rapid melting of snowpack, combined with a sudden increase in temperature is equivalent to an heavy rain that, concurrent to rainfall, increases the amount of meteoric inflows in river basins by generating a rapid rise in water levels along the river courses and triggering landslides . Instead a slow and gradual melting of the snowpack helps to feed ground water, maintaining high water levels in rivers and high degree of saturation of the soil. SWE data therefore can be useful for improving the evaluation of the initial state of soil moisture for forecasting the level of hydrological and hydraulic risk estimated according to the DPCM's law (Direttiva Presidente Consiglio dei Ministri) signed in February 27, 2004 and applied by the Civil Protection at regional level.
The SWE allows to estimate the presence of snow and by comparing two maps of SWE, prepared in two different time intervals, it is possible to quantify the snow melt. These data are important for assessing the initial conditions of the soil at the beginning of a meteorological event that may generate fluvial floods and landslides. Thus it is known that rain can cause different effects depending on the degree of saturation of the basin and the presence of snow. These kind of maps are published in web pages accessible by invitation only, called INFOMET. INFOMET is a web application developed by ARPA-SIMC able to visualize meteo data time series images. It is used by ARPA-SIMC, regional and national Civil Protection weather prediction personnel to visualized Numerical Weather Prediction (NWP) and other weather and hydro-geological data. The role of remote sensing data is fundamental because it can correct errors in snow distribution with a more realistic detection of presence/absence of the snow. During winters with a high number of snowfalls, SAM model can make mistakes due to instrumental errors, for example because of an heated rain gauge broken, or a few number of meteorological sensors useful for snow monitoring and model errors. As the temperature rises and the snow is replaced by rain, the SAM model may overestimate/ underestimate snow distribution due to underestimation/overestimation of snow melt.
Since we can manage more satellite sensors, we can have good chances to acquire a cloud free scene. During the season 2009/10, the integration of MGS images was very helpful because it was possible to have an image every 15 minute. In this case the temporal aspect was more prominent than the spatial resolution. An example of the importance of the estimate of SWE and snow melt is provided by the weather event of January, 20-21, 2009 which affected the Apennine basins in the centralwestern region. The period was characterized by the combination of heavy rainfall, an increase in temperature begun January, 18 th and a consequent snow melt that started flood events at the lower altitude. The comparison between the maps of SWE of January, 17 th at 23:00 and January, 20 th at 23:00 shows the contribution of snow melt in terms of inflows ( Fig. 9), especially in the catchments of Trebbia, Taro, Parma, Enza and Secchia. From a quantitative point of view, the snowpack reduction can be estimated at about 30-80 mm of water equivalent, averaged over a large part of the Apennines, with a peak of 80 mm in the Apennines between Piacenza and Reggio Emilia (Fig. 10). The amount of water from snow melt and from rain caused a large increase in inflows in the basins of rivers Trebbia, Taro, Parma, Enza and Secchia with consequent formation of floods and high volumes of runoff. This meteorological event have triggered numerous small landslides, involving the superficial layer of soil with deterioration and collapse. The highest numbers of disruptions were recorded in the basins Trebbia, Taro, Enza and Parma where, in addition to precipitation, there was a higher melting snowpack which was added to the precipitation contribution. During winter 2008/2009 these products were also used to analyzed which administrative municipalities where more affected by snow cover and to define contributions to snow's street cleaning.

Conclusion
Emilia-Romagna, unlike Alpine regions, has an inadequate snow monitoring network. The methodology here presented allows to monitor the snow cover and to estimate the SWE. It is based on data automatically collected by satellite and by weather stations and it needs particular expertise in snow detection and real time hydrological state of study area.
The tests carried out confirm that the integration between model and remote sensing analysis gives a good estimation of snow cover as compared to field measurement. Satellite data can correct false classifications of snow presence estimated by the model, especially in plain where no heated rain gauges are present. On the other hand, the model may add information in unclassifiable areas, due to partially cloud cover images. Model limits are related to the inadequate coverage of heated rain gauges in non mountain areas. Moreover heated rain gauges need to be monitored during the winter season in order to check that they are correctly working. Civil Protection take advantage of this product because it receives information about the snow cover and the volume of water that can be released. With these data and the weather forecasts, it is possible to monitor the hydrological risk and to alert finally the population in case of floods and/or landslides hazards.