Sentinel-1 time series for mapping snow cover depletion and timing of snowmelt in Arctic periglacial environments: case study from Zackenberg and Kobbefjord, 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 depletion and snowmelt with high spatiotemporal resolution to capture their understudied small-scale heterogeneity. We use 97 dual-polarized S-1 SAR images acquired over northeastern Greenland and 94 over southwestern Greenland in the interferometric wide swath mode from the years 2017 and 2018. Comparison of S-1 intensity against SC fraction maps derived from orthorectified terrestrial time-lapse imagery indicates that SAR backscatter can increase before a decrease in SC fraction is observed. Hence, the increase in backscatter is related to changing snowpack properties during the runoff phase as well as decreasing SC fraction. We here present a novel empirical approach based on the temporal evolution of the SAR signal to identify start of runoff (SOR), end of snow cover (EOS) and SC extent for each S-1 observation date during melt using backscatter thresholds as well as the derivative. Comparison of SC with orthorectified time-lapse imagery indicates that HV polarization outperforms HH when using a global threshold. The derivative avoids manual selection of thresholds and adapts to different environmental settings and seasonal conditions. With a global configuration (threshold: 4 dB; polarization: HV) as well as with the derivative, the overall accuracy of SC maps was in all cases above 75 % and in more than half of cases above 90 %. Based on the physical principle of SAR backscatter during snowmelt, our approach is expected to work well in other low-vegetation areas and, hence, could support large-scale SC monitoring at high spatiotemporal resolution (20 m, 6 d) with high accuracy.

1 Introduction

1.1 Snow in Arctic environments

Snow cover (SC) has been identified as an essential climate variable (GCOS-WMO, 2020) covering about 40 % to 50 % of the Northern Hemisphere during winter (Dietz et al., 2012; Rees, 2005; Tsai et al., 2019b). It plays an important role for various components of the Earth system like hydrology, ecology, and climatology as well as global energy, water and carbon cycles due to the large seasonal variability in snow extent on the Northern Hemisphere (between 4 to 46 million square kilometers in summer and winter, respectively Rees, 2005) and the specific physical properties (Arslan et al., 2017; Box et al., 2019; Dietz et al., 2012; Pedersen et al., 2018). Its high albedo (0.8 to 0.9) compared to snow-free coverages strongly influences the energy balance. Moreover, its insulating properties limit heat exchange between soil and atmosphere and, thereby, regulate the seasonal active layer thickness (Rees, 2005; 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, and the snowpack provides water and, with it, nutrients for the plants. Moreover, SC blocks the sunlight needed for photosynthetic activity. Further, the metabolic activity of plants is directly linked 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 topography that affects transport of snow due to wind and gravitational redistribution. Thereby, snow is shifted from windward slopes, ridges, and 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 Oguma, 2013; Kepski et al., 2017; Pedersen et al., 2018). Ablation is driven by temperature, turbulent fluxes and solar radiation (Mott et al., 2011, 2013), which, in combination with redistribution, leads to a high small-scale heterogeneity of SC during the 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, 2005). The decrease in SC duration and extent has been documented for the 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 a 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.

1.2 SAR remote sensing of snow

As such, SC monitoring requires remote sensing products at very high spatial and temporal resolution. Until recently, freely available synthetic aperture radar (SAR) data could not fulfill these requirements; however, as the Sentinel-1 (S-1) mission with its twin satellites provides freely available C-band SAR imagery at a spatial resolution of 10–20 m and with a 6 d repeat cycle, SAR data have become an attractive alternative (ESA, 2012; Marin et al., 2020). While ground-based methods lack in spatial coverage, established SC mapping methods based on optical spaceborne earth observation suffer from reduced temporal resolution due to cloud coverage and sun illumination, e.g., during polar night (Dong, 2018; Portenier et al., 2020; Tsai et al., 2019b). SAR remote sensing can overcome these limitations, as it operates independent of sun illumination and atmospheric conditions (Marin et al., 2020; Tsai et al., 2019b; Ulaby et al., 2014). Previous studies applying SAR data have indicated their ability to detect wet snow, as liquid water in the snowpack decreases the dielectric constant, leading to a higher absorption coefficient (Marin et al., 2020; Nagler and Rott, 2000; Ulaby et al., 2014). Hereby, the penetration depth decreases to a few centimeters, and most of the signal is absorbed or reflected by the uppermost parts of the snowpack, resulting in a strong decrease in the backscatter (Marin et al., 2020; Ulaby et al., 2014). This effect is the basis for the bi-temporal approach referred to as Nagler’s method (Nagler and Rott, 2000), which compares backscatter during snowmelt with reference satellite imagery acquired during snow-free or dry-snow conditions (Nagler et al., 2016, 2018; Nagler and Rott, 2000; Snapir et al., 2019). However, this approach is only able to detect wet snow. More advanced methods additionally use digital elevation models (DEMs) (Nagler and Rott, 2000; Storvold and Malnes, 2004; Thakur et al., 2017) or optical remote sensing data (Nagler et al., 2018; Snapir et al., 2019; Thakur et al., 2018) to monitor dry snow. Nevertheless, these approaches cannot cover the small-scale heterogeneity of SC, because DEM-assisted approaches overestimate dry snow in steep terrain, while the inclusion of optical remote sensing data decreases the spatiotemporal resolution of the product (Solberg et al., 2010; Storvold and Malnes, 2004; Tsai et al., 2019b).

Figure 1 Schematic illustration of the pixel-wise seasonal evolution of synthetic aperture radar (SAR) backscatter intensity above snow with its phases of melt and relevant snow properties, i.e., liquid water content, runoff and snow cover fraction. This illustration (modified from Tsai et al., 2020) is based on the findings of Marin et al. (2020) and Tsai et al. (2020) and shows the parameters used in this study: SOR – start of runoff; SOD – start of snow cover fraction decrease which varies in time (indicated by the range of SOD and the different snow cover fraction curves) depending on snow pack properties and melt rates; EOS${}_{\mathrm{cam}/\text{S-1}}$​​​​​​​ – end of snow cover detected by camera/Sentinel-1; EOD – end of snow cover fraction decrease. The understudied interaction between Sentinel-1 SAR backscatter and snow cover fraction during the runoff phase and the potential of SAR time series to detect snow cover depletion during melt are investigated by this study.

