ESTIMATION OF SURFACE SNOW WETNESS USING SENTINEL-2 MULTISPECTRAL DATA

Snow cover characterization and estimation of snow geophysical parameters is a significant area of research in water resource management and surface hydrological processes. With advances in spaceborne remote sensing, much progress has been achieved in the qualitative and quantitative characterization of snow geophysical parameters. However, most of the methods available in the literature are based on the microwave backscatter response of snow. These methods are mostly based on the remote sensing data available from active microwave sensors. Moreover, in alpine terrains, such as in the Himalayas, due to the geometrical distortions, the missing data is significant in the active microwave remote sensing data. In this paper, we present a methodology utilizing the multispectral observations of Sentinel-2 satellite for the estimation of surface snow wetness. The proposed approach is based on the popular triangle method which is significantly utilized for the assessment of soil moisture. In this case, we develop a triangular feature space using the near infrared (NIR) reflectance and the normalized differenced snow index (NDSI). Based on the assumption that the NIR reflectance is linearly related to the liquid water content in the snow, we derive a physical relationship for the estimation of snow wetness. The modeled estimates of snow wetness from the proposed approach were compared with in-situ measurements of surface snow wetness. A high correlation determined by the coefficient of determination of 0.94 and an error of 0.535 was observed between the proposed estimates of snow wetness and in-situ measurements.


INTRODUCTION
Snow cover characterization and estimation of snow geophysical parameters are significant for understanding the hydrological budget of the glacial rivers and for studying the snow surface melt processes.Alpine snow and glaciers constitute a significant part of the cryosphere which yields through melt runoff one of the major usable water resources.Particularly, for countries like India which rely significantly on glacier and snowmelt runoffs for water resource, the understanding of snowmelt processes is significant.Forecasting of snowmelt runoff requires timely information on the spatial-temporal distribution of the geophysical parameters of snow such as the liquid water content (LWC), density of snow, and the depth of the snowpack (Shi and Dozier, 1995).The continuous monitoring of snowpack variables is also significant in avalanche forecasting.In the Himalayas, due to higher elevation and difficult terrains, often several areas remain inaccessible, especially during the winter season, imposing several constraints on the possibility of field campaigns for in-situ measurements of snow geophysical parameters.
With the advances in the remote sensing technology, monitoring of snowpack state variables (snow geophysical parameters) such as the extent of snow commonly known as the snow cover area (SCA) maps, density, depth etc. has been possible.These variables are widely used in the snowmelt runoff model (SRM) for seasonal forecasts on runoff (Rango and Martinec, 1979).Active microwave remote sensing data based on fully polarimetric synthetic aperture radar (SAR) systems have been significantly incorporated into physical scattering models for estimating snow wetness at C-band (Shi and Dozier, 1995;Surendar et al., 2015).Singh and Venkataraman (2009) used the dual-polarimetric ASAR C-band SAR data for the determination of snow dielectric constant using a SAR based inversion model and estimated the snow density using the Looyenga's formula (Looyenga, 1965).Polarimetric decomposition techniques have also been efficiently utilized for estimation of snow wetness * Corresponding author using fully polarimetric Radarsat-2 SAR data by Surendar et al., (Surendar et al., 2015).
With optical remote sensing, quantitative studies on snow geophysical parameters have been limited due to the constraints posed by the optical characteristics of snow and from those inherent to the sensors (Hall et al., 2006).Snow exhibits unique properties in the visible, infrared and thermal spectrum.The metamorphic state, age and level of contamination largely determine the spectral response of snow.Freshly deposited snow exists in a dry powdered form with a very small grain size <1mm and very low LWC <<1% by volume.Dry snow exhibits very high reflectance in the visible and near infrared (NIR) spectrum.The reflectance reduces as the snow ages and crystallizes after constant processes of melt and refreeze which results in an increase in snow grain size and compactness of the snowpack.In the NIR spectrum, snow grain size with respect to reflectance shows a higher sensitivity.This sensitivity is also attributed to the LWC in snow (Dozier et al., 2009).A relation between the spectral reflectance of Landsat Thematic Mapper and the snow grain size was defined by Dozier (1989).The snowpack is wet when LWC increases (LWC>1%) and becomes slushy when the LWC is much higher (LWC>10%) for the snow remains in crystallized form.
The near surface moisture index (NSMI) was proposed by Lampkin and Yool (2004) for the qualitative assessment of snow surface wetness.NSMI incorporates the normalized differenced snow index (NDSI) and thermal brightness temperature and was represented as a function of snow grain size based on the Moderate Resolution Imaging Spectroradiometer (MODIS).Previously, a similar concept was used in the assessment of soil moisture based on the triangle method, where a triangular feature space was developed using land surface temperature and the normalized differenced vegetation index (Nemani et al., 1993).However, the application of these approaches is dependent upon the availability of the thermal channel in the multispectral data.
These approaches based on the thermal channel are also dependent upon the atmospheric composition.The NIR reflectance has been effectively used by Gupta et al. (2005), for mapping wet snow using the Indian Remote Sensing Satellite (IRS-1C/1-D) with Linear Self-scanning Sensor (LISS-III).The hard threshold defined for dry snow was found greater than or equal to 0.5 of the NIR reflectance, attributed to the lower dynamic range of the LISS III sensor.However, quantitative assessment of surface snow wetness using optical data remains limited.
In this paper, we elaborate on the methodology presented by Varade and Dikshit (2017), for the estimation of snow wetness using a generalized linear physical model which relates the snow wetness and the NIR reflectance.We present a triangular feature space generated by incorporating NDSI and NIR reflectance.The longer edges of the triangular feature space represent the dry and wet edges corresponding to the upper and lower edges of the space, respectively.In the proposed methodology these edges are identified based on visual perception.These edges are used to develop an empirical model for the estimation of snow wetness.

