Sentinel-1 time series for mapping snow cover and timing of snowmelt in Arctic periglacial environments: Case study from the Zackenberg Valley, Greenland

Snow cover (SC) and timing of snowmelt are key regulators of a wide range of Arctic ecosystem functions. Both are strongly influenced by the amplified arctic warming and essential variables to understand environmental changes and their dynamics. This study evaluates the potential of Sentinel-1 (S-1) synthetic aperture radar (SAR) time series for monitoring SC and snowmelt with high spatiotemporal resolution to capture their understudied small-scale heterogeneity. We use 97 dual-polarized S-1 SAR images acquired over north-eastern Greenland in the interferometric wide swath mode from the years 5 2017 and 2018. Comparison of S-1 intensity against SC fraction maps derived from orthorectified terrestrial time lapse imagery indicates an increase of the SAR intensity before a decrease of SC fraction is observed. Hence, increase of backscatter is related to changing snowpack properties during the runoff phase as well as decreasing SC fraction. We here present a novel approach using backscatter intensity thresholds to identify start and end of snowmelt (SOS and EOS), perennial snow and wet/dry SC based on the temporal evolution of the SAR signal. Comparison of SC with orthorectified time lapse imagery indicate that HV 10 polarization outperforms HH when using a global threshold. With a global configuration (Threshold: 4 dB; polarization: HV), the overall accuracy of SC maps was in all cases above 75 % and in more than half cases above 90 % enabling a large-scale SC monitoring at high spatiotemporal resolution (20 m, 6 days) with high accuracy.

layer thickness (Rees, 2006;Tsai et al., 2019b). The snowpack plays an important role for water storage and supply (Dietz et al., 2012;Marin et al., 2020) and it is a key factor for arctic phenology, ecology and the distribution of flora (Assmann et al., 2019;Ide and Oguma, 2013;Kepski et al., 2017). On the one hand, the thermal insulation protects plants from frost damages 25 and the snowpack provides water and, with it, nutrients for the plants. On the other hand, SC blocks the sunlight needed for photosynthetic activity. Further, plants are linking their metabolic activity directly to the timing of snowmelt (Assmann et al., 2019;Ide and Oguma, 2013;Kankaanpää et al., 2018;Kepski et al., 2017;Pedersen et al., 2018). However, the timing of snowmelt is highly variable in space and time and influenced by snow accumulation, redistribution, and ablation. The former two depend on the climatic conditions, e.g. latitudinal and altitudinal position and continentality, as well as on the local 30 topography that affects transport of snow due to wind and gravitational redistribution. Thereby, snow is shifted from windward slopes, ridges, steep and high terrain to wind sheltered leeward slopes, sinks and low-lying terrain (Elberling et al., 2008;Farinotti et al., 2010;Lehning et al., 2008;Mott et al., 2018;Pedersen et al., 2016). Therefore, the redistribution effects generate a similar pattern of snow accumulation at the end of the winter primarily driven by the local topography, while the overall amount of snow depends on the amount of solid winter precipitation (Buus-Hinkler et al., 2006;Farinotti et al., 2010;Ide and 35 Oguma, 2013; Kepski et al., 2017;Pedersen et al., 2018). Ablation is driven by temperature, turbulent fluxes and solar radiation (Mott et al., 2011(Mott et al., , 2013, which, in combination with redistribution, leads to a high small-scale heterogeneity of SC during melt season (Mott et al., 2018). Knowledge about SC is important because it is decreasing with rising temperatures, resulting in a negative self-strengthening feedback between temperature and SC, which might partially drive the arctic amplification (Hock et al., 2019;Meredith et al., 2019;Rees, 2006). The decrease of SC duration and extent has been documented for the 40 northern hemisphere over the last 40 years (Box et al., 2019;Brown and Robinson, 2011;Meredith et al., 2019)) and could have a negative impact on species richness (Niittynen et al., 2018), but changes on local scale do not indicate clear trends (e.g. Pedersen et al., 2016;Young et al., 2018). Moreover, Hock et al. (2019) state that knowledge about SC distribution is still limited, especially at small spatiotemporal scales. 2 Study site and datasets

Study area
The Zackenberg research area (ZRA) is located in Northeast Greenland (Fig. 2a) approx. 40 km west of the outer coast and 90 70 km east of the inland ice sheet (Meltofte and Rasch, 2008). ZRA covers most parts of the Zackenberg valley floor and the surrounding slopes of Zackenberg Mountain (west) and Aucellabjerg (north and east) (Fig. 2c). The mean annual temperature in the valley is about -9°C and mean daily temperatures are usually above 0°C from early June until mid-September Pedersen et al., 2018). The mean annual precipitation, with 80 to 85 % falling as snow, is about 260 mm but varies largely from year to year (150 -400 mm) . Glacial ice occurs mostly at higher altitudes (>1000m) due 95 to low precipitation (Mernild et al., 2007). The radiation budget is dominated by polar night and day, which have a length of 89 days and 106 days, respectively (Pedersen et al., 2016). The climate is rather continental with high temperature fluctuations and low humidity due to the build-up of sea ice during winter in the Young Sound (Westergaard-Nielsen et al., 2017). The topographic setting and predominant northern winds during winter cause similar patterns of snow accumulation every year (Elberling et al., 2008;Hinkler et al., 2008;Kankaanpää et al., 2018;Pedersen et al., 2016). Vegetation formations below 200 100 m a.s.l. are dominated by small shrubs and grasses ( Fig. 2b) (Elberling et al., 2008;Westergaard-Nielsen et al., 2017). With increasing altitude, the percentage of bare ground and rock increases, while above 600 a.s.l. only scarce vegetation is found (Buus-Hinkler et al., 2006;Elberling et al., 2008). The transition from SC to snow-free ground in ZRA begins in late May in years with low snow accumulation, but can be prolonged until early July in years with high snow accumulation and even until late July in the very extreme year of 2018 (López-Blanco et al., 2020). A peak in photosynthetic activity is usually reached in late July to early August (Meltofte and Rasch, 2008;Westergaard-Nielsen et al., 2017) and the vegetation period usually ends in early September. Physical snow properties (e.g. SC fraction using time-lapse cameras) are measured regularly within the GeoBasis program (Skov et al., 2019;Westergaard-Nielsen et al., 2017).