According to Tsai et al. (2019b), SAR time series have rarely been investigated for SC mapping as most studies apply their method only on a few scenes and mostly for one single year, hence not taking advantage of the abundant data available. Only Tsai et al. (2019a, c) developed a SC monitoring approach using the entire temporal information of a S-1 SAR time series by incorporating interferometric, polarimetric, and backscatter features as well as elevation and land cover information into a supervised classification approach. On the contrary, recent studies developed methods for snow depth estimation (Lievens et al., 2019) and snowmelt phase detection (Marin et al., 2020) based only on the distinct seasonal SAR backscatter signal of snow described in the following (Fig. 1): during winter, the liquid water content is close to 0 %, and the snowpack is nearly transparent to C-Band SAR, except for larger snowpacks (>1 m), where a sensitivity of C-band cross polarization is indicated by Lievens et al. (2019, 2022). The returned signal is, therefore, predominantly scattered from ground surface underneath the snowpack with backscatter similar to snow-free conditions (Marin et al., 2020; Nagler et al., 2016). With increasing temperatures snow is wetting, leading to a strong decrease in backscatter during the moistening and ripening phase (Marin et al., 2020; Nagler et al., 2016; Rees, 2005). In the last phase of snowmelt, which is called runoff, backscattering increases again until an intensity level comparable to snow-free ground conditions is reached (Marin et al., 2020). This steady increase is not fully understood yet, but three possible explanation are described by Marin et al. (2020): (i) the increasing surface roughness of the snowpack; (ii) an increased number and size of intrusions like ice lenses or snow grains; (iii) subpixel parts turning into a snow-free state, leading to patchy SC towards the end of the melting period causing a mixed pixel signal response (Marin et al., 2020).

We here propose a novel approach, which adapts Nagler’s method (Nagler and Rott, 2000) to the seasonal evolution of the SAR signal using thresholds based on the seasonal minima of the SAR time series as well as the backscatter derivative for fast, simple, but effective SC mapping during snowmelt. With this new methodology, we can derive start of runoff (SOR) (based on the method from Marin et al., 2020), end of snow cover (EOS), the extent of end-of-season SC and SC extent maps for each S-1 observation date during melt. We validate EOS and SC derived from S-1 time series (2017 and 2018) with a reference dataset, generated from time-lapse photography available for Zackenberg Valley (northeast Greenland) as well as Kobbefjord (also known as Kangerluarsunnguaq Fjord, southwest Greenland), and assess the interaction between SC fraction and backscatter intensity in detail.

2 Study site and datasets

2.1 Study area

The Zackenberg (ZRA) and Kobbefjord (KRA) research areas are part of the Greenland Ecosystem Monitoring program (GEM; https://www.g-e-m.dk, last access: 4 July 2021), which has performed long-term monitoring of ecosystem components since 1995 (ZRA) and 2007 (KRA). The ZRA is located in high-Arctic northeast Greenland (Fig. 2a) approx. 40 km west of the outer coast and 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 (Hansen et al., 2008; 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) (Hansen et al., 2008). Maximum snow depths vary also considerably from year to year (0.4–1.4 m) (López-Blanco et al., 2020). Glacial ice occurs mostly at higher altitudes (>1000 m) due to low precipitation (Mernild et al., 2007). The radiation budget is dominated by polar night and day, which have a length of 89 and 106 d, 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 m a.s.l. are dominated by small shrubs and grasses (Fig. 2e) (Elberling et al., 2008; Westergaard-Nielsen et al., 2017). With increasing altitude, the percentage of bare ground and rock increases, while above 600 m 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).

The KRA is located in low-Arctic southwest Greenland at the bottom of the 16 km long fjord Kangerluarsunnguaq Fjord/Kobbefjord, approx. 20 km southeast of Greenland’s capital, Nuuk (Fig. 2b). The KRA drainage area, which covers 32 km², is characterized by three valleys surrounded by steep mountains up to 1375 m (Fig. 2d). Three major lakes (5 % of total area) dominate the valley system (Abermann et al., 2021). The GEM ClimateBasis main weather station (KOB) is located at the eastern end of the largest lake and provides high-quality meteorological data. The mean annual air temperature in KRA at KOB is −0.1 ^∘C (2008–2020), with the highest monthly temperatures in July (10.6 ^∘C) and the lowest in February (−8.6 ^∘C). The mean annual precipitation is around 830 mm and about one-third of it falls as snow (Abermann et al., 2021). Both total annual precipitation and maximum snow depths vary considerably from year to year (470–1170 mm and 0.3–1.3 m, respectively), and the winters are often characterized by shorter warmer periods causing episodic snowmelt (Pedersen et al., 2015). The predominant winter wind directions (eastern sector) and the topography cause similar, yet highly heterogeneous, patterns of snow accumulation every year (Myreng et al., 2020). The transition from snow-covered to snow-free ground in KRA normally begins in late April or early May, and the lower elevated parts of the valley system will usually be snow-free around the beginning of June. The rapid snowmelt often concurs with strong discharge peaks in the monitored rivers (Abermann et al., 2019). The KRA valleys are dominated by dwarf shrub heath, dry south-facing slopes and smaller fen areas (Bay et al., 2008), whereas the steep mountain slopes and higher-elevated areas are dominated by bare ground, snow patches, rockslides and small hanging glaciers (Pedersen et al., 2015).

Physical snow properties (e.g., SC fraction using time-lapse cameras) are measured regularly within the GeoBasis program at both ZRA and KRA (Skov et al., 2019; Westergaard-Nielsen et al., 2017; Rasmussen et al., 2020) following generally the same protocols but being adapted to local conditions and logistics.

Figure 2 (a, b) Overview of the Zackenberg/Kobbefjord study area and Sentinel-1 scene located in northeastern/southwestern Greenland (see red dots in inlet of panel b). (c, d) Study sites of Zackenberg (c) and Kobbefjord (d) with camera field of view and location of time-lapse camera (green dot) and Zackenberg research station (yellow dot) in panel (c). (e) Master image of time-lapse camera with used field of view and location of Zackenberg research station (yellow dot). Background image from GeoBasis Zackenberg (Greenland Ecosystem Monitoring Secretariat, 2021). (f) Master image of Kobbefjord study area. Background image from GeoBasis Kobbefjord (Greenland Ecosystem Monitoring Secretariat, 2021).

2.2 Datasets

For this study, we investigated all S-1 single look complex (SLC) data interferometric wide swath mode (IW) acquisitions from the ascending orbit (relative orbit 74) over ZRA and one descending orbit (relative orbit 54) over KRA between the 1 January 2017 and 31 December 2018 from the European Space Agency (ESA) Copernicus Open Access Hub (ESA, 2021) and the Alaska Satellite Facility (Alaska Satellite Facility, 2021) at a repeat cycle of 6 d. As dual-polarimetric mode (HH + HV) started being operational for Zackenberg (orbit 74) only from 21 May 2017 onwards, all previous acquisitions were neglected, leading to a limited time series in 2017 (ESA, 2020). In total, 38 scenes in the year 2017 and 59 scenes in the year 2018 were downloaded for ZRA as well as 41 (2017) and 53 (2018) scenes for KRA. The local acquisition time is 18:30 UTC in the afternoon in ZRA and 6:45 in the morning (09:45 UTC) in KRA.

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 distant ranges and close foreslopes. It covers about 45 km² of the valley floor and the west-facing slopes of Aucellabjerg (see Fig. 2c, e). 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). In KRA, daily images cover two lakes, the valley floor and the surrounding slopes. The selected field of view excludes hidden slopes in far range as well as the lake surface and covers about 7 km². Images are acquired by a 6-megapixel HP Photosmart E427 camera system in JPG format (Rasmussen et al., 2020). We manually excluded images of bad quality due to cloud cover, fog or rain prior to the processing. In ZRA, 89 out of 137 images in 2017 and 91 out of 127 in 2018 were found useful (Fig. 3). In KRA, images before 1 May and after 30 September were also excluded due to mountain shadow, but start-of-season SC extent was visually cross-checked to ensure the correct SC area is covered. A total of 93 out of these 151 images in 2017 and 84 out of 153 images in 2018 were found useful (Fig. 3). The Greenland Ecosystem Monitoring program GeoBasis (Greenland Ecosystem Monitoring Secretariat, 2021) provides 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.

