Multi-Sensor Analysis of Snow Seasonality and a Preliminary Assessment of SAR Backscatter Sensitivity to Arctic Vegetation: Limits and Capabilities

Snow melt timing and the last day of snow cover have a significant impact on vegetation phenology in the Svalbard archipelago. The aim of this study is to assess the seasonal variations of the snow using a multi-sensor approach and to analyze the sensitivity of the Synthetic Aperture Radar (SAR) backscatter to vegetation growth and soil moisture in an arctic environment. A combined approach using time series data from active remote sensing sensors such as SAR and passive optical sensors is a known technique in snow monitoring, while there is little knowledge of the radar Cband’s response pattern to vegetation dynamics in the arctic. First, we created multi-sensor masks using the HV backscatter coefficients from Sentinel-1 and the Normalized Difference Snow Index (NDSI) time series from Sentinel-2, monitoring the snow dynamics in Adventdalen (Svalbard) for the season from 2017 to 2018. Second, radar sensitivity analysis was performed using the HV polarized channel responses to vegetation growth and soil moisture dynamics. (1) Our results showed that the C-band radar data are capable of monitoring the seasonal variability in timing of snow melting in Adventdalen, revealing an earlier start by approximately 20 days in 2018 compared to 2017. (2) From the sensitivity analyses, the HV channel showed a major response to the vegetation component in areas with drier graminoid dominated vegetation without water-saturated soil (R = 0.69). However, the temperature was strongly correlated with the HV channel (R = 0.74) during the years with delayed snow melting. Areas of frozen tundra with drier vegetation dominated by graminoids had delayed soil thawing processes and therefore this may limit the ability of the radar to follow the vegetation growth pattern and soil moisture.


Introduction
The timing of snow melt and the first day free of snow are considered indicators of Arctic climate and ecosystem status in response to global warming [1]. Moreover, the depth of snow and the period of snow melt contribute significantly to defining the phenological phases of vegetation [2,3] and plant biomass [4]. Due to an increase in winter precipitation [5], as well as in the frequency of extreme weather events in the Svalbard archipelago [6,7], snow cover monitoring is highly important.
Along with the ability to detect snow seasonality, there is a growing interest in studying the relationship between the timing of snow melt and vegetation phenology (late snow, onset of the growing season, and advanced snow melting). Malnes et al. (2016) [8] used satellite remote sensing (MODIS) over northern Norway in order to detect the snow covered season with significant accuracy due to start (r = 0.51, p > 0.05) and end (r = 0.79, p > 0.05) of the snow covered season for most meteorological stations for the 2000-2010 period. However, in some areas and some years, the snow covered season could not be detected due to long overcast periods.
Since cloud cover often dominates the Arctic environment and especially in Svalbard, it would be useful to monitor phenology with sensors that are unaffected by clouds. Satellite remote sensing imagery can be used in the analysis and interpretation of the snow-pack variations of Svalbard [9]. Multi-sensor analysis has proven to be effective to monitor snow coverage and its relation to vegetation phenology [10,11]. In the optical domain, the Normalized Difference Snow Index (NDSI) has been successfully used to detect snow cover [12,13]. This index takes advantage of the high reflectance of snow in the visible wavelengths and low reflectance in the short wave infrared wavelengths. While optical sensors are based on surface reflectance, SAR sensors allow evaluating change that is occurring within the surface. This is because the strong dielectric contrast between the solid and liquid phases of the water generates changes in the backscatter coefficients. Indeed, even a small amount of liquid water reduces the depth of penetration of the radar signal, allowing identifying the start of the snow melting process [14].
The most common approach to snow mapping with Synthetic Aperture Radar (SAR) is based on multi temporal comparison of images of the same area in snow-free/dry and wet snow conditions [15][16][17][18]. Attenuation due to snow in the frequency range 1-12 GHz is very low, and dry snow covers cannot be discerned from the bare ground. Conversely, there is a signal attenuation in wet snow conditions, due to a change in the dielectric properties of the surfaces [16,19,20].
Due to polar night, SAR data are essential in the arctic during the winter period. However, because of the high cloud coverage that limits optical images, SAR data are also important during the summer period. Once the snow has disappeared, the radar beam can simultaneously penetrate both the vegetation and the soil to a depth that is difficult to determine [21,22]. Innovative approaches are based on the use of SAR and optical sensors to follow the phenological phases of vegetation [23]. These approaches may be useful in areas where cloud cover is widespread during the growing season. Moreover, the sensitivity of backscattering coefficients to the vegetation of high arctic not yet been investigated to our knowledge.
The Copernicus program allowed free access to time series of different sensors, which can be used in synergy. Sentinel-1A and 1B SAR sensors (centre frequency of 5.405 GHz) provide medium and high-resolution time series of C-band data, while Sentinel-2A and 2B optical sensors acquire 13 spectral bands in the optical domain [24,25]. The two sensors used simultaneously can improve the characterization of the snow season and vegetation growth, allowing us to derive multi-sensor products with a high temporal and spatial resolution. The central aim of this study was to derive wet snow maps and snow maps from the time series of S-1 and S-2, in order to assess the seasonal variations of the snow. Then, a sensitivity analysis of the backscattering coefficient σ 0 to vegetation and soil moisture was made and related to snow dynamics. The final phase enabled an overall understanding of the sensitivity of the SAR signal to vegetation in relation to snow dynamics. The novelty of our study is to assess the impact of vegetation growth and soil moisture on the SAR signal in the Arctic, taking into account the variability of the snow seasonality.
Compared to the studies presented above, our main objectives of the study are: 1.
to apply a combination of radar and optical satellite data (Sentinel-1 and Sentinel-2) to map the spatial and temporal pattern of wet and dry snow conditions, and its relation to the vegetation growth season; 2. to understand how polarized radar data (HV-horizontal transmit and vertical received) can contribute to detect the pattern of arctic vegetation growth and soil moisture.