METHODOLOGY
2.1 Selection of suitable band for the proposed feature space.
The significant variation in the spectral response of snow in the visible and shortwave infrared (SWIR) spectrum is utilized in the NDSI for the discrimination of snow.Higher NDSI corresponds to pure snow pixels with finer grain size and vice versa.However, NDSI is a qualitative measure, which depends upon several factors such as the illumination angle, existence of soot etc. especially in case of mixed pixels which are often plenty in medium resolution spaceborne multispectral data.The spectral reflectance of snow in the visible spectrum is high and nearly uniform for both wet and dry snow.However, in the NIR spectrum dry snow has high reflectance and wet snow has lower reflectance depending upon the LWC of snow surface (Lavan Kumar et al., 2017).Thus, the proposed methodology is based on utilizing the NIR reflectance instead of the brightness temperature as proposed by Lampkin and Yool (2004).
Although a similar observation is seen in the spectral response of snow in the SWIR spectrum, the range of reflectance for each class corresponding to dry snow, wet snow and snow-free areas is rather compressed.It is evident that for selection of the specific band for the proposed triangular feature space, we must confirm its capability in discrimination of different land cover classes especially the subclasses of snow i.e. dry and wet snow, which is essential for the determination of snow LWC. Figure 1 delineates the spectral response of different classes of pixels observed in the study area including dry and wet snow.The pixels marked 'D0' and 'W0', correspond to the verified pixel with in-situ measurement.The other pixels were selected based on visual interpretation and other auxiliary data.Figure 1 illustrates that the spectral response of dry snow is the highest amongst all the classes, while that of wet snow varies.The reflectance for different classes is separable in bands B8 and B8A corresponding to 842 nm and 865 nm channels of Sentinel-2.At the SWIR bands, the reflectance significantly reduces, where the wet snow exhibits very low reflectance due to higher absorption in SWIR spectra (Hall et al., 1995;Nolin and Liang, 2000).and B11 and B12 with respect to NDSI