Further, we used the ArcticDEM with a spatial resolution of 2 m (Porter et al., 2018; Noh and Howat, 2015; Polar Geospatial Center, 2021) 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, 2021) satellites with a spatial resolution of 3 m are additionally used for the orthorectification process of the time-lapse imagery. We selected for both sites a single cloud-free PlanetScope acquisition during the snowmelt season of 2018 with many distinguishable snow patches and patterns for which a corresponding time-lapse camera image of good quality was available. The used scene in ZRA (satellite IDs: 1008 and 1014) was acquired on 12 August 2018 at about 11:00 local time. The scene for KRA (satellite ID: 0f2d) was acquired on 15 June 2018 13:00 UTC (10:00 local time).

Figure 3 Acquisition dates of available data: blue – Sentinel-1 (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. The developed approach also uses S-1 acquisitions before May and after October, but the acquisitions are not shown here for better visualization of the relevant observation period.

3 Methodology

3.1 Sentinel-1 snowmelt and snow cover products

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 (ZRA: IW2, bursts 3 to 5; KRA: IW3, bursts 6 to 7), calibration to backscatter coefficient β₀ and debursting. Thereafter, multi-looking (three looks in range, one look in azimuth), speckle filtering using a boxcar filter with a window size of 3 by 3 pixels and calibration to backscatter coefficient γ₀ using the terrain-flattening approach similar to Small (2011) (Small et al., 2021) with the ArcticDEM (Porter et al., 2018) was applied. The data were 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 γ₀ backscatter intensities in decibel (dB) for both HV and HH polarizations (Ullmann et al., 2019).

Our approach uses the γ₀ backscatter intensity time series and a simple set of thresholds as input to (i) identify day of year (DOY) of start of runoff (SOR) and end of snow cover (EOS), (ii) detect start-of-season snow-free areas and end-of-season snow-covered patches, and (iii) derive a SC extent map for each S-1 observation date during melt (Fig. 4). The threshold setting is based on the characteristic seasonal backscatter behavior above snow, which is described by Marin et al. (2020) and Lievens et al. (2019) as a strong decrease in intensity in early spring followed by an increase in intensity during melt in late spring and summer. To detect the seasonal development of intensities, we compute the seasonal pixel-wise backscatter minimum min(γ₀) during the melt period (1 March–31 August) as an adaptive reference. Thereafter, SOR_S-1 is determined (Eq. 1) as the day of year (DOY), where the backscatter intensity (γ₀) reaches the seasonal minimum min(γ₀) (i.e., in accordance with runoff phase detection of Marin et al., 2020):

$$\begin{array}{}\text{(1)}& {\mathrm{SOR}}_{\text{S-1}}(x,y)=\left\{\begin{array}{ll}\mathrm{DOY}& ,\text{if}{\mathit{\gamma }}_{\mathrm{0}}(x,y,\mathrm{DOY})=min\left({\mathit{\gamma }}_{\mathrm{0}}\right)\\ -\mathrm{1}& ,\text{otherwise}.\end{array}\right.\end{array}$$

Following SOR_S-1, EOS_S-1 is identified using two different approaches: the threshold-based approach or the derivative-based approach. First, we describe the former: EOS_S-1 is determined as the DOY, where the backscatter exceeds the seasonal minimum min(γ₀) 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 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 start-of-season snow-free if no distinct temporal signal is found and the backscatter does not exceed the threshold t for three consecutive acquisitions (Eq. 2).

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& {\mathrm{EOS}}_{\text{S-1}}(x,y)=\\ & \left\{\begin{array}{ll}\text{first DOY}& ,\text{where}{\mathit{\gamma }}_{\mathrm{0}}(x,y,\mathrm{DOY})\\ & >min\left({\mathit{\gamma }}_{\mathrm{0}}\right)+t\\ & \left(\text{three consecutive times}\right)\\ \text{start-of-season snow-free}& ,\text{otherwise}\end{array}\right.\end{array}\end{array}$$

Besides this threshold-based approach, we also implement an approach using the derivative of the backscatter time series to detect EOS. The derivative (Δγ₀_DOY(x,y)) is defined as the difference of backscatter in dB between two acquisitions (Eq. 3). We believe that the derivative will adapt automatically to different environmental and seasonal conditions and, hence, does not require an manual parametrization of the threshold. As the derivative is more prone to speckle-induced variations, we additionally filter the time series in the temporal domain with a window of four acquisitions (one before, two after). By using an equal number of images, we reduce effects of systematic backscatter variations between S-1A and B. We do not apply a spatial filter to keep the same high spatial resolution as in the threshold approach.

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{\Delta }{\mathit{\gamma }}_{\mathrm{0}}(x,y,\mathrm{DOY})& ={\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{filtered}}(x,y,\mathrm{DOY}+\mathrm{\Delta }t)\\ & -{\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{filtered}}(x,y,\mathrm{DOY})\end{array}\end{array}$$

Using the filtered time series, we first exclude areas that show less than 3 dB variation during the melt period (March–August) and classify them as start-of-season snow-free as we assume that such areas do not show a distinct backscatter behavior related to snowmelt. Then, EOS_S-1 is determined as the DOY, where the derivative is below 1 dB after the derivative has shown an increase larger than 1 dB (Eq. 4). If, thereafter, Δγ₀_DOY(x,y) exceeds 1 dB at a higher backscatter level than the previously detected EOS, EOS_S-1 is reassigned as DOY, where the derivative falls again below 1 dB.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathrm{EOS}}_{\text{S-1}}(x,y)=\\ & \left\{\begin{array}{ll}\text{start-of-season snow-free}& ,\text{variation below 3\hspace{0.17em}dB}\\ \text{first DOY}& ,\text{where}\mathrm{\Delta }{\mathit{\gamma }}_{\mathrm{0}}(x,y,\mathrm{DOY})<\mathrm{1}\mathrm{dB}\\ & \left(\text{after}\mathrm{\Delta }{\mathit{\gamma }}_{\mathrm{0}}\right(x,y,\mathrm{DOY})>\mathrm{1}\mathrm{dB})\\ \text{start-of-season snow-free}& ,\text{otherwise}\end{array}\right.\end{array}\end{array}$$

Due to the high-Arctic environment, end-of-season snow-covered patches can occur. However, the algorithms presented in Eqs. (2) and (4) 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 end-of-season snow-covered patches based on the following assumption: end-of-season snow cover has undergone several melt and refreeze cycles and has strong structural similarity to firn which contains large snow grains and ice lenses that act as targets causing strong backscattering, especially volume scattering (Marin et al., 2020; Nagler and Rott, 2000). Therefore, the HV backscatter in autumn and early winter, once the 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 confirms that this effect is occurring. Accordingly, pixels are classified as end-of-season snow-covered if their maximum HV backscattering max(γ_{0_HV}) in autumn (1 October–31 December) exceeds min(γ₀) by 9 dB and the detected EOS is after 15 August (i.e., DOY 227) (Eq. 5). The 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 but might not apply for other sites.

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \mathrm{EOS}\mathit{\_}{\mathrm{final}}_{\text{S-1}}(x,y)=\\ & \left\{\begin{array}{ll}\text{end-of-season snow-covered}& ,\text{if}\left(\begin{array}{c}max\left({\mathit{\gamma }}_{\mathrm{0}\mathit{\_}\mathrm{HV}}\right)>(min({\mathit{\gamma }}_{\mathrm{0}})+\mathrm{9})\text{and}\\ {\mathrm{EOS}}_{\text{(S-1)}}(x,y)>\mathrm{227}\end{array}\right)\\ {\mathrm{EOS}}_{\text{(S-1)}}(x,y)& ,\text{otherwise}\end{array}\right.\end{array}\end{array}$$

As a final step, the following Eq. (6) is used to derive SC maps from the EOS, the start-of-season snow-free and the end-of-season snow-covered products:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& \mathrm{SC}\mathrm{\_}{\mathrm{DOY}}_{\text{S-1}}(x,y)=\\ & \left\{\begin{array}{ll}\mathrm{0}\left(\text{no snow}\right)& ,\text{if}\mathrm{DOY}\ge \mathrm{EOS}\mathit{\_}{\mathrm{final}}_{\text{S-1}}(x,y)\\ & \text{or start-of-season snow-free}\\ \mathrm{1}\left(\text{snow-covered}\right)& ,\text{if DOY}<\mathrm{EOS}\mathit{\_}{\mathrm{final}}_{\text{S-1}}(x,y)\\ & \text{or end-of-season snow-covered}\end{array}\right.\end{array}.\end{array}$$

Figure 4 Workflow chart of proposed approach to detect (i) snow cover (SC) and timing of snowmelt as (ii) start of runoff (SOR) and (iii) end of snow cover (EOS) from Sentinel-1. EOS (iv) start-of-season snow-free and (v) end-of-season snow-covered areas are visualized together as the so-called melt layer (Figs. 11 and 12), whereas analysis is conducted separately.

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 could be derived from SOR/EOS and delineate where SC is contributing to runoff or not. For the further assessment of the S-1 SC products, this separation is neglected, as the terrestrial camera is not able to detect SOR.

3.2 Time-lapse imagery snow cover products

Orthorectified SC maps are generated from time-lapse camera images to validate the S-1 EOS and SC products (Fig. 5). (i) All time-lapse camera images are aligned to the master image (12 August 2018 – ZRA; 15 June 2018 – KRA) using discrete Fourier transformation (implemented in the Python package imreg_dft (Týč and Gohlke, 2014)). (ii) The stacked images are classified using the histogram minimum thresholding approach of Salvatori et al. (2011) in the following way: all pixels above and below the minimum of a bimodal histogram distribution are assigned as snow-covered and no snow, respectively. However, we use the brightness histogram instead of the blue band to avoid misclassification of lakes and ponds reflecting the sky (e.g., Westergaard-Nielsen et al., 2017). (iii) The master image is then orthorectified using the Python package georef_webcam (Buchelt, 2020) based on the approaches of Corripio (2004), Härer et al. (2016) and Portenier et al. (2020). The ArcticDEM (Porter et al., 2018) is used for the projection and ground control points (GCPs) of remarkable landscape features are derived from same-day high-resolution PlanetScope imagery to optimize the orthorectification procedure resulting in a projection with high geospatial accuracy. (iv) Thereafter, the projected master is coregistered to the PlanetScope image with additional GCPs. (v) Based on the georeferenced master, all aligned and classified images are projected to SC maps in the coordinate system of S-1 with a spatial resolution of 2.5 m and then aggregated to 20 m resolution by calculating the SC fraction (SC_Fraction_DOY).

Figure 5 Workflow chart of proposed approach to detect (i) snow cover (SC) and timing of snowmelt as (ii) end of snow cover (EOS) from time-lapse images. EOS (iii) start-of-season snow-free and (iv) end-of-season snow-covered areas are visualized together as so-called melt layer (Figs. 11 and 12), 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.

For comparison with the S-1 products, start-of-season snow-free areas, end-of-season snow-covered patches and EOS are derived from the SC fraction map using Eq. (7) (see also Fig. 5):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\mathrm{EOS}}_{\mathrm{camera}}(x,y)=\\ & \left\{\begin{array}{ll}\text{start-of-season snow-free}& ,\text{if}\mathrm{SC}\mathrm{\_}{\mathrm{Fraction}}_{\mathrm{DOY}\mathrm{\_}\mathrm{1}}<\mathrm{0.5}\\ \mathrm{DOY}& ,\text{where}\mathrm{SC}\mathrm{\_}{\mathrm{Fraction}}_{\mathrm{DOY}}<\mathrm{0.5}\\ & \left(\text{three consecutive times}\right)\\ \text{end-of-season snow-covered}& ,\text{if}\mathrm{SC}\mathrm{\_}\mathrm{Fraction}>\mathrm{0.5}\\ & \left(\text{for all acquisitions}\right)\end{array}\right.\end{array}.\end{array}$$

(i) Pixels which are less than 50 % snow-covered in the first observation in spring (SC_Fraction_{DOY_1}) are classified as start-of-season snow-free, whereas (ii) end-of-season snow-covered pixels 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, (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) is 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 only used for the backscatter–SC-fraction interaction analysis.

Figure 6 (a) Seasonal development of snow cover (SC) deducted from the orthorectified time-lapse imagery compared with Zackenberg in situ snow cover fraction data. Their level of correlation is represented by R² and RMSE. (b) Orthorectified master image and PlanetScope scene in the background with their respective ground control points (GCPs). The center of the yellow dot indicates the projected location of the GCP, while the red dot shows the true location of the GCP. The distance between both is the orthorectification error (RMSE = 9.4 m).

We tested the accuracy of the orthorectification procedure using 25 (26) GCPs to assure the correct location of the projected maps (Fig. 6b). The root mean square error (RMSE) of the orthorectification was 9.4 m (5.3 m) in Zackenberg (Kobbefjord), hence less than 1 S-1 pixel (20 m). Besides, SC fraction maps from the time-lapse imagery in Zackenberg are validated with not-orthorectified in situ data which used the same time-lapse images as source (Greenland Ecosystem Monitoring Secretariat, 2021; Skov et al., 2019; Westergaard-Nielsen et al., 2017) to assure that the alignment, classification and orthorectification of the time-lapse imagery do not induce errors in the seasonal development of SC (Fig. 6a). The correlation of in situ SC fraction and mean SC fraction of the orthorectified SC fraction maps is tested by calculating the coefficient of determination (R² = 0.98) and the RMSE (= 3.7 % SC fraction). With such high geolocation accuracy and high agreement in reproducing in situ SC fraction, we can assure minimized systematic errors in the generated SC fraction maps (SC_Fraction_DOY). The advantage compared to state-of-the-art validation procedures based on other satellite data is the high spatiotemporal resolution of the generated reference product (2.5 m with a temporal resolution of 1–10 d dependent on weather conditions). Thereby, we can evaluate the response of the SAR signal directly to changes in small-scale SC and the SC fraction.

3.3 Evaluation and assessment of the products

3.3.1 Backscatter–snow-cover-fraction interaction analysis

The development of S-1 backscatter intensity is compared to the SC fraction data derived from the time-lapse imagery (SC_Fraction_DOY). We use the SC fraction derived from the time-lapse imagery as an indicator to examine the timing 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 development during time-lapse imagery observation. Thus, areas with EOD after 15 August and start-of-season snow-free areas are excluded from analysis. Original S-1 γ₀ intensity data (in dB) as well as intensities, which were rescaled from 0 (melt season minimum min(γ₀) intensity) to 1 (mean snow-free summer intensity), are compared to same-day SC fraction data of the time-lapse imagery. Furthermore, we identified the temporal distance to SOD and EOD for areas with 100 % and 0 % SC fraction, respectively, in order to capture the development of backscatter intensities before SOD and after EOD. Negative days indicate the number of days before the first observable decrease in SC fraction (SOD). Positive days indicate the number of days after the observed snow cover fraction has reached 0 % (EOD).

3.3.2 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 of the threshold-based approach. 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 due to the different threshold basis (snow-free/dry-snow backscatter vs. seasonal minimum) and the use of different levels of preprocessed SAR data (σ⁰ vs. terrain-corrected γ⁰). Therefore, we investigated a threshold range from 2 to 8 dB for HV and from 2 to 10 dB for HH polarization. In addition, we carried out the same assessment for the derivative approach to compare the results with the best threshold approach. For each threshold configuration as well as the derivative, 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 analyses as our approach using the EOS to derive SC is not capable of detecting new snowfall events in autumn. The datasets of start-of-season snow-free areas and end-of-season snow-covered patches are included in the SC accuracy analysis as follows: the class start-of-season snow-free is converted to SC at the start of the season with start-of-season snow-free pixels assigned as no snow and not start-of-season snow-free pixels as snow-covered. The end-of-season snow-covered class is handled likewise as SC at the end of the season. 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, d) and derive the receiving operator characteristic (ROC) from the overall TP and FP rate. Further analysis is carried out for the derivative approach as well as for the global polarization–threshold configuration with the best ROC (highest TP rate + lowest FP rate). The threshold-based product assessment is considered weaker as it was conducted with the same reference used for identifying this configuration, whereas no a priori knowledge is required for the derivative-based products. (i) The seasonal development of the confusion matrix metrics is assessed. (ii) We compare the EOS layers based on the areas, where the S-1 product (EOS_final_S-1) and the georeferenced time-lapse imagery (EOS_camera) detect an EOS, by calculating R², RMSE and mean absolute error (MAE). (iii) We calculate the difference in EOS DOY between the S-1 product and the 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) d. The former gives the percentage of pixels, which were assigned to the temporally closest S-1 acquisition, and the latter two correspond to one or two S-1 revisit cycles, respectively. (iv) Additionally, mean and median of the difference dataset are calculated.

4 Results

4.1 Interaction between backscatter increase and snow cover fraction

We assess the distinct seasonal backscatter behavior of S-1 over snow in ZRA and KRA in comparison to the SC fraction maps based on the orthorectified time-lapse imagery (SC_Fraction_DOY). For ZRA, the backscatter intensity shows similar trends in both years and both polarizations. Before the SOD, we observe a longer period of very low backscatter values (20–30 d) in 2018, whereas this period is shorter (5–10 d) in 2017 (Fig. 7a). The seasonal minimum is reached in ZRA about 10 to 30 d before the SOD (Fig. 7a). A sharp increase in intensity within the last 10–15 d before SOD (Fig. 7a, b) 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. 7b). The remaining increase to 1 (mean snow-free summer backscatter intensity) occurs mostly during the decrease in SC fraction (Fig. 7b, c). The comparison of the absolute HH and HV intensities shows differences (Fig. 7a). With HH polarization, a higher overall intensity 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. 7a). The rescaled intensities are stable for both years and polarizations, and we observe a linear increase in rescaled backscatter with increasing SC fraction (Fig. 7c). In KRA, the period of very low backscatter values is very short and distinct (10 d) in 2017 with a sharp decrease beforehand but shows a prolonged period of reduced backscatter intensity for 2018 (Fig. 7e). The seasonal minimum is reached in KRA in less than 10 d before SOD (Fig. 7d). Almost no increase in intensity is observed within the last 10–15 d before SOD (Fig. 7d, e) for both years unlike in ZRA. The increase in the rescaled backscatter intensity shows again a linear behavior and occurs almost entirely during the decrease in SC fraction in 2017 and 2018 only that in 2018 the averaged backscatter value is not reaching values lower than 0.4 (Fig. 7e, f). The absolute backscatter values indicate generally a lower seasonal variability in KRA than in ZRA, especially for 2018 due to higher backscatter minima. The snow-free HV backscatter intensity in Kobbefjord is generally higher than in Zackenberg, whereas no such difference is visible in HH (Fig. 7d).

Figure 7 Development of Sentinel-1 backscatter intensity over the melting season compared to same-day snow cover (SC) fraction from orthorectified time-lapse imagery (SC_Fraction_DOY) for ZRA (a–c) and KRA (d–f). Negative days indicate the number of days before the first observable decrease in SC fraction (SOD – start of decrease). Positive days indicate the number of days after the observed snow cover fraction has reached 0 % (EOD – end of decrease). The x-axis range (grey area) between SOD and EOD in panels (a, d) and (b, e) is modified to the observed average period between these dates (10 d) and does not represent the actual time range. (a, d) Seasonal development of absolute backscatter intensities (γ₀) with standard deviation and (b, e) of rescaled intensities (0: seasonal backscatter minimum min(γ₀); 1: mean snow-free summer backscatter intensity). (c, f) Close view of the development of rescaled backscatter between SOD and EOD during the decrease in SC fraction.​​​​​​​

4.2 Selection of threshold, derivative and polarization for Sentinel-1 snow products

We assessed the seasonal SC development for each threshold t in combination with each polarization to identify the parameter configuration that best fits the seasonal SC development in the orthorectified time-lapse imagery (Figs. 8, 9 and 10). We further compared these results with the derivative-based approach. For the determination of t, two opposite developments influence the accuracy of the resulting datasets in ZRA: the higher t is, the more areas are mistakenly classified as start-of-season snow-free (Fig. 10a), 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. 10b). The constant increase in the 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. 8). In KRA a similar development for overestimation of start-of-season snow-free areas is detected (Fig. 10c), but the melt offset is generally moved towards positive values shifting the optimum towards lower t (Fig. 10d). Thereby, the highest SC accuracies are reached with lower thresholds independently of the used polarization. Comparing the two polarizations, we observe a higher negative offset in EOS detection for HH (Fig. 10b) but slightly lower overestimation of start-of-season snow-free areas (Fig. 10a) in Zackenberg, and we observe the opposite trend for Kobbefjord (Fig. 10c, d). The derivative approach adapts to differences in optimal thresholds (especially for ZRA in HH), resulting in a reduced overestimation of start-of-season snow-free areas for ZRA and KRA, while offset in EOS detection remains low (±3 d) for ZRA (Fig. 10b) and stable around +5 to +7 d for KRA (Fig. 10d). New SC due to snowfall events in autumn is not detectable by the proposed approaches (Figs. 8 and 9).