Study Area
The study area is the Adventdalen valley and the surrounding plateaux, close to the town of Longyearbyen (78 • 13 N 15 • 33 E), in the Svalbard archipelago ( Figure 1). The study area lies approximately between latitudes 78 • 20 and 78 • 07 N, and longitudes 15 • 10 and 17 • 10 E. The periglacial landscape is characterized by vast plateaux intersected by wide glacial valleys and alluvial plains. Precipitation is low and the maritime environment strongly influences the snowpack characteristics [26]. Adventdalen is characterized by a polar-tundra climate [27] and is located in the Middle-Arctic and Northern-Arctic tundra zone. The Middle-Arctic zone is characterized by dwarf-shrub heaths, where Cassiope tetragona often dominates, and in the study area small patches of Betula nana are found, whereas in the Northern Arctic Tundra zone, the genus Luzula is characteristic with Salix polaris, Saxifraga oppositifolia and Dryas octopetala [28,29] also common. The meteorological station located in Adventdalen (Station number SN99870, Norwegian Meteorological Institute) recorded an average air temperature in July of 6.8 • C and 7.2 • C for the years 2017 and 2018, respectively. Within the projects SnoEco (NRC ref. 230970), Sentinels Synergy Framework (EC FP7 collaborative project), and SIOS (www.svalbard-sios.org, accessed on 18 September 2019), ten ground stations were set up in Adventdalen. Each station [30] was equipped with data loggers, environmental sensors for soil temperature and moisture at 10 cm depth and time-lapse ordinary RGB cameras and NDVI sensors positioned at 2 m above the ground. Since the location of some stations changed between 2017 and 2018, we selected only five stations; their characteristics are listed in Table 1, while an example of images is shown in Figure 2. Camera images covered an area of approx. 1.4 square meters.

