Snow Cover Area Estimation Using Radar and Optical Satellite Information

Obtaining the seasonal variation of snow cover in areas of the Argentinian Andes is important for hydrological studies and can facilitate proper planning of water resources, with regard to irrigation, supply, flood attenuation and hydroelectricity. Remote sensors that work in the visible and infrared wavelength range are operational tools for monitoring the snow in clear skies. However, microwave satellites are able to obtain data regardless of atmospheric conditions. The advantage of using radar images is that they are very useful to obtain highly accurate parameters such as snow moisture depth, density and water equivalent resulting in improved forecasting models. In this paper, we analyze an ERS-2 image of the Andes mountain range in the northern region of the Neuquén province, Patagonia, Argentina. The objective was to obtain the spatial distribution of wet and dry snow and to compare these results with data from optical sensors (LANDSAT) in order to understand the topographic variables that influence the spatial distribution of wet snow. Optical information from sensors like LANDSAT TM 5 was analyzed to obtain fractional and binary snow indexes during a passage simultaneously with radar data. Surface temperature is used to study the association between the different types of snow altitudinal ranges and surface temperature. In this paper, we selected a scene on October 8th 2005. The entire methodology was systematized in a code implemented in IDL language.


Introduction
Snow is composed of ice crystals, liquid water and air.The ice crystals are deposited on the Earth's surface as a result of atmospheric precipitation.One of the characteristics that best describes a package of snow is the density ρs, which is typically found in a range of 0.2 to 0.6 mg•m −3 .Freshly fallen snow can have a density around 0.1 mg•m −3 .Density of older snow packages increases as a result of the compaction through the effects of wind and gravity or thermal metamorphism.
The internal structure of a snow pack is determined by the crystal grain size, generally defined as the mean radius or the equivalent radio of ice crystals.Usual grain sizes tend to be between 0.1 and 3 mm [1].The temperature of the snow is an important variable to determine whether snow is dry or wet.Dry snow does not contain liquid water.The presence of water occurs when the temperature is higher or close to 0˚C.With temperatures below zero the snow remains dry, whereas with temperatures greater than or equal to 0˚C will be wet.The factors that influence the snow reflectivity are the size of the grain, the depth and density of the snow layer and the amount of impurities it contains [2].Fresh snow presents more reflectivity than frozen snow, with the lowest reflectivity corresponding to dirty snow [2].
In the field of microwave remote sensing, Shi and Dozier [3], have described that the mayor contribution of the scattered signal from a wet snow package in C-B and is the volume backscattering, and the surface backscattering from the air-snow interface.These two kinds of backscattering are sensitive to snow properties such as liquid water content, density, shape and size of ice crystal and surface roughness.According to the mentioned authors, the liquid water content in snow caused a huge dielectric loss, which increased the absorption coefficient.Therefore, volume backscattering is proportionally inversed to snow wetness, and also proportionally inversed to density.
The potential of radar systems in retrieving wet snow has been investigated.Naglerand Rott [4] [5] have developed an algorithm for mapping wet snow in mountainous terrain using repeat pass Synthetic Aperture Radar (SAR) images from ERS SAR.Pettinato et al. [6], have developed a multi-temporal analysis of C-band SAR images of ERS SAR and ENVISAT ASAR combined with results of model simulations, to monitor the extent of snow cover in Italian Alps.An operational algorithm to produce snow cover maps from remote sensing data in the Italian Alps has been developed by Pettinato et al. [7].The algorithm can generate maps by combining optical data from MODIS and SAR data from ENVISAT/ASAR.
Salcedo [8] analyzed a sequence of ERS-2 images in the 2005-2006 period, and evaluated the influence of the snow melting process on an extraordinary flood of the Neuquén River during July of 2006.
The snow covered area (SCA) is an important parameter for understanding the seasonal variation in mountainous areas.This parameter is useful for hydrological studies, water resource management and climatology.Current approaches include methodologies and algorithms for the estimation of SCA based on the developments of Naglerand Rott [4] [5].This paper presents the estimation of wet and dry snow areas in the studied watershed using radar and surface temperature derived from LANDSAT TM 5.

Study Area Description and Data
The study area is comprised of the Neuquén river watershed and headwaters located in an area of the Andes mountain range and the Cordillera Del Viento range.It includes Varvarco river watershed, this zone is a complex terrain area with altitudes ranging between 2000 and 4500 m with steep slopes (see Figure 1).
Neuquén river basin has a hydrological regime which depends on rainfall (especially in winter season) and snow melt during springtime, and defined the behavior of flow along the hydrological year.Generally flashflood events take place during summer due to intense convective rainfall and caused several damages to people and buildings.However, during July 2006 a huge flood event occurred which almost caused the overtopping of the Cerros Colorados Dam.
Several meteorological, geomorphological and hydrological factors have played different rolls in order to explain the phenomena.Two year (2005 and 2006) of extraordinary snowfall, high temperature during winter and intense rainfall concentrated in July (above the mean value), have caused the flashflood event.In the framework of the hydrological situation explain above, this paper described a methodology to map snow cover area, then validate with surface temperature obtain from optical imagery, and contribute to clarify some details about the flood of 2006.
The methods discussed in this paper have been applied to ERS-2 images and altitude data obtained from Shuttle radar topography model digital elevation model (SRTM).Surface temperature values were obtained from theLandsat5 TM thermal channel.The period of data included the hydrological 2005-2006 cycle, which corresponds to a year characterized by extreme events of precipitation and snow accumulation.

