Interactive comment on “ Snow density retrieval using SAR data : algorithm validation and applications in part of North Western Himalaya ”

General comments The subject of the paper (i.e. the retrieval of snow density) is quite interesting and is one of the open issues on the table of the scientific community. Also, it is very interesting to look at the performances of the retrieval approaches used in the work since, up to now, the methods have been used in very few works. In particular: the tested approaches make use of backscattering data at frequencies which are poorly affected by the presence of dry snow. In fact, a common assumption is that dry snow can be considered transparent for frequencies below C-band (this is one of the reason because the proposed ESA mission CoReH2O was designed at X and Ku-band). This because the volume scattering of dry snow at Cand L-band is negligible. Overall the work seems to have been done in a good way, however, there are at least a couple


Introduction
The Himalaya holds one of the largest reservoirs of fresh water in form of glaciers and snow outside the Polar region.About 10 % of total area of the Himalaya is covered with glaciers and additional area, nearly 30 % support the snow cover (Singh and Singh, 2001;Singh et al., 2011).According to estimation, there are about 9575 glaciers, covering an area of about 38 000 km 2 in the Indian part of the Himalaya (Raina, 2009).The snow physical parameters, such as snow wetness, which shows the degree of liquid water content in snow pack, along with snow density and Snow Water Equivalent (SWE) are most important parameters for many water resources related studies such as snowmelt runoff and snow avalanche modeling (Rees, 2006).The traditional survey of these parameters is very expensive and difficult in rugged Himalayan mountains.The optical remote sensing has been used effectively to map the snow cover area (Rango and Salomonson, 1975;Dhanju, 1983;Dozier, 1984Dozier, , 1989;;Hall et al., 1995), qualitative snow wetness (Gupta et al., 2005) and snow grain size with hyperspectral data (Dozier, 1998;Dozier and Painter, 2004).The mountain shadow, clouds and retrieval of physical properties of snow using optical data remains a major gap area, and in these cases, SAR offers a better alternative to estimate snow and glacier dynamics parameters, as shown in various studies (Braun and Rao, 2001;Rau et al., 2001;Storvold et al., 2006;Thakur et al., 2009;Venkataraman and Singh, 2009), and is discussed in detail in next section.The current research paper is divided in 4 sections.In the first section introduction, basics concept of backscattering from dry snowpack and literature review on SAR for snow density retrieval is given.Study area details, data and software used are given in Sects.1.3 and 1.4.The main methodology used is given in Sect. 2. The SAR data pre-processing, ground truth and details of snow density retrieval algorithm are given in Sects.2.1, 2.2 and 2.3.Results and discussion are given in Sect.3, with detailed discussion on results of snow density inversion model, validation and accuracy of retrieved snow density, variations of snow density with elevation and aspect zones, limitations of the present study and its possible applications are given in Sects.3.1 to 3.5.Conclusions are presented in Sect. 4.

Backscattering from dry snowpack
Radar back scattering from a snow surface depends upon the dielectric constant of the surface (Tiuri et al., 1984;Ulaby and Stiles, 1980;Rott, 1984), its roughness properties, and the geometry of the scattering (Evans, 1963;Ambach, 1980;Nyfors, 1982;Hallikainen et al., 1986;Matzler, 1987).Many mathematical models of surface scattering have been developed, some based on physical laws, some on empirical data fitting, and some on a combination of the two (Rees, 2006;Strozzi, 1996;Strozzi and Matzler, 1998).Backscattering from a typical snow pack can be appropriately modeled by a sum of contributions from the snow surface, the volume scattering from the snow layer and the contribution from the ground (Fig. 1).The total backscatter for the wet 1929 1.2 SAR for snow density retrieval Shi and Dozier (2000a, b) demonstrated the use the SIR C-L-X datasets of Mammoth Mountains, CA, USA for snow water equivalence estimation by inferring snow density, subsurface properties and snow depth.Tadono et al. (2002) developed an algorithm for estimating snow hydrological parameters in wet snow regions using combined C and L-band satellite based SAR data for a part of the Hokuriku District of Japan, well known for having wet snow from the beginning of the winter season.They developed an algorithm to estimate snow density and roughness distributions using Radarsat images acquired in February.Furthermore, this study showed the possibility of inferring snow depths from combined parameters estimated with an L-band SAR image.Snehmani et al. (2010) used C-band SAR data of ASAR-APS dual polarisation for snow density estimation.In this approach for developing an algorithm for snow density estimation, the volume scattering model and the small perturbation model have been used with an exponential correlation (Ulaby et al., 1986) function for the surface backscattering contribution from the snow-ground interface.The final inversion model is a function of only dielectric constant of the snow.Niang et al. (2007) used new inversion method for snow density and snow liquid water content retrieval using C-band data from ENVISAT/ASAR alternating polarization in alpine environment.
In context of Indian Western Himalayans, small research has been done using polarimetric SAR data to infer snow density, except for the recent work reported using dual polarization ENVISAT-ASAR and Advanced Land Observing Satellite (ALOS)-Phased Array type L-band Synthetic Aperture Radar (PALSAR) data (Singh and Venkataraman, 2007, 2009, 2010;Snehmani et al., 2010;Thakur et al., 2008Thakur et al., , 2012)).None of these studies had used fully polarimetric data for inferring the snow parameters; therefore, 1930 the present study is an attempt to address this gap area as well as to evaluate the existing methods of snow density estimation using C and L-band SAR data in the part of North Western Himalaya (NWH).