Datasets
The analyzed data sets are composed of time series from Sentinel-1, Sentinel-2, and ground station data. 115 S-1A and S-1B images acquired from February 2017 (4 February 2017) to December 2018 (26 December 2018) were processed. With a revisit frequency of six days, the C-band (wavelength, λ = 5.5 cm) images were acquired in Interferometric Wide swath mode (IW), with one relative orbit (track 014), and ascending pass. The images were available in cross polarized HV ('Horizontal transmit' and 'Vertical receive') channel, with a spatial resolution of 10 m. The S-2 Multi Spectral Instrument (MSI) acquires 13 spectral bands in the Visible, Near-Infra-Red (NIR) and Short Wave Infra-Red (SWIR) domains, with a spatial resolution of 10 to 60 m. Time series of S-2A and S-2B were processed starting from Level-1C Top-Of-Atmosphere reflectance. The high latitude results indeed in low solar elevation angles, which generate an under-correction of the atmospheric signal for the Bottom-Of-Atmosphere (BOA) processing level 2A [31] Due to polar night, images were acquired from April to September of each year. The polar orbit of S-2 enables daily data of Adventdalen to be obtained. However, only 43 images were processed in this study because of fog and frequent clouds, which reduced the number of images available. The ground stations were equipped with Meter Group SRS-NDVI Sensor [32] and soil temperature/moisture sensors (5 TM [33]). Water content and temperature of soil were measured at 10 cm depth. Station ST1 and ST3 additionally included infrared thermometers (SI-421 [34]) to monitor surface temperature. Data were recorded every four hours from mid-April to the end of September. The 5TM soil moisture and temperature sensors were not installed before the beginning of June, as the active layer was frozen until then. Time lapse cameras (WingScapes, model WCT-00122), with 8 MP of resolution, recorded images three times a day (9 am, 12 am and 3 pm). April to the beginning of October for three of the stations (6, 7 and 9), while station 1 recorded from mid of August to October, and station 3 from end of July to October)

Methodology
In this study, we primarily applied known techniques to define the melting of snow and the last day of snow. Following snowmelt, a study on the radar response to vegetation growth and soil moisture was carried out. The scheme of the suggested methodology is shown in Figure 3 and briefly described in the following sections.

SAR, Optical and Ground Data Processing
The pre-processing of S-1 data includes several standard steps to derive geocoded intensity images from Level-1 GRD (Ground Range Detected) data. Each scene was geocoded in the Norut software package GSAR [35] and stored as geotiff files. To reduce speckleaffecting backscattering values and obtain a more homogenous snow cover pattern, a Frost filter with a 7 × 7 window was applied to the images [36]. Top-of-atmospheric Sentinel-2 products (L1C) were processed to obtain clear-sky time-series. The cloud detection provided with the Sentinel-2 (L2A) product was often not sufficiently robust, and so clouds were detected by a visual analysis of the images combined with different cloud-detection algorithms for cloud removal. Specifically, we used algorithms from the literature [37,38], along with our own developed algorithms based on our experiences with time series processing in Svalbard [39,40]. After this step, in the period from late April to late September, 16-22 days were cloud-free in 2017, and 15-21 days in 2018. Then, to reduce the time discrepancy with SAR data, the cloud-free pixels were interpolated to daily data, by using a Kernel Ridge Regression machine learning method [41]. Finally, time series from ground sensors were filtered, selecting only the acquisition between 12 and 4 pm. The values within this interval of time were averaged in order to reduce the noise.