NIR-NDSI feature space
From Figure 3 it is evident that dry snow exhibits higher NDSI values and vice versa for wet snow.There are two major issues that need to be dealt with while defining a model for the estimation of snow wetness.Firstly segregation of snow pixels with pixels in the snow-free area and secondly segregation of dry and wet pixels in the snow cover area where NDSI represents mixed values.In general, for snow covered area NDSI exhibits values greater than 0.4.However, depending upon numerous factors such as the illumination angle, level of contamination etc. the threshold may change.A simple method to decide the threshold value for snow is to identify patches of snow and snowfree areas and analyze their histograms.NDSI greater than 0.7 represents snow with finer grain sizes which are often the case with dry snow (fresh powdered or crystallized refrozen).
The analysis in the previous section can be summarized into the statement that spectral response in the visible and infrared regions may be used for modeling surface snow wetness (LWC).
To facilitate this, we propose the development of a feature space based on the observations of NDSI and NIR reflectance, which we refer to as the "NIR-NDSI feature space".A generalized layout of the NIR-NDSI feature space is shown in Figure 4.For this feature space, we utilize the B8A band at 865 nm wavelength of Sentinel-2 which was observed to be a better choice as observed from the analysis in the previous section.The feature space is developed using a scatter distribution of NDSI and B8A reflectance on the horizontal and vertical axes, respectively.
For modeling the snow wetness, we utilize a similar approach as proposed by Sadeghi et al., (2017).The significant components of the NIR-NDSI feature space are the dry and the wet edges.In case of this feature space, the edges are in the reverse order as compared to the conventional feature space used in the estimation of soil moisture (Nemani et al., 1993;Tang et al., 2010;Sadeghi et al., 2017).The cluster of pixels corresponding to snow-free, dry snow and wet snow are also shown in Figure 4.The pure pixels are identified at the rightmost part of the feature space.In this case, the pure snow pixels are characterized as dry and wet snow pixels existing at the top and bottom, respectively, of the rightmost part of the feature space.This follows from the fact that, although pure snow pixels have higher NDSI value, the wet snow pixels exhibit lower NIR reflectance than dry snow pixels.
In the case of mixed pixels, the NIR reflectance is reduced due to contributions from other adjacent classes such as trees and vegetation or due to contamination in the snow.As a result, not only the NIR reflectance is reduced, but melting is also enhanced, leading to an increase in the LWC of the snow.These mixed pixels are defined by mid-range NDSI with relatively lower NIR reflectance.The snow-free areas are characterized by lower NIR reflectance and lower NDSI values and appear at the leftmost end of the NIR-NDSI feature space.
Figure 4.The generalized layout of dry and wet snow pixels in the NIR-NDSI feature space.

Selection of the edges
Due to the topography in alpine terrain such as in the Himalayas, shadows often exist.In this case, shadowed pixels largely affect the scatter distribution between the NIR reflectance and NDSI.
The distortion that may arise due to the existence of shadows is also attributed to the presence of clouds.These pixels contribute to enlarged wet pixel zones where the NIR reflectance would be lower and the NDSI would still show these pixels within the snow cover area.A solution to this would be to mask these pixels.
In a similar manner, standing water in the study area would not ideally absorb the incident radiation in the NIR spectrum completely.This results in corresponding pixels in the wet zone with very low NIR reflectance, leading to an oversaturated wet edge.Due to the existence of oversaturated wet zones, automatic edge detection techniques (Tang et al., 2010) techniques are not applicable for the proposed feature space, which also implies that the proposed feature space is sensitive to oversaturated pixels.The preferred method, in this case, is the manual selection of wet and dry edges.A study suggested that visual inspection of pixel distributions provides the optimum edges (Toby N. Carlson, 2013).