Figure 8 (a, b, d, e) Seasonal development of overall snow cover (SC) in Zackenberg of the Sentinel-1 products (SC_DOY_S-1) for the threshold–polarization configurations as well as the derivative (der) compared to orthorectified time-lapse imagery (SC_DOY_camera). Autumn SC is not detected by the proposed approach as only the seasonal decrease in SC is observable. (c, f) Receiver operator characteristic (ROC) derived from the overall true positive (TP) rates and false positive (FP) rates of the SC maps of (c) 2017 and (f) 2018 based on same-day S-1 acquisitions and orthorectified time-lapse imagery. The used threshold is depicted as a label next to each point.

Figure 9 (a, b, d, e) Seasonal development of overall snow cover and (c, f) receiver operator characteristic (ROC) of Kobbefjord.

Figure 10 The key snow cover (SC) mapping accuracy parameters influenced by the threshold and polarization setting are (a, c) overestimation of start-of-season snow-free areas for higher thresholds leading to an underestimation of SC in early melt season and (b, d) a negative offset in mean EOS for low thresholds leading to underestimation of SC in late melt season.​​​​​​​

As the ROC analysis in Fig. 8c and f and Fig. 9c and f indicates, the t values for the best results depend on the used polarization and the observed year and site. Higher TP rates in 2018 compared to 2017 are observed for both sites (Fig. 8c, f and Fig. 9c, f). In ZRA, 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. 8c, f). The increase in the optimal threshold is higher in HH than in HV (Fig. 8c, f). Hence in Zackenberg, the performance of a global threshold is more robust for HV. In Kobbefjord, the ROC shows no large differences between the years. The best accuracies are reached for both polarizations with a threshold of 2 to 3 dB. HV tends to produce higher TP rates than HH in Kobbefjord (Fig. 9c, f). In Zackenberg, the derivative approach works equally well for both polarizations, and the difference to the optimum threshold configuration is low. In Kobbefjord, the derivative approach using the HV polarization outperforms HH. The derivative generates lower accuracies in 2018, whereas in 2017 the reached accuracies are equal to the optimum threshold. We use the global threshold configuration that shows the best ROC response for both sites and years (HV – 4 dB) as well as the derivative based on HV for further analysis. The analysis using the threshold-based products is considered weaker than the derivative method, as a priori knowledge is required to identify the optimal global threshold. The resulting layer composites of EOS, start-of-season snow-free and end-of-season snow-covered, for this configuration are shown in Fig. 11 for ZRA and in Fig. 12 for KRA.