Snow-Melting Detection
To detect temporal changes in snow structural and dielectric properties due to the presence of water with S-1 C band backscatter, we used the ratio method [42]. Starting from the climate database of the Norwegian Meteorological Institute (https://klimaservicesenter. no/ Norwegian Centre for Climate Services (NCCS), accessed on 5 August 2019,) we selected winter intervals with specific characteristics of snow depth, air temperature and time between extreme events in days (snow depth ≥ 5 cm, mean and max of temperature ≤ 0 • C, ±5 days before and after). The same procedure was applied to the summer period (snow depth = 0 cm, mean and min of temperature ≥ 5 • C, ±5 days before and after, no precipitation). The meteorological station is located at 15 m a.s.l in the central part of Adventdalen (latitude 78.2022 • -longitude 15.831 • ). Using this information, it was possible to select images during the winter as a reference for dry snow/snow covered, as well during the summer period as snow free references. The specification of the selected images are shown in Table 2. As a first step, the 2017 time series was used to optimise a wet-snow threshold (Th). Three images acquired during the winter period were averaged, representing the snowcovered conditions (σ sc 2017 ).
The same procedure was applied to three summer images, in order to obtain a reference image for the snow-free surface (σ sf 2017 ). The threshold for defining wet-snow surfaces was quantified by calculating the ratio (R Th ) of these two images (Equation (1)): Then, a sample of 30,000 snow-covered and snow-free pixels was extracted from the ratio, respectively. RGB S-2 images from the corresponding period were used as a guide to define snow-covered and snow-free surfaces. From the intersection of the frequency distribution of pixels values, a threshold was obtained. As shown in Figure 4, the intersection between the distribution of snow-free and snow-covered S-1 pixels was −2.8 dB.
This threshold was applied to discriminate wet-snow in the 2017 and 2018 time series, for a total of 115 wet-snow masks. Simultaneously, a stack of images was created between the winter reference images of the years 2017 and 2018 (σ sc2017 , σ sc2018 ) and each image of the time series, calculating their ratio according to Equation (2): where R TS (Equation (2)) is the ratio of each image of the time series (σ time series ) and the reference images ( σ sc reference ) of 2017 (σ sc2017 ) and 2018 (σ sc20178 ), respectively. After rationing each image with the reference one, we applied the previously fixed threshold. Finally, a binary classification of wet-snow (Mask wet-snow ) was produced for the years 2017 and 2018 (Equation (3)).

Snow-Free Mapping and Accuracy Assessment
The Normalized Difference Snow Index (NDSI) [12] was calculated from S-2 time series for the years 2017 and 2018 according to Equation (4) : where b03 (Equation (4)) corresponds to green band (central wavelength 0.560 µm) and b11 (Equation (4)) to shortwave infrared band (central wavelength 1.610 µm). A threshold of 0.6 was then applied to create a binary snow/snow-free map for the season 2017 and 2018 [36]. The threshold was chosen based on our experience from previous work on defining snow-free surfaces in Arctic environments [8].
The comparison between the wet-snow maps derived from Sentinel-1 and the snow maps derived from Sentinel-2 was carried out.
At the same time, the ground stations and the meteorological station n. SN99870 (Norwegian Meteorological Institute) were used to validate the Sentinel-2 masks. The validation of the optical snow cover/snow free mask was carried out using photos from the PhenoCams. When images were not available (ST1, ST3 season 2017), the NDVI from the ground stations was used [43] as: To obtain a more robust validation, air temperature at 2 m height and snow depth from the meteorological station (snow depth = 0 cm, mean and min of temperature > 0) were incorporated into the validation process.

Ground Sensor Data Analysis
The next phase focused on understanding the dominant components in the backscatter signal in relation to vegetation development. S-1 σ 0 time series were extracted in homogeneous areas of around 50 × 50 m corresponding to the ground stations. To remove the noise, a linear interpolation and a moving average filter [44] were applied to the time series.
Correlation analyses were performed between time series of S-1 and ground sensor data to understand the contribution of vegetation, and soil moisture on the SAR signal. Subsequently, in order to measure the importance of the factors in a quantitative basis, a dominance analysis [45] of the variables was carried out.

Results
First, a validation of the masks were carried out. Subsequently, the S-1 melting season was calculated for the year 2017 and 2018, together with the snow season from S-2. Finally, a sensitivity analysis between S-1 backscatter and ground information was performed.

Snow Masks, Inter-Satellite Cross-Comparison and Ground Validation
Applying the −2.8 dB threshold, wet snow masks were created for the 2017 and 2018 seasons. The masks were selected for the May-September period, when optical images were also available. Next we created a binary daily snow cover mask by applying the S-2 NDSI threshold. The S-1 wet-snow masks and S-2 snow masks were overlaid and visually verified together with RGB images. An example of S-1 and S-2 masks is illustrated in Figure 5; using a RGB image as the base (Figure 5a), the mask obtained from the NDSI (Figure 5b) and the mask obtained from the coefficients of backscatter (Figure 5c) are overlaid.
Then, the first snow-free day of the S-2 masks was compared with the data of the ground sensors. From the regression model, the coefficient of determination R 2 between S-2 masks and ground data was 0.73. The results are summarized in Table 3.

Snow Seasonality
The total number of pixels affected by both by the presence of snow and melting snow during the two seasons was comparable. However, a significant difference (quantified in terms of km 2 of wet-snow) was observed in the temporal distribution of the process.  Finally, to obtain an overview of the snowmelt and snow seasons for both years we created maps of 'wet snow', 'dry snow', and 'snow-free' areas. An example is shown in

Multi Sensor Analyses and Vegetation
To understand the evolution of the HV signal in respect of the NDVI, the Soil Water Content (SWC) and the temperature of soil, a temporal analysis was computed. Figure 9 shows an example of the S-1 time series related to ground information; the time series are from 2017, station ST3.
The Pearson correlation coefficient (R) and p-value between time series are shown in Table 4. In the 2017 season, except for the ST3 and ST6 areas, the R between HV ∼ NDVI and HV ∼ SWC was below 0.6, with p-values less than 0.05 (excluding ST3 and ST1). Conversely, the R between HV and temperature in the ST6, ST7, and ST9 stations were significant (R mean of 0.74). Through the Pearson correlation coefficient, the positive correlation between the time series was defined, except for the station ST7 (HV ∼ SWC). In 2018, the results of the R were significant for HV ∼ NDVI and HV ∼ SWC. Instead, the HV ∼ temperature correlation showed a low value in all five areas. Again, the ST7 station was negative correlated with HV and SWC, along with the HV and temperature of the ST9 area. To determine the order of factors soil temperature and moisture content, dominance analyses of the linear models were performed. The results in Table 5 illustrate the dominance, expressed as a percentage, of the NDVI and SWC variables in linear regression with the HV channel. However, in 2018, for the stations ST1 and ST3 the SWC influences the HV channel more than the vegetation.

Discussion
The threshold for detection of wet snow was calculated and applied to 115 S-1 C band images. To minimize the threshold dependency on the reference images, it was essential to average three images during the winter (σ sc 2017 ) and summer (σ sf 2017 ) periods. The threshold identified in this study was in line with the thresholds found in the literature [15,16,46]. However, the time lapse between two acquisitions was decisive in the definition of wet snow, because there were strong changes from one day to the next. For this reason, together with the low amount of snow, it was not possible to detect wet snow in some areas of the valley floor. A specific threshold analysis for vegetated areas, rocks, and sediments could improve our wet snow masks [17].
In terms of optical images, the validation of the NDSI masks by PhenoCams/NDVI was not possible for the 2018 season. Due to the early snow melting, the start of the RGB image acquisition occurred only in a snow-free period. For this reason, it was only possible to use the information from meteorological station. The NDSI snow-mask/field data correlation, considering both years, was found to have an R 2 of 0.73. This is significantly better than the results obtained by the optical sensor MODIS over northern Norway [8].
The two seasons investigated showed a variability in the melting period. The 2017 season experienced late melting, with peaks around the end of June. On the contrary, in 2018 the snow began to melt much earlier. On average, the discrepancy in the valley between the two years was around 20 days.
Once the masks were obtained and validated, it was possible to get an overview of the snow season 2017 and 2018. The SAR and optical time-series were fused into a multi-sensor and multi-temporal snow cover masks. At this latitude, due to low incident angle, after September it was not possible to detect the first day of snow covered with optical data in the valley. The discrepancy between snow melting and the last day of snow corresponded in both cases to about 20 days at the lowest elevations in the two seasons under consideration. Therefore, in 2018, in these areas, snow melting was about 20 days earlier, as well as snow-free surface in comparison with the 2017 season. Considering that the vegetation phenology is strongly influenced by the dynamics of snow [2,3], it would be interesting to estimate the impact of the two seasons on a large scale [10], with a spatial resolution equal to 20 m (the spatial resolution of the S-2 short wave infrared band).
An exploratory investigation was conducted on the sensitivity of backscattering coefficients to snow-free surfaces in the high arctic environment. The structure of the vegetation and its dielectric properties, soil roughness and moisture influence the backscattering coefficients. For dense vegetation and small leaves, the C band has proved useful in detecting the dynamics of vegetation compared to other frequencies [47,48]. Moreover, in the cross polarized channels the contribution of vegetation is predominant [49]. Our results confirmed at most stations a dominance of the vegetation factor, expressed by the NDVI, on soil moisture (SWC). Nevertheless, the R between HV and NDVI showed, especially in 2017, a low significance and a variation depending on the areas under consideration. In the ST3 area, in both years, soil moisture was the prevalent component of σ 0 HV channel. The area has cryoturbated soil and is dominated by polygons of frost patterned ground, with little vegetation in the center. Soil moisture in this area averaged 0.30 m 3 /m 3 in both years, with maximum values not exceeding 0.34 m 3 /m 3 . Furthermore, the dominant vegetation of this area are graminoid that grow in tussocks; an example is the genus Luzula. These plants have short leaves (1-6 cm and 3-5 mm wide) and keep withered leaves and sheaths from previous years [50,51]. These factors could explain the dominance of the soil component over the vegetation in this area. On the contrary, station ST6 has a plain of sandy sediments dominated by herbs, such as the genus Festuca. In this case, the leaves are up to 10-20 cm long, with 0.7-1.0 mm broad when rolled up, up to 2.5 mm broad when flat [50,51]. For this reason, the sensitivity of the HV channel is greater with respect to the vegetation component than to the soil component. A similar argument applies to station ST1 and ST7, which have graminoid dominated cover. A further element in the analysis of the results of the station ST7 is the saturation of the soil. The SWC reaches maximum values of 0.6-0.8 m 3 /m 3 in 2017 and 2018, respectively. This could limit the penetration of the radar signal into the soil, resulting in limited sensitivity to this component. Another important point is the relevance of the correlations with temperature in the two years under consideration. In 2017, the year with delayed melting and disappearance of the snow, the SAR signal is more related to temperature at stations ST6, ST7, and ST9, than to vegetation and soil. On the contrary, in 2018, the relevance of temperature is almost negligible at all the stations. The change in snow cover determines the thawing rates of the soil, and thus controls the temperatures in the soil and soil surface, as well as the thickness of the active layer [52].
As the soil temperature decreases, there is lower liquid water content, which causes the backscattered signal to reduce [53,54]. The amount of frozen soil/tundra (with drier vegetation dominated by graminoids) may therefore limit the radar's capability to follow the vegetation growth pattern and soil moisture. In addition, links were observed between thawing of tundra soils and rainfall, especially from May to June [55]. For a phenological state estimation model, the incorporation of the soil temperature and precipitation could improve the outcomes. A better and deeper understanding of the backscatter dependence on soil and vegetation properties could be achieved by using data from boreholes in Svalbard (e.g., Adventdalen and Kapp Linne') in combination with fieldwork procedures such as those conducted by Bergstedt et al. (2018) [53] in their study covering circumpolar regions including Scandinavia. As a final consideration, the variation in the correlation between the vegetation and the backscatter signal could be determined by the presence of reindeer and migratory birds on the vegetated areas. The grazing of Svalbard reindeer (Rangifer tarandus platyrhynchus), whose population is increasing with an estimated mean population size of 22,435 [56], as well as grazing and grubbing by an exponentially increasing population of geese (e.g., Anser brachyrhynchus) cause a severe loss of plant biomass [57] that should be taken into account in understanding the response of backscatter coefficients to vegetation dynamics.