Snow Cover Mapping Using Radar Data
The estimation of SCA and the snow classification as a function of moisture content was performedusingERS-2 images and elevation data obtained from a digital elevation model.
The estimation of the extension of the snow cover area and the classification of snow into wet, dry and bare soil classes were accomplished using radar data and the SRTM elevation model.The method involved the determination of a threshold classification, considering the Backscatter Coefficient (σ 0 ) versus height, to discriminate wet snow, dry snow and bare soil.Wet snow can easily be separated from dry snow or bare soil due to the high absorption of microwaves by the wet snow.
To reduce radiometric resolution errors, the backscattering coefficient can be calculated by simply using the average intensity value in the Precision Image (PRI) image of a distributed target using equation i.e. by averaging the intensity values of a certain amount of pixels within the group of pixels that corresponds to the target.
Then, the backscatter intensity is expressed as the average of the effective back-scattered section per unit of area referred to as backscattering coefficient (σ 0 ) (dimensionless coefficient): where Σσ t is sum of effective individual sections y A is the target area [9].
First, a band ratio was performed using two ERS-2 scenes: a reference image from the summer or winter season and the other dating from spring.This method was originally developed by Naglerand Rott [4], and is based on the changing behavior of snow over the seasons (Figure 2).During the spring season, when snow melting begins, the Backscattering Coefficient (σ 0 ) is low, whereas in winter or summer (with dry snow conditions or no snow cover respectively), σ 0 is high.Nagler & Rott [4] defined a threshold <−3 dB for wet snow.
Secondly, the Backscattering Coefficient values, obtained from the ratio image and the elevation were plotted into a 3D histogram (Figure 3).Three distinctive zones with particular pixel frequencies can be observed.The  Backscattering Coefficient (σ 0 ) and elevation thresholds are interactively selected in order to discriminate wet from dry snow and non snow covered surfaces.For the October 2005 ERS-2 image, the final threshold for the wet snow zone were <−3 dB and <2300 m.Areas without snow below 1500 m, showing wet snow around 2000 m and dry snow at higher altitudes.
Based on the conclusions of Nagler and Rott [4] and supported by the findings of Caves et al. [10] there could be find three situation in each radar images scene: • Situation A, where the reference backscattering coefficient ( ) is greater than spring case backscattering coefficient ( ) σ .This class is associated with presence of wet snow.
• Situation B where the spring case backscattering coefficient ( )

Surface Temperature
Backscatter coefficient thresholds are those proposed by Nagler and Rott [4], Caves et al. [10], Pettinato et al. [6], Pettina to et al. [7], however the representation of the relationship between σ 0 and height using isoliness implifies the definition of new thresholds between classes.In this sense, it is possible to search for specific thresholds within each image thus contributing to the improved discrimination of classes.
The Planck function is used to calculate the radiance emitted from blackbody objects (Liang et al., 2012): ( ) where L (W•m −2 ) is the spectral radiance, λ(m) is the wavelength, T (K) is the temperature of the object, h is Planck's constant (h = 6.62606896 ×10 −34 J•s), c is the speed of light (29979245.8m/s), e is the base of the natural logarithm and k is the Boltzmann Constant (1.380 6504 × 10 −23 J /K) [11].On the other hand, the radiance can be calculated from the digital levels, according to the following equation ( ) where L s is the spectral radiance in the thermal band (W/(m 2 sr µm)) [12], (obtained from the header file of images; band 6) and D n and D n max are the digital levels and maximum digital level of each scene and L max y L min are the upper and lower spectral radiance limit, with units of W/(m 2 sr µm)) Then, considering Equations ( 1) and ( 2), the brightness temperature of an object whose radiance has been measured by the sensor (Equations ( 1) and ( 2)) can be expressed as (Li et al., 2004): where T s is the effective temperature of the satellite brightness temperature (˚K) and K 1 and K 2 are constant (for Lands at 5 TM, K 1 = 607.76W/(m 2 mm•Sr) and K 2 = 1260.56K) [12] and L s radiance is integrated in the wavelength (Figure 5).Lands at satellites do not provide the surface temperature in an operational way due to the limitations of having a single band in the thermal spectrum to correct atmospheric effects and emissivity.To obtain surface temperature from satellite radiance, atmospheric profiles of temperature and water vapor are needed at the time of the passage of the satellite.In spite of the presence of weather stations in the study area, the information and data available are not simultaneous with the passage of the satellite.
Radiance emitted by an object on the ground is also dimmed and augmented by the atmosphere so this effect must be considered in the estimation of the ascending and descending transmittance and atmospheric radiance.In this study, we have used an operational tool available online for atmospheric correction, which allows you to calculate the atmospheric transmittance and radiance using the MODTRAN model 4.
The space-reaching radiance is converted to surface-leaving radiance by the radiative transfer budget equation [13]: ( ) where τ is the atmospheric transmission; ε is the emissivity of the surface; L T is the radiance of a blackbody target of kinetic temperature T (W/m 2 /Sr/μm), L u is upwelling radiance, L d is the down welling radiance (W/m 2 /Sr/μm) and L TOA is top of atmosphere radiance measured by the instrument.Barsi et al. [14] estimated that the surface temperature values can be obtained with an error of ±2 K in places where the emissivity is known and the atmosphere is relatively clear.
The method used in this paper assumed that the thermal emissivity was highly correlated positively with the normalized difference vegetation index (NDVI).
As a first approximation, it is possible to obtain the NDVI values from reflectivities obtained at the sensor, or at the top of the atmosphere (NDVI TOA ).However, it is more accurate to correct these values to obtain the surface reflectivity and estimate values representative of natural surfaces (NDVI S ), as the NDVI arises from a normalized difference, and small differences between NDVI TOA and NDVI S can be expected.The band 3 and 4 reflectances used in the calculation were atmospherically corrected using the Fast Line-of-sight Atmospheric Analysis of Spectral.
Hypercubes, (FLAASH) module of ENVI geographical information system.FLAASH is a first-principles atmospheric correction tool that corrects wavelengths in the visible through near-infrared and shortwave infrared regions.FLAASH incorporates the MODe rate resolution atmospheric TRANs mission radiation transfer code (MODTRAN4) in calculation.MODTRAN is a computer program designed to model atmospheric propagation of electromagnetic radiation for the 100 -50,000 cm −1 (0.2 to 100 um) spectral range.
Figure 5 presents the surface temperature map from October 2005, the borders where color turns to blue indicates the 1˚C isotherm (above 1600 m as shown in Figure 6).Coolest temperatures coincide with highest altitudes and warmest temperatures are located mainly inside the valleys (see Figure 1).