Snow wetness modeling
The surface snow wetness can be modeled from the dry and wet edge values of NIR reflectance as shown in equation 1.  between snow wetness and NIR reflectance is based on the assumption that the sensitivity of NIR reflectance is a linear function of the snow wetness.
The parameters in equation 1, can be obtained based on the dry and wet edges of the feature using their corresponding slope and intercept.The wet and dry edges can be represented by the following equations of straight lines.
where Rw and Rd represent the reflectance at θw and θd, respectively.The parameters id and iw indicate the intercepts of the dry and wet edges respectively, and the parameters sd, and sw indicate the slopes of the dry and wet edges, respectively.The snow wetness can be estimated by combining equations 1, 2 and 3 as follows.
( ) ( ) In contrast to NSMI (Lampkin and Yool, 2004), the proposed model in equation 4, characterizes a much stable model for ascertaining the LWC in the top layer of the snowpack and is expected to remain nearly time and atmosphere invariant, since, the reflectance is a property of the incident surface.

Study area and data description
A study area from Dhundi to Palchan was selected along the Beas river in Himachal Pradesh, India, geographically centered at a latitude and longitude of 31° 59' 32" N and 76° 44' 49" E, respectively.A corresponding Sentinel-2 imagery was available through the Copernicus program of the European Space Agency (ESA), with the date of acquisition on 9 th February 2017, as shown in Figure 5.The band digital numbers of the Sentinel-2 imagery were converted to bottom of atmosphere reflectance using atmospheric correction based on the Sen2cor L2A processor (Louis et al., 2016).The multispectral bands were resampled to a common spatial resolution of 20 m.A field campaign was conducted in the study area during 11 th -12 th February, 2017, where in-situ measurements using SnowFork (Techel and Pielmeier, 2011) were obtained at different locations.
The geographical location of in-situ measurement points is shown in Figure 6 with pins, where green and red pins correspond to dry and wet snow observations, respectively.The average values for the various in-situ measurements are shown in Table 1.To maintain consistency in the observations redundancy was maintained by collecting simultaneous 5-7 observations at each point.

NIR-NDSI feature space
The feature space generated using the band B8A and NDSI from the test data is shown in Figure 7.As discussed previously the feature space contains oversaturated pixels which result in an overall distortion of the shape of the space.The distortion is less significant in case of manual selection of the dry and wet edges which was the method used in this study.The slope and intercept of the dry edges were selected to be 0.08 and 0.7444, respectively.The slope and intercept of the wet edge were selected to be 0.08 and 0.06, respectively.The corresponding points for determining the extended edges were identified to begin near the vertex where the two edges meet and on the right side by manually inspecting the density of the points.From experiments, it was identified that extended edges tend to introduce less ambiguity in the estimation of snow wetness.
Figure 7. NIR-NDSI feature space using B8A reflectance in the NIR spectrum.

Qualitative assessment of wet and dry snow wetness
As discussed earlier due to the topography in mountainous regions, many shadowed areas are observed as shown in Figure 8 (red ellipses).The shadowed areas exhibit very high LWC in the generated snow wetness map from the proposed methodology.Due to the similarity in reflectance between shadowed areas and the lower valley areas where snow is soiled or contaminated, masking of shadows was not possible.However, in the snow wetness map, they can be easily removed as observed from Figure 8.In practice, we can simply ignore the shadowed areas as identified in Figure 8, since the observed wetness is understood to be erratic within these areas.However, for some regions, the partially shadowed areas show a similar magnitude of LWC as some of the wet areas (4-6%, as observed from in-situ observations).The snow wetness map for the proposed feature space shows good compatibility with the field conditions, indicating a good in wetness for dry and wet snow regions, with relatively higher values for wet snow and lower wetness values for dry snow.This is evident from Figure 8, where within the valley areas the wetness observed is higher and much lower for the snow deposits on the mountain slopes.
Figure 8. Shadowed areas overlaid on the false color composite of the test data and on the snow wetness map generated from the proposed methodology.