Study area
The Beas River Basin upto Manali town with area of 350.21 km 2 has been selected for this study (Fig. 2).• C in June (Thakur et al., 2012).

Data and software used
The data used in this study are listed in

Methodology
The overall flow chart of SAR pre-processing, snow density retrieval and validation is shown in Fig. 3a.The details of SAR pre-processing, snow density retrieval and validation are given in subsequent sub-sections.

SAR data preprocessing
The SAR preprocessing was done with SARSCAPE software for converting raw SAR images into geo-coded backscatter images.The raw ASAR images (*.n1) and RS2 quad polarization (*.xml) were first read and converted into single look complex (SLC) images.These SLC images were multi looked to create the power (*.pwr) images.In mountainous areas, the topography significantly affects the geometry and radiometry in the SAR image, thus geo-coding and calibration using Digital Elevation Models (DEM) are required in order to relate the data to signature data (Frei et al., 1993;Holecz et al., 1993;Rosich and Meadows, 2004).The radiometric and geometric calibration of these multi-looked images was done using ASTER Global Digital Elevation Model (GDEM).The local incidence angle map (Fig. 3b), layover and shadow maps generated during calibration step have also been used in radiometric calibration and masking of layover/shadow areas.All ASAR backscatter images were re-sampled to pixel size of 25 m, RS2 images to 10 m and ALOS-PALSAR images to 15 m, with Universal Transverse Mercator (UTM) projection, zone 43, datum World Geodetic System (WGS) 84.
In this study, four masks have been used, (a) mask for local incidence angle, (b) mask for layover and shadow (Fig. 3c), (c) mask for forest (Fig. 3d), and (d) mask for no snow area using DEM and MODIS/LISS-III SCA.Out of these, forest mask (∼ 15 % of total area, which has been derived from Landsat ETM data (Table 1), is common to all the datasets, and other masks varies as per SAR image characteristics and time of imaging.These masks are needed, as the present inversion models are not formulated to account for all incidence angles, layover/shadow areas and compensating backscatter from forest areas.

Ground truth
The extensive ground truth campaigns were carried out during winters of 2007-2010; these campaigns were synchronized with con-current satellites overpasses, to collect the data of snow density, snow wetness, snow depth, snow-pack temperature and snow grain size parameters at various locations in Manali sub-basin.Three permanent snow and hydro-meteorological stations (refer to Fig. 2) at Bahang, Solang and Dhundi, of Snow and Avalanche Studies Establishment (SASE) at Manali, were also used for snow data collection (Instruments include automatic weather station, snow gauge etc).This ground data of hydro meteorological stations along with three snow pits at each location, total of nine samples, were used for validation of retrieved snow density.

Algorithms for snow density retrieval
The current study has used two approaches for dry snow density estimation from SAR data.In first approach, modified algorithm of Shi and Dozier (2000) has been used.The main inversion model for snow density estimation, originally made for L-band copolarized data is given below, but has been modified for C-band SAR data at the level of coefficients calculations, by updating the wave number, k, as that of the C-band SAR data and is explained in detail in this section.The C-band SAR data based inversion model for snow density retrieval has been named as, Modified Shi Snow Density Inversion Model, MSSDIM.Assuming the dry snow pack has no significant volume backscattering at L-band (as well as C-band) as snow grains are much smaller than incident L and C-band wavelengths, no significant volume scattering can be generated by a snow pack (Shi and Dozier, 2000); therefore, total backscatter can be written as: Where, σ t pp is the total backscatter in a given polarization (pp), k 0 is the incident wave number at the air-snow interface in cm ground interface in cm −1 for the given radar frequency, θ i is the incidence angle at airsnow interface, θ r is the refractive angle in snow pack or the incidence angle at snowground interface.Note that the dielectric contrast (ε g /ε s ) should be used instead of ε g (relative dielectric constant of ground) and that k 1 should be used instead of k 0 when calculating the backscattering σ g pp at the snow-ground interface.In the present study, the wave number k 0 has been used as 1.13, corresponding to ENVISAT ASAR C-band frequency of 5.33 GHz, wavelength 5.6 cm and k 0 as 0.29 for L-band PALSAR data to calculate the coefficients of Eq. ( 2).In this study snow pack is assumed to be of single layer equivalent density, instead of multi-layer snow pack, which has density variations.Since most of natural terrain considered in study area has a small surface roughness (dry winter grassland and barren areas) and random surface slope, single scattering would dominate over multiple scattering in most situations.In this study, it should be noted that at C-band (ASAR frequency), snow density affects the magnitude of the volume scattering and the surface scattering properties at the snow-ground interface, this would affect the overall accuracy of estimated snow density when using C-band data, but accuracy should be better with L-band PALSAR data.Finally, the algorithm for estimation of snow density, which uses only σ t hh and σ t vv SAR measurements, is given as In above equation, T pp is Fresnel transmission coefficients and depends on the polarization pp, the incidence angle θ i at the air-snow interface, and the dielectric constant of the snow-pack ε s , ε s is the only unknown in Eq. ( 2) and θ i can be calculated from a combination of the ASAR orbital data and a digital elevation model.The coefficients a to e have been calculated for full image by using the modified wave number, k 0 as 1.13 corresponding to C-band ASAR data, as compared to 0.29 value for L-band, as used in original formulation of Shi and Dozier (2000) and also used in present case with PALSAR data.Therefore, for a given C and L-band SAR measurements of σ t hh and σ t vv , ε s can be numerically estimated by varying the coefficients a, b, c, d and e to find the root of Eq. (2).It does not require a priori knowledge of the dielectric and roughness properties of the soil under the snow, but accuracy of this inversion model will be higher for L-band as compared to C-band.Furthermore, snow density can be estimated from Looyenga's semi-empirical dielectric formula (Looyenga, 1965), which provides a good fit to Polder and van Santen's physical formula (Matzler, 1996).
The second approach for estimating the snow density from C-band SAR data uses the ASAR-APS dual polarization data (Snehmani et al., 2010).In this approach for developing an algorithm for snow density estimation, the volume scattering model and the small perturbation model have been used with an exponential correlation (Ulaby et al., 1986) function for the surface backscattering contribution from the snow-ground interface.The final inversion model (Snehmani et al., 2010) is given in Eq. ( 4) and it is a function of only dielectric constant of the dry snow.
Where, α vv and α hh are the same as in Small Perturbation Model (Cloude, 1992).T vv and T hh are Fresnel transmission coefficients.By knowing the incident angle, Eq. ( 4 can be used to estimate the dielectric constant of the snow, which can be directly related to snow density using Looyenga's semi-empirical formula (Looyenga, 1965).The accuracy of snow density results with respect to ground measurements has been done using the coefficient of determination or R 2 (Steel and Torrie, 1960).The coefficient of determination is a measure of how well the regression line represents the data (Nagelkerke, 1991).If the regression line passes exactly through every point on the scatter plot, it would be able to explain all of the variation.

Results of snow density inversion model
The snow density maps derived using modified and original Shi and Dozier (2000) algorithm (MSSDIM, approach-I) are shown in Fig. 4a-d, and the approach-I and II results are shown in tabular form in Table 2.The dry snow area estimated using approach-II for 20 January 2008, comes out to be 8 % of total area, with snow density ranges from 0.01 to 0.2 g cm −3 , as compared to approach-I, where it is estimated as 33 %.Similarly, the high snow density (> 0.51 g cm −3 ) area in approach-II was 32 to 39 %, but for approach-I it was nil for all the analyzed images (Table 2).The ground observations show that mean snow density of snow pack at Dhundi during 20 January 2008, 25 January 2008 and 20 January 2009 varied from 0.09 to 0.15 g cm −3 and snow depth varied from 235, 196 and 120 cm, respectively.Therefore, the snow density retrieved from approach-I match better with ground observations, as compared to approach-II, where snow density is over-estimated.The results of ALOS-PALSAR for 9 December 2009 are also comparable with ground observations, but due to low incidence angle (25.65 • ) of PALSAR data and partial coverage, large area has been masked out in this map.Overall, the ASAR-APS data based snow density is overestimated using approach-II as compared to approach-I (table 2).The over-estimation is also evident from the fact that it has the large area in density class > 0.4 g cm −3 and also large area under wet snow class (snow density > 0.5 g cm −3 ) at higher elevations, whereas air temperature in these areas is less as compared to lower elevation area.Therefore, this approach needs to be refined to work at entire basin scale, mainly by improving formulation of snow-ground scattering and volume scattering.
The results of approach-I are less similar as compared to the original model results of Shi and Dozier (2000).In case of approach-I (MSSDIM), as C-band has been used, it has added to slightly more snow-ground scatter and volume scatter, which has resulted in less accuracy, but the fully polarimetric L-band PALSAR data (9 December 2009) results of snow density are comparable to original results of Shi and Dozier (2000).Still, the accuracy of approach-I for snow density estimation is relatively higher as compared to that of approach-II, as it is also evident from the Table 2, where approach-I has majority of snow density classes in the range of 0.1 to 0.2 g cm −3 .

Validation and accuracy of retrieved snow density
The snow pit based ground sample of Dhundi, Solang and Manali are used for accuracy assessment of snow density (Fig. 6).The overall coefficient of determination, R 2 for approach-I is 0.85 and for approach II is 0.72.The accuracy in terms of root mean square error (r.m.s.e.) using approach-I comes out to be 0.03 g cm −3 .The R 2 calculated here only represents the ground observation vs. derived values at three sites, but maps have been created for full watershed area or part of image.The snow density derived from SAR based inversion models could not be verified at higher elevation areas.This is due to the fact that there is no accessibility at these higher altitude areas in winters due to closure of roads in heavy snowfall conditions, and there are no ground observation sites at elevations more than 3500 m above mean sea level.This is one of the limitations of this study.For approach-II, the accuracy of results are less as compared to accuracy obtained by Snehmani et al. (2010).However, the later study does not reflect overall picture, as entire model was run on full ASAR image, without forest/glacier mask and ground points were taken only in few selected areas, which are same as shown in Fig. 1.Therefore, the results of both the approaches of dry snow density further need to be checked with more ground data uniformly distributed at all elevations and land use classes to give more accurate results and its further use in snowmelt runoff and snow avalanche models.Some of the area is masked out due to the presence of high relief and steep slopes present in study area, which makes this watershed area outside the permissible range of local incidence angle (10-70 • ) of snow density inversion model.The relative contribution of layover, shadow and forest masks is 19 to 27 % in case of snow density maps (approach-II, Table 2), where no mask for local incidence angle has been used.
For snow density inversion model, MSSDIM, the mask area varies from 28 to 48 % with 56 % for PALSAR data (as only medium incidence angle image is available).Thus, there is a strong need for creating a separate forest backscatter model for study area (to map snow and retrieve snow density below forest, e.g., Koskinen et al., 2010), along with inversion models, which have higher range of applicability in terms of local incidence angle.The layover, shadow effects are minimized by selecting the higher incidence angle (> 40 • ) SAR data.Overall, the snow density is under-estimated at higher elevations and over-estimated at lower elevations, mainly by inversion model based on approach-II.

Variations of snow density with elevation and aspect zones
The variations of snow density with various zones of elevation and aspect are shown in sample Fig. 5a and b for 20 January 2008.It has been observed from this analysis that majority of snow density at all elevations and aspects zones comes in classes 1 and 2, i.e., 0.06-0.2g cm −3 .The image of 25 January 2008 shows higher snow density classes (0.3-0.5 g cm −3 ) at all elevations zones, mainly in higher areas with northern

Limitations of the present study
The present study has used masking for excluding the areas of forest and glaciers.This has been done as the current inversion models have not been formulated to work in all surface conditions.This is one of the limitations of the present study.In future, if a separate forest and glacier backscatter models for the study area can be created in C, L or X-band SAR data, then a backscatter compensation for these Land Use Land Cover (LULC) categories can be applied to improve the range of applicability and accuracy of snow parameter inversion models.The snow depth could not be retrieved from SAR data in this study due to insufficient ground data on snow depth and nonavailability of multi frequency SAR data of same date.
Another limitation of the current inversion models is that it considers average properties of entire snowpack, assuming it to be a single snow layer model.The results of current models for snow density give average value for the snowpack.In reality, snowpack has variable snow properties from top of layer to bottom of snowpack (Wold, 1986).The density of snowpack varies with depth of snow layer, thereby also changing the SWE of snowpack.Therefore, improvement needs to be made in backscattering models of wet as well as dry snowpack, by including the multi-layer snow properties as part of surface and volume backscatter models.

Applications of derived snow density maps
The snow density is one of the most useful snow parameter which is of great importance for snow avalanche and snowmelt runoff modeling studies.Snow density is one of the inputs for the estimation of SWE along with snow depth.The area from Solang to Dhundi snow observatory has 11 known avalanche sites (Sharma and Ganju, 2000) and also this area has variable snow cover varying from 15 % to 100 % which results in significant snowmelt runoff in spring season.The optimal snow density for slab avalanche formation is in the range of 0.10 to 0.25 g cm for loose snow avalanches.Similarly the derived snow density can be used to find the degree day factor which is main parameter in snowmelt runoff models (Martinec, 1960;Rango and Martinec, 1995).When the snow density increases, the albedo decreases and the liquid water content in snow increases.Thus the snow density is an index of the changing properties which favor the snowmelt (Rango and Martinec, 1995).Overall, the results of present study can be used for snowmelt runoff modeling and snow avalanche modeling in selected areas, where model validation has been done.

Conclusions
For estimation of dry snow dielectric constant and snow density in Himalayan region, the inversion model is applied on radar backscattering coefficient.For solving complex equations, on each incidence angle, codes have been developed in MATLAB, Mathematica and C++.No incidence angle mask was used for snow density model with approach-II.This study concludes that ASAR-APS based inversion models tends to overestimate the snow density, mainly in areas of higher elevations, and further ground truth and refinement of inversion models are needed for their operational use by user's agencies and government departments.The present study also concludes that the modified model of Shi and Dozier, (2000), named as MSSDIM, along with original model (approach-I) can be used for estimation of snow density using dual polarization (HH/VV) C-band and L-band SAR data respectively.The MSSDIM model gives better results as compared to semi-empirical approach (approach-II).In the MSSDIM method of snow density estimation, L-band data gives slightly better results than C-band data.
The MSSDIM further needs to be tested with more L-band data of study area, as well as refinements and more modifications in original coefficients of the inversion model with C-band SAR data.These changes are required to improve the accuracy of estimated snow density and its operational usability in hydrological and snow avalanche models.In future, the effect of vegetation cover and underlying glacier on mapping of

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | snow is mainly from surface and volume scattering only, whereas snow-ground surface scatter dominates in the dry snow pack.The terms A, B and C of Fig. 1 are, (A) backscattering from the snow-air interface, (B) volume scattering from the snow pack and, (C) backscattering from the underlying ground surface.

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | −1 , k 1 is the incident wave number at the snow-1933 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ε s ) Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

) 1935 Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper |

1937
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | aspects and relatively medium density for other areas.The other images have majority of snow density class in 0.1 to 0.2 g cm −3 .Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

− 3 .
Field data show that slabs also occur outside this density range.Low-density dry snow can be responsible 1939 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 4 :
Fig. 4: Snow density maps derived from approach-I (MSSDIM method) for study area for (a): Snow density map for 20 th January 2008

1957Fig. 5 :
Fig. 5: Snow density variation with Elevation Zones (EZs) of study area The Beas River rises from Rohtang pass 4350 m above mean sea level, 51 km north of Manali.The Beas Kund glacier is the main source of water for Beas River upto Manali.Climate of the study area is cool and dry, with cold season