Discussion
The analysis of σ 0 versus elevation, allowed the identification of the best thresholds for estimation of snow classes and the differentiation of snow from bare soil.The thresholds favored mapping and final classification of areal snow cover.The whole procedure was systematized in an algorithm that automated the process.
Figure 6 presents the scatter diagram between temperature and height above the ground for image acquired on October 8th 2005.The LANDSAT 5TM image and the digital elevation model.Values of surface temperature below 277 K are presented in altitudinal levels above 1600 m, and surface temperature below 265 K only occur at levels above 3700 m.Surface temperature versus altitude, shows an inverse relationship with greater variability in the lowest levels.Figure 7 presents backscattering coefficient versus surface temperature, the maximum frequency of wet snow pixels are associated with surface temperature of 274 K (1˚C), which is around fusion temperature, and confirms the location of wet snow.As can be seen in Figure 5, Surface temperature reaching 274 K is located at levels above 1600 m.
The surface temperature isotherm of 274 K (1˚C) ± 1 K is associated with an average height of 2800 m above ground level (snow line proposal is the surface temperature isotherm corresponding to 274 K).Wet snow distribution as a function of the height demonstrated the greatest dispersion.The lowest observed level with wet snow presence at the considered case reached up to 2000 m above ground level.

Conclusions
The threshold methodology for the classification of snow, depending on its moisture content was implemented in ERS-2 images (band C).However it is hoped that in future researches it will be able to apply the method using L-band and X-band data, in particular taking into account the development of the SAOCOM argentine radar satellite (L-band).The versatility of this method lies in the possibility of searching for thresholds according to the characteristics of the study area and the conditions of the snow cover.
The maximum frequency of wet snow pixels corresponded to a surface temperature of 274 K (1˚C).These temperature values were associated with an average height of 2800 m above ground level.The snow line proposal was 274 K ± 1 K surface temperature isotherm.
Wet snow distribution as a function of the height demonstrated the greatest dispersion.The lowest level that wet snow was observed in our study area was 2000 m.
Specific thresholds for dry, wet snow and bare soil applied to each image, improve the discrimination of classes.

Figure 1 .
Figure 1.LANDSAT image of the study area, including the northern area of the Neuquén river watershed.
class could be associated with re-frozen processes.• Situation C where the backscattering coefficients has similar values in both images, this situation could be associated with bare soil (without snow) or dry snow.Taking into account the different situations explained above and with the study of the morphology of the isolines of the 3D histogram, allows differentiating the thresholds which clearly indicate the three coverages.The total extent of the snow-covered area (wet + dry snow) was compared with snow indexes obtained from LANDSAT images.Furthermore slope and orientation of the terrain was analyzed to detect wet and dry snow variation associated with topography.The snow cover area map for October 2005 is shown in Figure 4.

Figure 4 .
Figure 4. Snow cover map for October 2005 using LANDSAT 5 TM and ERS-2 radar images.

Figure 6 .
Figure 6.Scatter diagram of the soil temperature as a function of height on October 8th 2005.