Validation of estimated snow wetness
A comparison of the estimated snow wetness was carried out with respect to the in-situ measurements collected during the field campaign.Figure 9 shows the correlation plot for the estimated and the observed snow wetness.A well-clustered distribution of points is observed in the correlation plot shown in Figure 9.The dry snow pixels are clustered around the lower left corner of the plot, while the wet snow pixels are clustered around the upper right corner of the plot indicating a good agreement with LWC aspects of snow.The linear fit between the estimated and observed values of snow wetness is modeled as follows.
1.213 1.08 oe ww  (5) where wo and we represent the observed and estimated values of snow wetness.The linear fit also confirms the assumption that NIR reflectance varies linearly with respect to the LWC in the snow.A coefficient of determination (R 2 ) of 0.94 was observed with a root mean square error (RMSE) of 0.535 between the estimated and observed snow wetness at the in-situ sampling sites.Although a good correlation is evident from the segregation of the clusters, the intrinsic correlation between the points within the cluster also shows a good agreement with R 2 0.732 for the rightmost cluster and 0.432 for the leftmost cluster, representing wet and dry/refrost snow respectively.Due to the unavailability of the NIR bands in the 980 nm -1030 nm range, the discrimination between wet and refrost snow is not possible with multispectral sensors (Dozier et al., 2009).The refrost snow should contribute to the cluster with lower wetness (leftmost cluster) and may have caused a reduction in the observed correlation for the dry snow measurements.In this study, although limited in-situ observations were possible for the study area, the estimated snow wetness was observed to comply with the expected field conditions around Dhundi and Solang, in Himachal Pradesh, India for the month of February.
During the collection of possible in-situ measurements, snow wetness around 3%vol was not observed.Thus, Figure 9 comprises a gap in points around 3% vol wetness.A reduction in the observed correlation is certainly possible, given a larger number of in-situ measurements.However, since the results comply with the overall field conditions, we believe the proposed approach has potential in delivering significant results.We must also mention that in-situ collection of data depends upon numerous factors including availability of funds, manpower, the complexity of the terrain, accessibility and the immediate weather conditions.

CONCLUSION
In this study, we present a linear physical model based on a feature space derived from the observations of NIR reflectance and NDSI.The proposed approach utilizes the linear variation of spectral reflectance in the NIR spectrum with respect to snow wetness.The proposed method is invariant to atmospheric factors since it utilizes the NIR reflectance instead of other variables such as the land surface temperature which is significantly affected by atmospheric composition.The proposed methodology was evaluated with Sentinel-2 multispectral data.
However, the proposed model can be easily extended to several other multispectral sensors including unmanned aerial systems (UAS), which may play a key role in monitoring snowpack conditions for avalanche forecasting.The experiments based on the proposed approach indicated sensitivity towards natural obstruction as observed from the effect of shadows in the results.However, in the case of UAS surveys, the approach may find its true potential since UAS based collection of data is expected to be relatively less obstructed by the presence of clouds and shadows.Overall in the unobstructed area, the results from the proposed approach showed good agreement with the expected field conditions.A good correlation determined by R 2 of 0.94 and RMSE of 0.535 was observed between the in-situ measurements of snow wetness and the estimated snow wetness.

Figure 2 .
Figure 2. Statistical parameters for the different bands.Due to overlapping response, another assessment for suitability of bands is carried out using statistical measures such as variance, skewness, kurtosis as shown in Figure 2. In Figure 2, SSp denotes the root mean squared sum of the parameters skewness, kurtosis and variance, and NSSp denotes the normalized values of SSp.Based on the values of NSSp B8A is apparently a better choice for selection.Although the SWIR bands B11 and B12 show values capable of discriminating dry and wet snow based on

Figure 3 .
Figure 3. Spectral response of different classes in bands B8A and B11 and B12 with respect to NDSI the snow wetness or the LWC in the top layer of the snowpack, θd is the local minimum LWC, θw is the local maximum LWC and θ is the local LWC in the top layer of the snowpack.The snow wetness can be related by the NIR reflectance as shown in equation 1, where R NIR , represents the local reflectance in the NIR spectrum.NIR d R and NIR w R represent the corresponding NIR reflectance at the dry and wet edge respectively in the NIR-NDSI feature space.This relation 0

Figure 5 .
Figure 5. Study area and test dataset

Figure 9 .
Figure 9. Correlation between estimated snow wetness from the proposed methodology and observed snow wetness.