Figure 11 Melt layers (composites of day of year (DOY) for end of snow cover (EOS), start-of-season snow-free areas and end-of-season snow-covered patches) of Zackenberg: (a–b) reference product by the time-lapse imagery with 2.5 m resolution for (a) 2017 and (b) 2018; subset of the generated S-1 product with 20 m resolution bounded by the camera field of view with global parameter configuration (polarization: HV; threshold t: 4 dB) for (c) 2017 and (d) 2018; (e–f) the same area of the S-1 product with the derivative approach for (e) 2017 and (f) 2018.

Figure 12 Melt layers of Kobbefjord: (a–b) reference product by the time-lapse imagery reference; subset of the S-1 product (c–d) with global parameter configuration and (e–f) with the derivative approach.

The seasonal development of confusion matrix parameters (Fig. 13, two left columns) 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. In ZRA, FP and FN occur during early melt season in 2017, whereas in 2018 mostly FN during late melt is observed. The comparably high proportion of FP responses (14.6 %) in 2017 on DOY 177 (Fig. 13a, e) is caused by a late snowfall event in about the beginning of July undetected in the S-1 dataset. As we observe predominantly FN almost throughout all observations, SC is generally rather underestimated by the S-1 product in ZRA. An underestimation of early-season SC induced by an overestimation of start-of-season snow-free areas (Fig. 10a, e) is observed in 2017, whereas in 2018 the observed underestimation of late-season SC (Fig. 10b, f) as well as end-of-season snow-covered is caused by the temporal offset in EOS detection (Fig. 10b). This indicates that a global threshold of 4 dB might be slightly too high for 2017 and slightly too low for 2018 in ZRA, which is consistent with the ROC analysis (Fig. 8c, f). In KRA, predominantly FP is observed during the melting period, resulting in a slight overestimation of SC by the S-1 product (Fig. 13).