Datasets
For this study, we investigated all S-1 single look complex ( Figure 3. Acquisition dates of available data: blue: S-1; grey: camera time lapse imagery excluded due to bad quality; orange: used camera images; green: coinciding S-1 and camera dates, which were used for comparison and assessment of the products.
In ZRA, daily time lapse imagery of the central part of the Zackenberg Valley has been taken from the east facing slope of the Zackenberg mountain at solar noon since 1997 (Skov et al., 2019;Buus-Hinkler et al., 2006). The selected camera field of view excludes far distant ranges and close foreslopes. It covers about 45 km2 of the valley floor and the west-facing slopes of Aucellabjerg (see Fig. 2cd). The time lapse imagery from 2017 and 2018 was taken with a 10 Megapixel Canon EOS 1000D camera system and stored in a 24-bit JPG format (Skov et al., 2019;Westergaard-Nielsen et al., 2017). We manually excluded 120 images of bad quality due to cloud cover, fog, or rain prior to the processing. In total 89 out of 137 images in 2017 and 91 out of 127 in 2018 were found useful (Fig. 3). The Greenland Ecosystem Monitoring program GeoBasis (Greenland Ecosystem Monitoring Secretariat, 2020) provide measures on SC fraction in ZRA derived from the above-mentioned camera (Skov et al., 2019;Westergaard-Nielsen et al., 2017), which is used as additional validation of the reference SC maps generated for this study.

125
Further, we used the ArcticDEM with a spatial resolution of 2 m (Porter et al., 2018;Noh and Howat, 2015;Polar Geospatial Center, 2020) for the processing of the S-1 scenes and the camera images, as it provides sufficiently high spatial resolution and accuracy needed for the calibration and orthorectification (Candela et al., 2017;Meddens et al., 2017). High resolution optical data from PlanetScope (Planet Team, 2019) satellites with a spatial resolution of 3 m are additionally used for the orthorectification process of the time lapse imagery. We selected a single cloud-free PlanetScope acquisition during the snowmelt season  The processing of S-1 products followed the workflow developed in Ullmann et al. (2019): The scenes were processed in IDL (version 8) and ENVI (version 5) using the functionalities of the ESA SNAP software (version 7). The processing included the application of the most recent orbit-file, split, calibration to backscatter coefficient β0 and debursting. Thereafter, multilooking (1 in range, 3 in azimuth), speckle filtering using a boxcar filter with a window size of 3 by 3 pixels and calibration to applied. The data were then terrain corrected using the Range-Doppler approach (Richards, 2009;Ulaby et al., 2014) and the ArcticDEM (Porter et al., 2018). Areas of shadow and layover were masked out. Datasets were processed to a spatial resolution of 20 by 20 m and all scenes were stacked and resampled to a common grid using bi-linear interpolation. The final product of S-1 preprocessing is the temporal stack of γ 0 backscatter intensities in decibel (dB) for both, HV and HH polarizations (Ullmann et al., 2019).
Our approach uses the γ 0 backscatter intensity time series as input and a simple set of thresholds, which are based on the characteristic seasonal backscatter behaviour above snow, to detect not only wet SC, but also ( Following SOS S-1 , EOS S-1 is determined as the DOY, where the backscatter exceeds the seasonal minimum min(γ 0 ) by more than a user-defined threshold t (Eq. (2)). Two further conditions are applied to filter events of melting followed by a freeze up of the snowpack, as these events cause a similar temporal evolution of the backscatter signal and, therefore, could lead to false detection: (i) In order to filter short-term snowpack refreeze events during the melt, the exceeding of t must apply for three consecutive acquisitions (Eq. (2)). (ii) If backscatter values with less than 2 dB difference to the seasonal minimum 160 are observed after a detected EOS and during early melt season (before 1 July), the pixel is considered as wet snow and a new EOS will be searched. Thereby, we filter early season short-term melt events followed by a refreeze of the snowpack. A pixel is classified as permanently snow-free, if no distinct temporal signal is found and the backscatter does not exceed the threshold t for three consecutive acquisitions (Eq. (2)).
permanently snow-f ree , otherwise (2) 165 Due to the high arctic environment in ZRA, perennial snow patches can occur. However, the algorithm presented in Eq. (2) to find EOS can only detect the end of the wet snow status, but cannot distinguish, whether the status changes to a snow-free or back to a dry-snow state. Therefore, we implement an additional condition to identify perennial snow patches based on the following assumption: Perennial snow has undergone several melt-and refreeze-cycles and has strong structural similarity to firn. Consequently, firn contains large snow grains and ice lenses that act as targets causing strong backscattering, especially 170 volume scattering (Marin et al., 2020;Nagler and Rott, 2000). Therefore, the HV backscattering intensity in autumn and early winter, once snowpack is completely refrozen, should be much higher than in areas with snow completely melted and now covered by new snow with small grains. Visual assessment of the HV backscatter intensity confirms that this effect is occurring in the ZRA. Accordingly, pixels are classified as perennial snow, if their maximum HV backscattering max(γ 0_HV ) in autumn (1 October -31 December) exceeds min(γ 0 ) by 9 dB and the detected EOS is after 15 August (i.e. DOY 227) (Eq. (3)). The 175 threshold of 9 dB is set in accordance with our observations that HV snow-free summer intensity does not exceed the seasonal minimum by more than this value. and As a final step, the following Eq. (4) is used to derive SC maps from the EOS, the permanently snow-free and the perennial snow products: The resulting binary SC map at a specific DOY is 1 (snow-covered) or 0 (no-snow), if the DOY is before or after EOS, respectively. Optionally, further information on the snow status can be derived from the SOS. The SC is dry or wet, if the selected DOY is before or after SOS, respectively. For the further assessment of the S-1 SC products, the separation of wet and dry snow is neglected, as the terrestrial camera is not able to detect SOS.  Validation with Ground Control Points DOY camera camera S1 S1 S1 S1 B a c k s c a tt e r -S C fr a c ti o n in te ra c ti o n a n a ly s is EOS, permanently snow-free areas and perennial snow are visualized together as so-called melt layer ( Fig. 9), whereas analysis is conducted separately. SOD/EOD: start/end of snow cover decrease, i.e. last DOY, where SC fraction is at 100 % and first DOY, where SC fraction reaches 0 %, respectively. SOD and EOD are used for the backscatter -SC fraction interaction analysis.

SOD & EOD
For comparison with the S-1 products, permanently snow-free areas, perennial snow patches and EOS are derived from the 200 SC fraction map using Eq. (5) (see also Fig. 4): (i) Pixels, which are less than 50 % snow covered in the first observation in spring (SC_F raction DOY _1 ), are classified as permanently snow-free, whereas (ii) perennial snow are at least 50 % snow covered in all orthorectified SC fraction maps. (iii) EOS is determined for the remaining area as the first of three consecutive dates, where SC fraction falls below 50 %. Besides,

205
(iv) Binary SC maps (SC_DOY camera ) for validating S-1 SC maps are generated using the boundary condition that pixels have to be more than 50 % snow covered in order to be classified as snow.
(v) Further, DOY of start and end of SC decrease (SOD/EOD) are defined for each pixel as the last observation with 100 % SC fraction and the first observation with 0 % SC fraction, respectively. SOD and EOD are used for the backscatter -SC fraction interaction analysis (see Fig. 4).
We tested the accuracy of the orthorectification procedure using 25 GCPs to assure the correct location of the projected maps 210 (Fig. 5b). The root mean square error (RMSE) of the orthorectification was 9.4 m and, therefore, less than one S-1 pixel (20 m). The development of S-1 backscatter intensity is compared to the SC fraction data derived from the time lapse imagery (SC_F raction DOY ). We use the SC fraction derived from the time lapse imagery as an indicator to examine the timing 225 of backscatter increase. We select all S-1 pixels with sufficient spatial data coverage by the time lapse imagery and include only areas which cover the entire snowmelt evolution during time-lapse imagery observation. Thus, areas with SOD before 1 June and EOD after 15 August are excluded from analysis. Original S-1 γ 0 intensity data (in dB) as well as intensities, which were rescaled from 0 (melt season minimum min(γ 0 ) intensity) to 1 (mean snow-free summer intensity) are compared to same-day SC fraction data of the time lapse imagery. Areas with 100 % and 0 % SC fraction were further segmented according 230 to the temporal distance to SOD and EOD, respectively.

Assessment of product accuracy dependent on selected threshold and polarization
According to our observations, the selection of polarization and threshold t (Eq. (2)) is crucial for the accuracy of the S-1 snow products and using the standard threshold (2 to 3 dB in Nagler's method (Nagler and Rott, 2000;Nagler et al., 2016;Snapir et al., 2019)) might not be suited. Therefore, we investigated a threshold range from 2 to 8 dB for HV and from 3 to 10 dB 235 for HH polarization. For each configuration, the accuracy of the SC maps (SC_DOY S-1 ; SC_DOY camera ) was assessed in the following way: The assessment is based only on same-day S-1 acquisitions and orthorectified time lapse imagery to avoid errors due to temporal offsets (see used dates marked in green in Fig. 3). Acquisitions after 1 September are excluded from all analysis as our approach using the EOS to derive SC is not capable of detecting new snowfall events in autumn. The datasets of permanently snow-free areas and perennial snow patches are included in the SC accuracy analysis as follows: The class 240 permanently snow-free is converted to SC at the start of the season with permanently snow-free pixels assigned as no-snow, and not permanently snow-free pixels as snow-covered. The perennial snow class is handled likewise as SC at the end of the season with perennial areas assigned as snow-covered and not perennial areas as no-snow. We calculate confusion matrix metrics (true positive (TP), true negative (TN), false positive (FP) and false negative (FN) rates as well as overall accuracy) for these SC maps for the selected view-field of the camera (Fig. 2c) and derive the receiving operator characteristic (ROC) 245 from the overall TP and FP rate. Further analysis is carried out for the global polarization-threshold configuration with the best ROC (highest TP rate + lowest FP rate) and is conducted with the same reference used for identifying this configuration: (i) The temporal evolution of the confusion matrix metrics is assessed. (ii) We compare the EOS layers based on the areas, where the S-1 product (EOS_f inal S-1 ) and the georeferenced time lapse imagery (EOS camera ) detect an EOS, by calculating R 2 , RMSE and mean absolute error (MAE). (iii) We calculate the difference in EOS DOY between the S-1 product and the 250 orthorectified time lapse imagery for 2017 and 2018 together and examine the percentage of pixels covered at the same date as well as within different time ranges of less than 3 (±2), 6 (±5) and 12 (±11) days. The former gives the percentage of pixels, which were assigned to the temporally closest S-1 acquisition, the latter two correspond to one or two S-1 revisit cycles, respectively. (iv) Additionally, mean and median of the difference dataset are calculated.

Interaction between backscatter increase and snow cover fraction
We assess the distinct seasonal backscatter behaviour of S-1 over snow in ZRA in comparison to the SC fraction maps based on the orthorectified time lapse imagery (SC_F raction DOY ). The backscatter intensity shows similar trends in both years and both polarizations. Before the SOD, backscatter intensity is rather low, reaching a minimum about 10 to 30 days before the SOD (Fig. 6a). In 2018 we observe a longer period of very low backscatter values (20 -30 days), whereas this period is very 260 short (5 -10 days) in 2017. A sharp increase in intensity within the last 10 -15 days before SOD (Fig. 6ab) is observed for both years. Neglecting the logarithmic scaling of intensity in the rescaled backscatter intensity, we observe an increment from about +0.2 above the seasonal minimum to +0.5 to +0.6 above the seasonal minimum before SOD (Fig. 6b). The remaining increase to 1 (mean snow-free summer backscatter intensity) occurs mostly during the decrease of SC fraction (Fig. 6bc). The comparison of the absolute HH and HV intensities show differences (Fig. 6a). With HH polarization, a higher overall intensity 265 and a larger difference between seasonal minimum and mean snow-free intensity is observed compared to HV. This difference is lower in 2017 than in 2018 due to higher backscatter minima (Fig. 6a). The rescaled intensities are stable for both years and polarizations and we observe a linear increase of rescaled backscatter with increasing SC fraction (Fig. 6c). We assessed the SC evolution for each threshold t and polarization to identify the parameter configuration that best fits the SC evolution in the orthorectified time lapse imagery (Fig. 7). For the determination of t, two opposite developments influence the accuracy of the resulting datasets: The higher t, the more areas are mistakenly classified as permanently snow-free (Fig.   8a), as intensities of these locations do not exceed the defined threshold. In contrast, EOS is detected better using higher t, whereas low t values cause a negative offset and EOS is detected earlier than observed (Fig. 8b). The constant increase of the 275 intensity during the melt period causes this earlier detection of EOS with lower t values. Thereby, higher values of t result in an underestimation of SC in the early melt period, whereas lower values of t lead to an underestimation of SC during the later melt period (Fig. 7). New SC due to snowfall events in autumn is not detectable by the proposed approach (Fig. 7). results dependent on the used polarization and the observed year. The ROC reaches higher values in 2018 compared to 2017 (Fig. 7cf). The optimal value for t differs between the years with higher values found in 2018 (e.g. for HV from 3 to 4 dB in 2017 to 4 to 5 dB in 2018; Fig. 7cf). The increase of the optimal threshold is higher in HH than in HV (Fig. 7cf). Hence, the performance of a global threshold is more robust for HV. We used the global parameter configuration that shows the best ROC response for further analysis. This is reached with HV polarization and a threshold of 4 dB (Fig. 7cf). The resulting layer 285 composites of EOS, permanently snow-free and perennial snow for this configuration are shown in Fig. 9.
The temporal evolution of confusion matrix parameters (Fig. 10ab) shows that overall accuracy is always above 75 % and in more than half of the cases above 90 %. The lowest overall accuracy occurs during the melt period, when melt is the strongest and the decrease in SC is the highest. Whereas in 2017, FP and FN are occurring during early melt season, in 2018 mostly FN during late melt is observed. The comparably high proportion of FP responses (14.6 %) in 2017 on DOY 177 (Fig. 10a) is 290 caused by a late snowfall event about the beginning of July undetected in the S-1 dataset. As we observe predominantly FN almost throughout all observations, SC is overall rather underestimated by the S-1 product. An underestimation of early-season SC induced by an overestimation of permanently snow-free areas (Fig. 9ab, 10a) is observed in 2017, whereas in 2018 the observed underestimation of late-season SC (Fig. 10b) as well as perennial snow (Fig. 9cd) is caused by the temporal offset in EOS detection (Fig. 8b). This indicates that a global threshold of 4 dB might be slightly too high for 2017 and slightly too low 295 for 2018, which is consistent with the ROC analysis (Fig. 7cf). For the EOS product, a R 2 score of 0.41 (p<0.001), a RMSE of 13.5 days and a MAE of 9.4 days is observed. The regression density plot in Fig. 10c shows the correlation and the histogram plot in Fig. 10d shows the range of temporal difference in days between the two datasets. Within ±2 days, ±5 days and ±11 days, 21 %, 47 % and 72 % of all pixel-wise EOS dates are detected. The mean (-3.1 days) and the median value (-4 days) of the histogram are slightly negative, which is consistent with the observed negative offset in EOS detection (Fig. 8b).

Influence of SC fraction and snow properties on SAR backscatter intensity
We observe a distinct seasonal behaviour in S-1 C-band backscatter intensity with a clear decrease during early melt, reaching a minimum just prior to the onset of melting, followed by a constant increase towards the EOS. This is consistent with findings in Lievens et al. (2019) and Marin et al. (2020). According to our observations, this development, which is found in both 305 polarizations (Fig. 6a), is driven by changes of the snowpack, e.g. increased surface roughness, larger size and number of snow grains, as well as by decreasing fractional SC. About half of the increase in intensity occurs within 10-15 days before the decrease in SC fraction starts (Fig. 6ab). This increase before SOD in HH is probably caused by the higher surface roughness, while larger and denser snow grains likely cause increased volume scattering and an increase in HV intensity. The further increase in intensity along with decreasing SC fraction is almost linear and potentially driven by the increasingly higher 310 proportion of the signal coming from snow-free parts of the pixel and possibly a further alteration of the snowpack (Fig.   6c). This linear increase seems suited to derive SC fraction from SAR backscatter intensity (e.g. Luojus et al. (2006) and Koskinen et al. (2009)); however, such an approach would need to address the strong pre-SOD increase in intensity and effects from changing surface properties underneath the snow, speckle and the viewing geometry, which result in variability of the SAR signal and make the discrimination between snow-free areas and areas with low SC fraction in threshold-based SAR 315 approaches like Nagler's method (Nagler et al., 2016;Nagler and Rott, 2000), as well as in our approach, challenging. These results are of comparably high validity, as the observations have been compared to a high-resolution reference dataset of same-day SC fraction maps derived from orthorectified time lapse imagery.

Influence of threshold and polarization on the products
For the selection of thresholds, we observe two main drivers for SC mapping inaccuracies: (i) increased underestimation of SC 320 during early melt linked to an increased overestimation of permanently snow-free areas with higher t (Fig. 7, 8a); (ii) increased underestimation of SC during late melt linked to an increased offset of earlier EOS detection with lower t (Fig. 7, 8b). For an accurate result, these contrary effects need to be balanced. While using a season-independent global threshold leads to a better performance in HV compared to HH due to the lower absolute seasonal changes in backscatter intensities in HV (Fig.   6a), season-dependent thresholds can produce accurate results in both polarizations. The better global performance of cross-325 polarization for SC detection is in accordance with other studies applying Nagler's method (Nagler et al., 2016;Thakur et al., 2018), which also indicated better performance of the cross-polarization compared to the co-polarization channel. The lower ROC performance in 2017 compared to 2018 (Fig. 7cf) could be caused by the limited length of the time series or lower overall snow depths. The increase towards higher values of the best fitting t in 2018, which is in accordance with the observed higher seasonal backscatter difference (Fig. 6a), is probably caused by higher overall snow depths observed in 2018 (López-Blanco 330 et al., 2020). Additional in situ reference data could further improve the accuracy of the products by finding a season dependent optimal threshold. However, the threshold can vary about ±1 dB around the optimum while still giving good results (Fig. 7cf) which indicates that used global threshold in HV is applicable. Further analysis of the products with the global polarization-threshold configuration might be critical, as it is conducted with the same dataset used for identifying this configuration. However, the degree of optimization for the threshold setting is 335 reduced to a minimum with a global threshold instead of training a season dependent threshold. Using HV and t = 4 dB , EOS is detected with a reasonable accuracy within two S-1 observations (Fig. 10cd). Potentially a denser time series incorporating different S-1 orbits could improve the accuracy in snowmelt detection. The SC maps reproduce the overall SC evolution with nine out of ten SC maps above 80 % overall accuracy and more than half above 90 % (Fig. 10ab) while most of the error is due to underestimation of SC. Thereby, our product generates similar accuracies as other latest SC mapping approaches with 340 similar spatial resolution based on optical remote sensing (e.g. Gascoin et al., 2019;Girona-Mata et al., 2019;Piazzi et al., 2019). However, it has to be pointed out that autumn and episodic SC due to snowfall events are not detectable by the proposed approach as only the seasonal decrease of SC is observable and dense vegetation in other study areas might cause increased inaccuracies due to the insensitivity of C-band SAR for snow in forested areas (Nagler et al., 2016;Tsai et al., 2019c).

Major advantages of the proposed approach compared to other recent SAR based snow cover studies 345
With this new approach multiple advances compared to other recent studies on SAR-based SC detection and the current standard, Nagler's method, are made: (i) we use the entire time series instead of only a few images per year unlike most previous studies (according to Tsai et al., 2019b); (ii) we avoid the manual selection of reference images for Nagler's method (Nagler and Rott, 2000) omits challenges like finding a completely snow-free or dry snow scene as well as a potential deterioration of the reference caused by altered backscatter signals due to perennial snow and firn. Due to the simple backscatter threshold 350 approach, we keep analysis fast and reduce processing capacity compared to the supervised classification approach by Tsai et al. (2019a, c), who additionally calculated interferometric and polarimetric features. Even though the spatial resolution is lower in their approach (100 m), resulting overall accuracies for likewise low vegetated areas are similar to the observed ones here (Tsai et al., 2019c). Threshold setting must be assessed in more detail to confirm, whether a global threshold is applicable also in other sites and years. However, in order to include mapping of dry SC, no spatial or temporal deterioration of the product 355 is required by incorporating optical data (e.g. like in Nagler et al. (2018) and Snapir et al. (2019)). Instead, it is possible to map wet and dry snow exclusively from a SAR time series. Important hydrological measures like start and end of snowmelt (SOS / EOS) can additionally be derived, whereas the former is not detectable with optical remote sensing data. Potentially, the snow phase detection algorithm by Marin et al. (2020) could be incorporated to further separate melt phases in more detail.
Different S-1 orbits could be used to increase the temporal resolution of the product. With the additional possibility to map 360 perennial snow at decent accuracy and with a guaranteed 6-day observation range, many relevant parameters for SC monitoring are detected at a weekly basis by the here proposed approach: SOS, EOS, wet and dry SC, perennial snow and permanently snow-free areas. With this information available at a spatial resolution of 20 m, hydrological models could further use this information to derive additional parameters like snow water equivalent (e.g. based on the approach by Kerr et al. (2013)) or assess the delay of snowmelt runoff due to melt water storage in the snowpack and the soil (Marin et al., 2020;Tsai et al., 365 2020). Thereby, S-1 SAR time series enhances monitoring of hydrological cascading effects and could be used for a holistic hydrological monitoring of SC from the scale of a single catchment up to pan-Arctic observations.
In this study, we present a fast and simple approach for mapping snow cover (SC) and timing of snowmelt based on Sentinel-1 (S-1) synthetic aperture radar (SAR) time series. Using the distinct seasonal signal of backscatter intensity above snow, 370 the approach employs user-defined thresholds based on the seasonal backscatter minimum to (i) identify start and end of the snowmelt (SOS/EOS) as day-of-year (DOY), (ii) detect permanently snow-free areas and perennial snow patches and (iii) derive a SC map of wet/dry snow. EOS and SC are compared to maps derived from aligned and orthorectified terrestrial time lapse imagery providing much higher temporal (2.5 m) and spatial resolution (1 to 10 days) than the S-1 product.
We compared the seasonal evolution of the SAR backscatter intensity to orthorectified SC fraction maps based on same-day 375 time lapse imagery. We observe that about half of the HH and HV backscatter intensity increase during snowmelt occurs within 10-15 days before the decrease in SC fraction starts. From then onwards, backscatter increases linearly with decreasing SC fraction. Hence, changes in the snowpack (e.g. grain size and number, surface roughness) as well as decrease of fractional SC are drivers for the observed backscatter increase.
The new approach to map SC and snowmelt was tested with HH and HV polarizations and different backscatter intensity 380 thresholds (for HV from 2 to 8 dB; for HH from 3 to 10 dB) indicating the following major error sources: (i) underestimation of SC during early melt due to an overestimation of permanently snow-free areas caused by parts which do not exceed the selected backscatter threshold (increasing with higher thresholds); (ii) underestimation of SC during late melt as well as perennial snow due to a systematically earlier detection of EOS (increasing with lower thresholds); (iii) neglection of episodic and autumn SC due to snowfall events as only the seasonal decrease of SC is detectable. The variation in the optimum threshold is higher 385 in HH, which causes HV to produce better results with a global threshold. The optimal seasonal threshold value increases in accordance with snow depth. Using a global threshold of 4 dB with the HV polarization, EOS is correctly assigned to the closest S-1 acquisition for 21% of the area, while 47 % and 75 % are correctly detected within a period of one and two S-1 repeat cycles (6 and 12 days). The resulting SC maps are generated with an overall accuracy of always more than 75 % and in more than half of the cases above 90 %. A SC product with this spatiotemporal resolution (20 m -6 days) is, to our best