Conclusions
To follow snow dynamics on Svalbard, Sentinel-1 wet snow masks and Sentinel-2 snow masks were created and validated using ground data. An optimized threshold, applied to SAR backscattering coefficients, was used to detect wet snow during the years 2017 and 2018. Using the NDSI index, a daily snow mask was extracted from the Sentinel-2 sensor. Then Sentinel-1 and NDVI were used in synergy to follow vegetation dynamics. Our results have shown a variability in the length of snow seasons of 2017 and 2018. In 2018, the snow melted and the surface was free of snow about 20 days earlier than in 2017, and this was most clearly visible in the valley bottoms. In future studies, it would be interesting to evaluate the impact of these snow seasons on the vegetation phenology. Since high cloud coverage limits optical satellite data in the Arctic environment, a sensitivity study of cross-polarized channels to phenology and soil dynamics was performed.
The results of our study confirmed that vegetation is best detected by the HV channel in the Arctic environment. However, the ability to detect is limited due to the structure of the vegetation and the saturation of the soil. The late melting and disappearance of snow causes a further challenge in the monitoring of vegetation and soil dynamics. Therefore, when studying vegetation and soil dynamics, the amount of frozen soil/tundra should also be considered to understand the response of backscattering coefficients. Furthermore, the grazing of reindeer and arctic geese should be taken into account to fully understand the data. in the red-orange colours of the vallies. The difference between the maps is smaller for areas at higher elevations.