Figure 13 (a, b) Seasonal development of true positive, false positive, true negative, and false negative as well as overall accuracy of the snow cover maps of the year (a) 2017 and (b) 2018. Overall accuracy is always above 75 % and in more than half of the cases above 90 %. Generally, false negative responses or underestimation of snow cover is dominant. False positive responses or overestimation occurs mostly in 2017 around DOY 177. (c) Correlation between EOS detection from S-1 time series approach and time-lapse imagery. (d) Distribution of temporal offsets in EOS. Highlighted are used time ranges corresponding to 0.5, 1 and 2 Sentinel-1 revisits. (a–d) ZRA 4 dB threshold. (e–h) ZRA derivative. (i–l) KRA 4 dB threshold. (m–p) KRA derivative.

For the EOS product, a R² score of 0.41 (p<0.001), a RMSE of 13.5 d and a MAE of 9.4 d is observed in ZRA with the global threshold. Using the derivative increases the accuracy to an R² score of 0.63 (p<0.001), a RMSE of 11.6 d and a MAE of 8.1 d. The observed accuracy measures in KRA are low (R²<0) due to the limited variation of EOS in the camera field of view (Fig. 12). The regression density plots in Fig. 13 show the correlation, and the histogram plots in Fig. 13 show the range of temporal difference in days between the two datasets. Within ±2, ±5 and ±11 d, 15 % to 26 %, 46 % to 49 % and 72 % to 80 % of all pixel-wise EOS dates are detected for both sites and both approaches. Using the threshold, the mean (−3.1 d) and the median value (−4 d) of the histogram are slightly negative in ZRA, whereas with the derivative in ZRA, mean (0.0 d) and median (+1 d) are close to zero. In KRA, the threshold-based values of mean (+3.1 d) and median (−1 d) are slightly higher than in ZRA. Using the derivative approach shows, like in ZRA, a trend towards higher values (mean: +6.8 d; median: +3 d). This is consistent with the observed offsets in EOS detection (Fig. 10b, d) and the distribution of EOS shown in Figs. 11 and 12.

5 Discussion

5.1 Influence of SC fraction and snow properties on SAR backscatter intensity

We observe a distinct seasonal behavior in S-1 C-band backscatter with a clear decrease during early melt, reaching a minimum just prior to the onset of melting, followed by a constant linear increase towards the EOS. This is consistent with findings in Marin et al. (2020). According to our observations, this development, which is found in both polarizations (Fig. 7a), is probably driven by changes in the snowpack, e.g., increased surface roughness, larger size and number of snow grains like suggested by Marin et al. (2020), as well as by decreasing fractional SC. The prolonged period of low backscatter values and the higher variance in absolute backscatter intensity in ZRA compared to KRA is probably caused by overall extended melting periods in ZRA and possibly also influenced by the different acquisition times in the afternoon and morning, respectively. The afternoon acquisition usually detects wet snow earlier in the melt season due to the diurnal melt and refreeze cycle of the snowpack. In ZRA, about half of the increase in intensity occurs within 10–15 d before the decrease in SC fraction starts (Fig. 7a, b), whereas no such development is observed in KRA (Fig. 7d, e). The increase before SOD in HH is probably caused by the higher surface roughness, while increased depolarization at this rougher snow surface increases the backscatter intensity in HV. The absence of a pre-SOD increase in KRA could be either due to the different acquisition time in the morning but is more likely caused by systematic differences in SC distribution between ZRA and KRA: (i) snow depletion is occurring faster in KRA than in ZRA, leaving less time for alteration of the snowpack before depletion. (ii) Snow distribution shows higher heterogeneity on small scales resulting in smaller snow patches in KRA in comparison to ZRA. Hence, subpixel parts are turning snow-free almost directly after start of melt is observed. Another sign for this systematic difference indicating more homogeneous snow distribution on a pixel scale in ZRA compared to KRA is the temporal difference between SOD and EOD, which is about 10 d in ZRA but more than 20 d in KRA. The remaining intensity increase in ZRA and almost all the increase in KRA occur along with the decreasing SC fraction and are almost linear. It is potentially driven by the increasingly higher proportion of the signal coming from snow-free parts of the pixel and possibly a further alteration of the snowpack (Fig. 7c). 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 varying strength of 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 patchy wet snow in threshold-based SAR 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.

5.2 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 during early melt linked to an increased overestimation of start-of-season snow-free areas with higher t (Figs. 8, 9 and 10a, c); (ii) increased underestimation of SC during late melt linked to an increased offset of earlier EOS detection with lower t in ZRA (Figs. 8 and 10b). 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. 7a), season-dependent thresholds can produce accurate results in both polarizations but require in situ reference data. The better global performance of cross 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. 8c, f) could be caused by the limited length of the time series in ZRA and gaps in the time series in KRA. The increase towards higher values of the best fitting t for ZRA in 2018, which is in accordance with the observed higher seasonal backscatter difference (Fig. 7a), is possibly caused by higher overall snow depths observed in 2018 by López-Blanco et al. (2020). The slightly lower optimal thresholds (2 to 3 dB) in KRA might be caused by the lower variance in backscatter intensity due to the different acquisition time. The seasonal defined threshold can vary about ±1 dB around the optimum while still giving good results (Fig. 8c, f), which indicates that the used global threshold in HV for 4 dB is applicable. The derivative-based approach generates results similar to the optimal seasonal threshold. However, no optimization or further in situ reference data are required. The systematic positive offset in EOS detection for KRA shows a later detection of EOS by S-1 than compared to the camera (Fig. 10d). Hence, the derivative-based EOS might be sensitive to lower SC fraction than defined by EOS_camera (SC fraction < 50 %) and possibly also detects low fractional SC. This could also be the cause for the increased FP rates of the derivative approach compared to the optimal seasonal threshold (Fig. 8c, f and Fig. 9c, f). The lower TP rates of KRA in 2018 (Fig. 9f) might be caused by an increased overestimation of start-of-season snow-free areas.

The analysis using the threshold-based products is considered weaker than the derivative method, as a priori knowledge is required to identify the optimal global threshold. However, the degree of optimization for the threshold setting is reduced to a minimum with a global threshold instead of training a season- and site-dependent threshold. The approach using the derivative is less susceptible to this issue. Using HV and t=4 dB, EOS is detected with a reasonable accuracy within two S-1 observations (Fig. 13). Potentially a denser time series incorporating different S-1 orbits could improve the accuracy in snowmelt detection; however, different orbits need to be analyzed separately, due to differences in local incidence angle and, more importantly, different acquisition times, which need to be considered. However, our case study shows that our approach works well with both settings as different orbits (ascending in ZRA and descending in KRA) and, hence, different acquisition times (afternoon and morning, respectively) were used. The SC maps reproduce the overall SC development with 9 out of 10 SC maps above 80 % overall accuracy and more than half above 90 % (Fig. 13), while most of the error is due to underestimation of SC in ZRA and due to overestimation in KRA. Thereby, accuracies comparable to other latest SC mapping approaches using optical remote sensing with similar spatial resolution (e.g., Gascoin et al., 2019; Girona-Mata et al., 2019; Piazzi et al., 2019) are generated but with the advantage of being independent from cloud cover. 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 depletion of SC during melt is observable. Further, dense vegetation in other study areas might cause increased inaccuracies due to the insensitivity of C-band SAR for snow in such areas (Nagler et al., 2016; Tsai et al., 2019c). The differences between the threshold- and derivative-based approach are rather low, which indicates that both might be applicable on other sites. However, the derivative has the advantage of automatically adapting to different seasonal and environmental conditions. A systematic difference is observed in the detection of EOS, indicating a higher sensitivity of the derivative for low fractional SC.

5.3 Major advantages of the proposed approach compared to other recent SAR-based snow cover studies

With this new approach multiple advances are made compared to other recent studies on SAR-based SC detection and the current standard, Nagler’s method: (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) and omit 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 end-of-season SC and firn. Due to the simple backscatter threshold and derivative approaches, we keep the 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 similar low-vegetation areas are comparable to the ones observed here (Tsai et al., 2019c). The threshold setting must be assessed in more detail to confirm whether a global threshold is applicable also in other sites and years and to analyze which effects different snow properties, vegetation, substrate and local incidence angle might have on the seasonal backscatter behavior of S-1 above snow, e.g., seasonal minimum and variation, as well as the resulting product accuracy. The derivative adapts well to different environmental settings (high and low Arctic) and different seasonal conditions (high and low snow depths). Hence, the derivative approach could be applied to another site without requiring prior optimization of the threshold. Further, both approaches are based on the physical principle of SAR backscatter during snowmelt. Hence, they are expected to work well in other low-vegetation areas as the characteristic seasonal pattern has been observed in the Alps (Marin et al., 2020) as well as around the globe (Lievens et al., 2019). Parametrization, e.g., of the threshold or the melt season, however, must be adapted and adjusted to local conditions. Moreover, further research needs to address the transferability to areas with denser vegetation or human-induced activities which could influence the SAR backscatter signal.

With this new approach, many relevant parameters for SC monitoring are detected on a weekly basis by the here proposed approach: state and extent of SC during melt, end-of-season SC and start-of-season snow-free areas. Further, important hydrological measures like start of runoff (SOR) and end of snow cover (EOS) are 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, and Nagler’s method (Nagler and Rott, 2000) could be used to identify the start of wet-snow phase. Provided at a spatial resolution of 20 m, hydrological models could further use this information to derive additional parameters like snow water equivalent (based on reconstruction approaches presented, e.g., by Molotch and Margulis, 2008; Kerr et al., 2013; Bair et al., 2016; Rittger et al., 2016) or assess the delay of snowmelt runoff due to meltwater storage in the snowpack and the soil (Marin et al., 2020; Tsai et al., 2020). If available, different S-1 orbits could be used to increase the temporal resolution of the product but need to be analyzed separately due to differences in local incidence angle and acquisition times. Further, extra wide swath mode (EW) data could be used with our approach for snowmelt detection and snow cover depletion mapping on larger scales with coarser resolution. Thereby, our approach enhances monitoring of hydrological cascading effects and could support in combination with other methods and datasets a holistic hydrological monitoring of SC from the scale of a single catchment up to pan-Arctic observations.

6 Conclusions

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, the approach employs user-defined thresholds based on the seasonal backscatter minimum as well as the derivative of the time series to (i) identify start of runoff (SOR) and end of snow cover (EOS) as day of year (DOY), (ii) detect start-of-season snow-free areas and end-of-season snow-covered patches, and (iii) derive a SC extent map for each S-1 observation date during SC depletion. EOS and SC are compared to maps derived from aligned and orthorectified terrestrial time-lapse imagery providing much higher spatial (2.5 m) and higher temporal (1 to 10 d) resolution than the S-1 product.

We compared the seasonal evolution of the SAR backscatter intensity to orthorectified SC fraction maps based on same-day time-lapse imagery. We observe that in ZRA about half of the HH and HV backscatter intensity increase during snowmelt occurs within 10–15 d before the decrease in SC fraction starts, whereas no such increase was observed in KRA due to faster melt and, hence, less time available for snowpack alteration. 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 the decrease in 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 thresholds (2 to 8 dB), indicating the following major error sources: (i) underestimation of SC during early melt due to an overestimation of start-of-season 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 end-of-season SC due to a systematically earlier detection of EOS (increasing with lower thresholds) in ZRA; (iii) neglection of episodic and autumn SC due to snowfall events as only the depletion of SC is detectable. The variation in the optimum threshold is higher in HH, which causes HV to produce better results with a global threshold. Using a global threshold of 4 dB with the HV polarization, EOS is correctly assigned to the closest S-1 acquisition for 15 % to 27 % of the area, while 45 % to 49 % and 72 % to 80 % are correctly detected within a period of one and two S-1 repeat cycles (6 and 12 d). 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 %. Using the derivative instead produced similar results and adapts well to different environmental settings and seasonal conditions. Hence, the derivative approach could be applied to another site without requiring prior optimization of the threshold. Both approaches work well with different orbits and acquisition times. Being based on the physical principle of SAR backscatter during snowmelt and its characteristic seasonal pattern, the approach is expected to work well in other low-vegetation areas around the globe, but further research is required to confirm the transferability of the approach to other settings and especially also areas with denser vegetation or human-induced activities which could influence the SAR backscatter signal. Finally, a SC product with this spatiotemporal resolution (20 m – 6 d) is, to the best of our knowledge, not presented with any other open-data remote sensing approach.

Further improvement could take advantage of the combined use of different S-1 orbits to increase the temporal resolution of the product or incorporate the snow phase detection algorithm by Marin et al. (2020) and use Nagler’s method (Nagler and Rott, 2000) to detect the start of the wet-snow phase. Thereby, continuous S-1 snow monitoring could not only improve hydrological and climatological models but also lead to an enhanced understanding of the complex interactions between climate change, SC and the Arctic ecosystem.