Identification of Seasonal Snow Phase Changes From C-Band SAR Time Series With Dynamic Thresholds

In cryospheric studies, the most critical variable is estimating the beginning and evolution of the seasonal snow accumulation and its melting process. We propose and evaluate a new method for identifying seasonal snow cover phase changes in the Argentinean Andes using time series of Sentinel-1 synthetic aperture radar data available in the Google Earth Engine platform. We used meteorological and optical Sentinel-2 data to validate snow presence. First, we investigated seasonal snow cover dynamics in different regions of interest (ROIs). We identified three land surface cover periods: bare soil, dry snow, and melting snow. This finding is significant because other studies show that bare soil and dry snow have similar backscattering responses in C-band. Our methodology uses time series derivatives and their positive and negative anomalies. Finally, we compare our results with those obtained with a fixed threshold change detection approach. Our method was able to detect the phase change between bare soil and dry snow period in 75% of the ROIs, while a fixed threshold of $-\text{2} \,\mathrm{dB}$ only detects it in 42% of the cases. Furthermore, the derivative method also detects in 92% of the ROIs time series the beginning of the melting period, showing that it is a promising methodology for operative systems.


I. INTRODUCTION
E VERY year, approximately one-third of the Earth's land surface may be covered seasonally by snow, which, together with sea ice, makes it one of the most dynamic components of the cryosphere [1]. The structure and dimensions of snowpacks are complex and highly variable, both in space and time. Hence, regular monitoring of snow cover (extent, depth, and water equivalent) is, according to the Global Climate Observing System (GCOS), an essential climate variable with a high priority for global climate monitoring [2].
Seasonal snowpacks in the Argentinean Andes are the primary source of surface runoff and water supply in the adjacent lowlands. They, therefore, represent a critical resource for local domestic consumption, irrigation, tourism, industries, and hydroelectric generation. In addition, the mountain snowpacks' socioeconomic, environmental, and recreational benefits are enjoyed by over ten million people yearly in the adjacency of the Andes and many tourists that visit the mountain range [3].
Most of the Andes' snowpacks locate in remote or inaccessible areas, which makes these snow-covered areas one of the relatively most minor studied mountain environments in the world [4]. Moreover, conventional in situ studies and measurements provide limited area information and may be dangerous, expensive, and time-consuming [5]. For the abovementioned reasons, considering the comprehensive areal coverage, temporal variability, inaccessibility, and remote location of many snow-covered regions, remote sensing is the ideal data acquisition technique for monitoring snow cover and its trends and developments on both spatial and temporal large scales [6].
The most important types of imagery for this application are optical (including the use of near-infrared), passive microwave, and also synthetic aperture radar (SAR) imagery since, in all of these regions of the electromagnetic spectrum, snow has distinctive properties that make it comparatively easy to detect [6].
Monitoring snow properties with SAR presents several advantages, including high spatial resolution, insensitivity to clouds, operability at night, and accurate distance measurements through interferometry and polarimetry analysis for object structure detection. On the other hand, the use of SAR for snow cover extent (SCE) presents limitations, especially in mountainous areas, such as shadowing, overlapping, and foreshortening, but also limitations regarding the difficulty of discriminating dry snow cover from the bare ground since dry snow is "transparent" to specific wavelengths such as X-, C-, and L-bands [7].
The literature reports that the backscattering coefficient of SAR signal (σ 0 ) from surfaces covered by dry snow is approximately the same as the σ 0 from bare soil of the same surface [5], [8]. When the snowpack's liquid water content (LWC) increases in Spring, the dielectric constant increases, leading to a higher absorption coefficient that also decreases the signal received by the sensor [8], [9]. These changes enable the discrimination between wet snow from bare soil or dry snow, and it is the basis for the bitemporal approach proposed by Nagler and Rott [8]. Such bitemporal wet snow mapping uses two SAR images sensed at dry snow-covered (reference image) and wet snow-cover conditions (wet snow). It derives the wet SCE by thresholding the ratio image of their backscattering coefficients. This method has been widely used to detect wet snow, and has also been applied in Argentina with ALOS PALSAR (Lband) [10] and in Sentinel-1 (S1) (C band) together COSMO-SkyMed (X-band) images [11].
However, using threshold methods when estimating snow parameters has limitations, such as selecting optimal reference images, as this is a critical factor. There is also the need to adjust the threshold values to better account for the variability of the SAR signal by surface cover type, surface roughness, soil moisture, and terrain variations. For example, the backscattered signal of an illuminated target area is highly dependent on the incident angle of the signal. The backscattered intensity is high at low incidence angles compared with that at higher incidence angles over the same illuminated area [12]. For C-band SAR sensors, incidence angles close to 45°are ideal for detecting very wet snow (>3% LWC), while it is less effective for lower values of LWC. Local incidence angles (LIAs) vary according to slope aspect, topography, SAR sensor looking angle, and orbit pass acquisition (ascending or descending) [13]. Lone et al. [14] described that differences greater than 3 dB in the backscattering signal may occur within a study area, depending on the LIA and whether or not the regions of interest (ROIs) are located in slopes oriented or not to the looking angle of the sensor. Lone et al. [14] classified the areas of study according to whether they are located in slopes oriented to the looking angle of the sensor (facing slopes), opposite to the looking angle (shadow slopes), and the perpendicular locations (perpendicular slopes).
Considering the number of parameters that affect the selection of a fixed threshold, Karbou et al. [15] examined different options to improve wet snow detection of a S1 SAR image time series in the French Alps. Other authors used a soft threshold to provide a probability map, instead of the binary dry-wet snow map [16]. This methodology was implemented for RADAR SAT images in Norway [17] and in the Alps [18]. In addition, Manickam and Barros [12] used LIA to determine the appropriate threshold values to detect wet snow based on LandSat-8 visible imagery and snow-covered area maps.
Another methodology seldom investigated for snow cover mapping is SAR time series analysis [7]. Time series analysis of wet and dry snow maps created with SAR images is fundamental in characterizing snowpacks since optical albedo changes with grain size, which, in turn, depends on wetness, and snow pollution, being the last issue of great scientific concern [19]. The successful launch into orbit of two of the S1 mission satellites made snow cover time series analyses possible, as in the Alps study by Notarnicola et al. [20]. In this framework, the temporal series analysis offers strong possibilities to perform change detection [21], [22].
According to Tsai et al. [23], it is possible to distinguish between dry snow, wet snow, and snow-free areas using supervised classification of an S1 SAR time series incorporating interferometric, polarimetric, and backscatter features, along with elevation and land cover information.
Dingman [24] proposed methods for automatic snowmelt phase detection (moistening, ripening, and runoff) with SAR time series considering two hydrological years in Italy, Germany, and Switzerland. The authors identify the start of the melting process by a decrease in the multitemporal SAR signal in the afternoon of 2 dB or more. Such a fixed threshold is the major shortcoming of these works.
In order to overcome the fixed threshold limitations, Buchelt et al. [25] proposed using thresholds based on the seasonal minima of snow SAR signal and the backscatter derivative to identify the start of runoff. The date for the end of the snow season (EOSS) is the day of the year when the derivative is below 1 dB after the derivative has shown an increase larger than 1 dB.
Google Earth Engine (GEE) [26] makes it possible to work with the S1 image collection with ground range detected (GRD) level dB products. This platform makes multitemporal datasets accessible and provides a powerful programming interface for processing and visualizing the information before its download.
Although the number of publications regarding the analysis of seasonal snow with SAR time series is increasing, there is scarce information about the Andes snowpacks. Most studies focus on phase detection during the melt snow period. We conducted a temporal SAR backscattering study in 2019 for 12 ROIs near Bariloche, Argentina, showing significant differences between bare soil and dry snow-covered soil [27]. The results suggested that applying a fixed threshold of −2 dB or −3 dB for snowmelt phase detection, most of the ROIs from Bariloche will be considered wet snow during the dry snow deposition period.
Since the dynamics of seasonal snow cover is a continuous process that manifests differently according to various factors, it is possible to analyze it through its derivatives. The function's derivative determines areas where the function is constant, increases or decreases, and the intensity of the changes. In the case of time series, the derivative considers the behavior before and after some point in "time" [28].
Extending the study presented by Beltramone et al. [27], [29], this study presents a characterization of the influence of aspect and slope in snow SAR signals and a novel method to identify phase change in seasonal snow cover time series.
This novel approach calculates the first derivatives of smoothed S1 SAR data on GEE to detect positive and negative anomalies times differential curves of σ 0 , associated with phase changes in seasonal snow cover. In addition, we compared our methodology with an automatic detection method based on a fixed threshold to assess the detection results.
The rest of this article is organized as follows. The study area, ROI, and satellite and ancillary data are described in Section II. The backscatter temporal series processing and analysis, together with the procedure for wet snow detection through applying a fixed threshold and first derivatives with the anomalies calculation, are presented in Section III. Section IV presents the results. In Section V, discussion, the results are analyzed and interpreted in the context of other studies. Finally, Section VI concludes the article by discussing the relevance and applicability of the methodology. Fig. 1 shows the study area, located near San Carlos de Bariloche city, on the South shore of the Nahuel Huapi lake, in the province of Río Negro, Argentina. It is one of the most important cities at the foothill of the Andes. The combination of its altitude, latitude, and predominance of west-northwest winds make a cool temperate climate with a dry season that presents a west-east precipitation gradient. The yellow circles in Fig. 1(a) are 12 ROIs in hills surrounding the city with seasonal snowpacks. Fig. 1(b) and (c) shows an S1 acquisition of the same study area as in Fig. 1(a), covered by dry snow and wet snow conditions, respectively. However, the snow information recorded in SAR imagery fundamentally differs from optical/multispectral imagery.

A. Study Area
SAR imagery records surface characteristics related to the roughness and dielectric properties. In contrast, optical imagery records the reflection/absorption of the incoming solar irradiation at the top layer of the surface [30].
As seen from the figures, since SAR-based snow detection and mapping is not as intuitive as it may be for optical images, analyzing changes in the dielectric properties throughout time is critical to characterize seasonal snow.

B. Regions of Interest
The ROIs selected in this work are the 12 ROIs identified by [27]. Irregular polygons represent these ROIs to best suit  Table II shows their characteristics according to the analysis based on an S1 image and the SRTM product stated in [27]. The average LIA vary between 26°and 76.6°. Three ROIs have LIAs higher than 70 • , which are values close to being discarded for being considered a SAR-image distortion. However, since the ROIs do not have values over 78 • nor fall in the masked-out areas defined at the layover and shadow mask, they have been included in the analysis [27]. We also considered the elevation, slope, and aspect in the analysis, whose average values are given in Table II. The mean range of elevations among the ROIs varies between 1690 tand 2061 m. The most considerable elevation difference between an ROI's maximum and minimum value is 324 m, corresponding to ROI number 12. ROIs 1, 3, 5, 7, 10, and 11 exceed 22 m elevation value range. ROIs 8, 10, 11, and 12 present mean slopes higher 30 • ; however, the range is larger in ROIs 1, 5, and 7.

C. Satellite Products
We obtained S1 backscattering time series from the COPER-NICUS/S1_GRD collection available at the GEE platform. The collection includes S1 GRD scenes in dB, which we processed using the S1 Toolbox [31] to generate a calibrated, orthocorrected product. We selected all S1 2019 descending pass (relative orbit 10), VV polarization, and interferometric wide (IW) swath mode images.

D. Ancillary Data
In order to validate the methodology proposed in this work, we used ancillary datasets and remotely sensed imagery.
We downloaded daily precipitation and temperature data from the GIOVANI-NASA web service. 1 In geographic coordinates, the area considered to request data is 71.94950 • W, 41.13830 • S, 71.78470 • W, and 40.92960 • S. The precipitation product GPM_IMERGDF contains daily values of accumulated precipitated data in mm for an averaged area. In contrast, the surface temperature product Daytime/Ascending, AIRS, contains daily values of skin surface temperature in K for an averaged area with 0.1 m × 0.1 m spatial resolution. The period considered was 1 January until 31 December, 2019.
We obtained Sentinel-2 Multispectral Instrument (Level-2A) images from the COPERNICUS/S2 SR collection available in GEE with less than 50% cloud cover in 2019.
We used precipitation and temperature time series and Sentinel-2 snow-covered maps to validate the phase changes in seasonally snow-covered areas.
We employed the Shuttle Radar Topography Mission (SRTM) digital elevation data (with 30 m spatial resolution), also available in GEE, to mask out areas with slopes steeper than 45 • . The product was also used to characterize the slope, elevation, and aspect, among other variables of the ROIs.
In order to avoid areas of radar-specific imaging geometric distortions, we created a layover and fort-shortening mask from the S1 image from 31 January, 2019, in the ESA-SNAP software [32].

A. Snow Fall Events Assessment
We used Sentinel-2 images as ground truth to detect snowfall events and satellite-derived meteorological data.
We applied a cloud mask based on the QA60 band to the images selected from the Sentinel-2 collection in GEE.
We analyzed optical Sentinel 2-MSI images, with a spatial resolution of 10 m, to identify dates on which surface cover started to change between different classes (bare soil, wet soil, dry snow, wet snow, and mixed classes). For this purpose, we calculated the normalized difference snow index (NDSI) This index represents the relationship of the reflectance subtraction between the green and the SWIR band. For Sentinel-2 MSI Instrument, band 3 corresponds to the green channel and 11 to the SWIR. In addition, band number 8, which represents the near infrared (NIR), was considered to differentiate water bodies from snow [33]. We applied the NDSI to all images obtained after filtering the GEE collection through the following criteria: date, tile, % percentage of cloud (in the scene), and area of interest (selected ROIs). After applying the indices, we assigned 1 to those pixels with NDSI > 0.4 and NIR > 0.1 to obtain a binary snow classification. 1 [Online]. Available: https://giovanni.gsfc.nasa.gov/giovanni For each Sentinel-2 acquisition, we counted the pixels classified as "snow" on the binary maps and calculated the snow cover percentage.
We counted liquid precipitation events for daily accumulated precipitation data greater than 5 mm. In addition, we recorded the number of snowfall events for simultaneous conditions in which precipitation data were abovementioned 5 mm, and temperature data were below 1 • C.
We used meteorological information, the Sentinel-2 snow cover percentage, and visual interpretation of the SAR time series to differentiate three surface cover conditions named as follows: bare soil, dry snow and melting snow. The rule for determining the snow periods using the NDSI optical index and land surface temperature is given as follows.
1) If the NDSI <0.4, then the surface cover type is bare soil.
2) If air temperature >1 • C, and the NDSI >0.4, then the surface cover type is melting snow. 3) If air temperature <1 • C, and the NDSI >0.4, then the surface cover type is dry snow.

B. Backscatter Temporal Series Processing
We carried out most of the processing and analysis of SAR data in GEE. The processing of S1 images and downloading of the spatial subset, the construction of time series by ROI, and the analysis of LIA, slope aspect histograms, and pixel numbers, among other results, were carried out in the GEE code editor interface. The data selected in GEE were exported for further analysis and interpretation with programming routines in QGIS [34] and R [35].
GEE implements the following preprocessing steps to derive the backscatter coefficient in each pixel: apply orbit file, GRD border noise removal, thermal noise removal, radiometric calibration, terrain correction (orthorectification), and conversion to dB. In this work, we despeckled all S1 images from 2019 with a Refined Lee Speckle filter using a 3×3 kernel. Once the images were processed, we created in GEE a time series for each of the selected ROIs [27]. The values in the time series represent the median backscattering for each date.
ROIs were characterized according to their relative position to the acquisition orbit and looking angle of the sensor considering what was stated by Lone et al. [14]. We used S1 images to calculate their LIA and backscattering range and the SRTM-DEM product to calculate the aspect and slope.
In order to characterize different land cover states, we analyzed sets of three dates representing bare soil, dry snow-covered soil, and melting snow. We used optical images to select three emblematic dates per surface cover period (soil, dry snow, and wet snow). Such a selection is critical in the presence of phase change between bare soil and snow-covered soil. We also referred to meteorological data to discard dates in which storms might have changed the soil dielectric constant. In addition, we calculated the median backscattering value for each surface cover period for each ROI. We selected the dates whose value for all ROIs was comparable to their mean backscattering response. Finally, we rejected dates previous or after surface cover phase changes. Then, the values from the pixels contained in each selected ROI for that date were analyzed by descriptive statistics, density curves, and boxplots. We used R software for this stage. The dates selected for bare soil, dry snow, and melting snow periods are shown in the line time of

C. Wet Snow Detection
In this work, we considered the proposal by Marin et al. [36], related to using a −2 dB fixed threshold based on the multitemporal SAR signal. We adapted that approach to the satellite data available for our study area. Since our study area in the Argentinean Patagonia does not have the exact temporal resolution of S1 ascending and descending images acquisition, and the time series considered in this work are those available from the GEE collection, the processing of images is different from the one proposed in [36].
Only S1 descending acquisitions were selected from the GEE collection. Then, time series were created according to the methodology mentioned in Section III-B and imported in R. For the multitemporal SAR signal analysis for wet snow detection, we considered a difference of −2 dB on two consecutive dates.
The 12-time series downloaded from GEE were imported in R and smoothed using a moving median of three consecutive points.
For the time series positive and negative anomalies approach, we calculated the first derivative on each ROI's smoothed time series. We used the standard deviation (std) of the derivative time series as a metric to determine positive and negative anomalies. Values greater than 1 std were considered positive anomalies, and values lower than 1 std were considered negative anomalies.

A. Snow Fall and Snow Cover Data
We analyzed daily precipitation and surface temperature time series to detect snowfall events. The blue circles in Fig. 3 show 19 The filtered Sentinel-2 collection presents 13 images with cloud cover less than 50% for the date range between 1 April, 2019, to 31 January, 2020. Fig. 4 presents the percentage of snow cover per ROI according to the methodology presented in Section III-A. On 21 April, all ROIs presented bare soil on the Sentinel-2 image analyzed. The subsequent acquisition without clouds is from 21 May. On this date, ROIs number 2, 3, 4, 5, and 11 present more than 60% snow-covered pixels. Only ROIs 9, 10, and 11 remain with bare soil land surface cover. By 30 June, all ROIs are fully covered by snow, except ROI 4, which presents a percentage of snow-covered pixels of 84%. According to the NDSI and binary maps calculated, the 2019 dry snow period reaches from early July to the end of September. On 28 September, the snow-covered percentage on the ROIs starts to decrease, but it remains higher than 80%. The December acquisitions show a bare soil land surface cover on the ROIs, except for ROIs number 4 and 5 which some pixels are still covered by snow. Their bare soil season starts in January 2020.

B. ROIs Backscattering Time Series
The resulting GEE collection of S1 images consists of 30 for 2019. Images of the same geometry acquisition for the study area occur every 12 days. As mentioned in Section III, the time series shown in Fig. 5 consist of each ROI's median σ 0 values for each date.
As it can be observed in Fig. 5, every ROI presents the highest backscattering values during Summer and Autumn in the southern hemisphere (December-June). This season represents the bare soil period described in the methodology. Then, there can be seen a decrease in the backscattered signal which corresponds with Winter season conditions (June-September), and the period when the ROIs are covered with dry snow. Finally, a period with the lowest backscattering value is registered, which represents  wet or melting snow in Spring (September-December). These observations agree with those reported by Beltramone et al. [27].
The black dotted vertical lines in Fig. 5 show the three visible breaks of the time series backscattering behavior (change from bare soil to dry snow-covered period, change of dry snow to wet snow-covered period, and the last phase change from wet snow to bare soil covered period again). The backscattering values for each of these seasons follow a decreasing order being greater for Summer, Autumn, and Winter than for Spring. It can also be seen that different ROIs group in different backscattering ranges. ROIs number 2, 5, 7, and 8 present higher values compared with the rest of the ROIs throughout the time series. On the other hand, ROIs number 1, 10, 11, and 12 present lower values of backscattering. Although the time series present a similar trend, each ROI presents different local maximums and minimums. For example, the red and blue circles indicate variations in the mean backscattering signal due to particular events such as precipitation or early snowfalls.
The boxplots in Fig. 6 display the summary of the three land surface cover types defined in the methodology for each ROI dataset described in Section III-B. Each subplot is an ROI in which the ordinates are the backscattering response in dB, and the three surface cover conditions appear in different colors (bare soil in pink, dry snow in orange, and wet snow in green). The ROIs' boxplots have been grouped according to their aspect location and organized in columns. Except for ROIs number 6, 9, and 12, it is possible to distinguish between the three land surface conditions since their median values are significantly different. Fig. 7 shows the smoothed histograms and boxplots of different surface conditions of ROIs 3, 5, and 6, respectively. ROI 6 was selected because the boxplot is as reported in the literature: the backscattering response of dry snow surfaces is similar to the response of bare soil. Furthermore, ROI 6 is located on a perpendicular facing slope to the acquisition orbit and looking direction of the sensor, so we chose ROIs 3 (also perpendicular) and 5 (facing-slope) to compare the seasonal evolution of the backscattering. The surface conditions represent the values of every pixel during three dates (see Section III-B) from each season (bare soil, dry snow, and wet snow). Table III shows the difference between consecutive dates of backscattering values for the 12 ROIs for the transition between bare soil and dry snow-covered period dates and dry snow to wet snow-covered period season. Five pairs of dates accounted for the heterogeneity of the ROIs' environments and snow evolution. The pair of dates were selected considering what is shown in Fig. 5. The beginning of the dry snow accumulation period is 31 May. The global minimum of the time series, indicating the beginning of runoff, is 28 September. In addition, there is a local minimum in the time series of some ROIs dated 24 June.

C. Fixed Threshold Method Applied to Wet Snow Detection
The results highlighted in red are the dates where a phase change is detected using a fixed −2 dB threshold. ROIs 9 and 12 detect the first change on 24 June, which coincides with the The backscattering difference from 24 June and 6 July in ROIs 9, 10, and 12 presents a positive difference of more than 2 dB as a consequence of the signal increase after the local minima.

D. First Derivative Anomalies Method Applied to Wet Snow Detection
Negative anomalies of the time series derivative represent a reduction in backscattering and represent two phase changes in the year. The first one is a change in the environmental conditions between bare soil and snow-covered soil. The second one is a change in the conditions between dry snow and wet snow in Spring. Positive anomalies represent the last phase of the wet snow period, the runoff, where the backscattering intensity increases until reaching a level comparable to bare soil conditions. Over the derivative curves, the vertical green and red lines are the negative and positive anomalies, respectively. In ROIs 3 and 5, the first negative anomaly corresponds to the beginning of snowfalls and the dry snow period on 31 May. The second negative anomaly is the beginning of the wet snow phase on 28 September. These dates coincide with those suggested for the phase changes, presented in Fig. 5.      Our derivatives method detects the phase change between bare soil and dry snow in 75% of the ROIs, while a fixed threshold of −2 dB only identifies 42% of those cases. Our method also detects the beginning of the wet snow period, corresponding to the second phase change, in 92% of the cases. The derivative method detected two phase changes in eight of the 12 ROIs, approximately 67% of the results. The fixed threshold method applied to detect the beginning of wet snow was only successful in 58% of the ROIs.

V. DISCUSSION
This section analyzes snowfall event assessment, snow cover percentage, and backscattering time series characterization. We present and discuss the results of applying different methods to   detect changes in surface cover (bare soil, dry snow, and wet snow) through time series analysis.

A. Snow Fall and Snow Cover Data Analysis
The methodology proposed for snowfall assessment indicates that there have been several snowfalls from April to October 2019. This information is also reflected in the snow-covered percentage analysis derived from Sentinel-2 images. With the available images, the snow-covered surface period ranges from 21 May to 27 November, as shown in Fig. 4. Between the last snowfall and the end of the snow-covered surface period, there is over a month for some ROIs. That depends on each ROI's particular conditions that allow the snowpack to remain on the surface longer. As it is well known in the literature, Sentinel-2-derived NDSI maps do not discriminate between dry or wet snow. With the time series presented in Fig. 4, we can only identify the beginning of the runoff period since there is a shift from full snow-covered ROIs to mixed land cover classes, as stated by Marin et al. [36]. However, a complementary analysis with SAR data is needed for detecting wet snow conditions [7]. In this study, we can discriminate the snow-covered period into dry and wet snow by analyzing the SAR time series presented in Fig. 5. On 28 September the backscattering time series reaches its minimum starting the runoff period, thus the wet snow or melting snow period.

B. Characterization of ROIs Backscattering Time Series
The snowfall events detected in Fig. 3 are also reflected in the ROIs time series presented in Fig. 5. During the bare soil cover period (from 7 January to 31 May), in ROIs number 1, 2, 3, 4, 6, 8, and 9, two local maximums can be seen on the dates 8 March and 13 April. Fig. 5 shows in red circles the backscattered signal from 8 March and in blue circles the backscattering from 13 April. For the mentioned ROIs, these dates show a higher backscattering value than the average bare soil. The first corresponds to a precipitation event of 40 mm. The second local maximum coincides with the first snowfall, also detected in Fig. 3. The increase in the signal can be attributed to the fast melting of the first snowfall due to the ground temperature that is usually not frozen at the beginning of Winter. It is possible to validate this information with the snow-covered time series presented in Fig. 4 where there is no snow accumulation on the 21 April acquisition. The second and successive snowfalls are registered from 29 May.
According to the NDSI time series analysis, the snow-covered period starts on 21 May until 27 November, comprehending the dry snow and wet snow period. The SAR time series for the 12 ROIs show lower backscattering values compared with the bare soil period, once the season changes from Autumn to Winter, and therefore, from bare soil to snow-covered land surface condition. According to several authors, a change in the backscattering signal, once the dry snow period begins, is expected due to the change in the dominant scattering process [8]. In contrast with bare soil, the dry snow scattering mechanism is the sum of all orders of multiple volume scattering, surface scattering, and volume-surface interaction. Because the dielectric contrast between dry snow and soil exceeds that between dry snow and air, the contribution of rough surface scattering arises primarily from the snow-soil rough interface and not from the air-snow interface (the air-snow interface may have a stronger scattering contribution when the snow is wet). The snow-soil interface contributes to radar observations by the rough soil surface scattering and in less proportion by attenuation through the snow [5], [37]. However, the backscattering change between bare soil and dry snow conditions is usually ignored since the values are similar. In this study, and according to to [27], the difference between bare soil and dry snow in some ROIs is larger than 2 dB, as they can be seen in Fig. 5. Therefore, even if the reduction of the backscattering occurs during the snow-covered period, it cannot be attributed to the accumulation of dry snow.
Snowpacks are complex multilayer structures. There is a large number of unknowns, upon which the SAR backscattering is dependent [36]. As it was established by Ulaby et al. [38], the scattering and attenuation processes are governed by two sets of parameters: sensor and target parameters. Sensor parameters, include wavelength (λ), angle of incidence (θ), (relative to normal incidence), and polarization configurations of the incident wave (or transmit antenna) (P t ) and receive antenna (P r ). Target parameters include the dielectric and geometrical properties of the air-snow interface, the snowpack characteristics [snow grain size (mm) and shape, stratigraphy, snow density (kg/m 3 ), and snow depth (cm)] presence of impurities, underneath vegetation], and the snow-ground interface.
A recent study presented the theoretical sensitivity to snow water equivalent in seasonal snow via volume scattering processes, including the influence of surface and ground contributions [37]. In the mentioned work, the full dense medium radiative transfer (DMRT) model was used to find regression formulas and solve different backscattering scenarios versus important geophysical variables. Soil surface scattering depends on the roughness and on the soil permittivity, which in turn depends on soil moisture (volume of water present per unit volume of soil) and texture (soil composition). The algorithm developed by Kim et al. [39] assumes that surface roughness remains constant over a time series period (wet and frozen soils are assumed to have the same rough surface and the surface rms height h) [39]. In this framework, Tsang et al. [37] found that backscattering at VV polarization for C-band simulations with an incident angle equal to 40°decreases almost two decibels from 10% moisture soil (rms height of 0.5 cm) to frozen soil moisture (rms height 0.5 cm) being circa −13.3 and −15.3 dB, respectively [37]. This difference between frozen and not frozen 10% moistured soil is comparable to some changes detected for ROIs 2, 3, 9, 10, and 12 during the dry snow deposition period, which could be associated with frozen soil events (see Table III, columns 3 and 4). In addition, ROI number 9 presents a similar local incident angle than the one used by Tsang et al. [37], 40.89°, and also shows similar backscattering values for 10% moisture soil and frozen soil which are around −13 and −15 dB, respectively. A rigorous evaluation of these assumptions will be done in future works by carrying out field campaigns to measure surface roughness, and soil moisture and to detect frozen soil events.
The SAR backscattering time series, of Fig. 5, shows a second reduction of values on 28 September (Spring in the southern hemisphere). On that day, the lowest value of the received signal is registered, starting the runoff phase according to [36], which also coincides with the increasing temperatures and liquid precipitations shown in Fig. 3. As mentioned by Marin et al. [36] and Dingman [24], the runoff phase starts when no more liquid water can be retained in the snowpack. After the time series σ 0 reaches its minimum, the backscattering signal increases again.
According to Snapir et al. [40], when analyzing the backscattering signal from different ROIs located in mountains, it is possible to distinguish between three groups of ROIs. The LIA influences that classification Lone et al. [14] categorized the three groups into "facing," "perpendicular," and "shadow," depending on their location relative to the acquisition orbit and looking direction of the sensor. Fig. 5 presents characterized in different colors the groups of ROIs classified according to [14]. The ROIs with higher backscattering are located on a facing slope and are colored in the range of pinks. The ROIs colored brown to orange are the ones located on perpendicular slopes. Moreover, finally, the ROIs with the lowest backscattering signal and colored blue are the shadow slope ones.
In this case study, the three groups of ROIs can be characterized as follows.
1) Facing slopes: The ROIs in this group present the highest σ 0 values. They are said to "face" the sensor's acquisition orbit and looking angle and correspond to numbers 2, 5, 7, and 8 with predominantly slope aspects ranging NE-E-SE. The LIA of these ROIs varies between 22°and 32°.
2) Shadow slopes: The ROIs in this group are said to "be hidden" from the sensor's acquisition orbit and looking angle and correspond to numbers 1, 10, 11, and 12. The predominant slope aspects are SW-W-NW, and the ROIs' LIA varies between 63°and 76°. These are the ROIs with the lowest σ 0 values throughout the series. 3) Perpendicular slopes: The ROIs in this group present shared attributes from both the facing and the shadow slope ROIs. They are ROIs numbers 3, 4, 6, and 9. According to the aspect range determined by the SRTM product, the slopes that face the sensor are the ones predominantly N and S. The LIA of the mentioned ROIs range between 35°a nd 57°. It is important to consider the slope aspect when characterizing the study area and selecting thresholds, as it was mentioned by Manickam and Barros [12] since there may be a difference of 7 dB between North-East and North-West slopes.
Further analysis of the ROIs backscattering, considering their relative location to the sensor acquisition orbit, consisted in identifying the median values of each land surface cover period. Fig. 6 shows that the breaks and changes of land surface cover are not only visible in the time series, but they also present significant differences in the median boxplots grouped by the surface land cover condition for each ROI. The boxplot diagrams from bare soil, dry snow, and wet snow period of ROIs 1,2,3,4,5,7,8,10, and 11 show the differences mentioned by Beltramone et al. [27]. This evidence suggests that the differences in the backscattering response among soil-covered surfaces do not depend on the aspect, elevation, or slope since they appear in facing-slope (2, 5, and 7), perpendicular (3 and 4), and shadow-slope ROIs (1, 10, and 11). ROIs 6 and 9 are exceptions to this rule: soil without snow and dry snow period are similar. In contrast, ROI 12 presents similarities between the median of the dry snow and wet snow season.
Through this analysis, we found evidence that before using the change detection methodology, one must consider that some areas of study may present significant differences between the bare soil and the dry snow period. Therefore, choosing a single image, combining bare soil images, or combining bare soil and dry show images as it was proposed by Teverovsky et al. [11] and Notarnicola et al. [20] to use them as a reference image may influence the wet detection map result for different ROIs.
Since the data show a significant difference between surface classes in the boxplot analysis for all cases, using thresholds for class discrimination will present poor discriminating potential.
When analyzing the smoothed histograms of different land surface covered periods, ROI 3 seasonal behavior shows a significant difference between the distribution of the backscattering values of bare soil, dry snow, and wet snow surface periods. Fig. 7(a) is an example of what is stated in [27], showing that in the study area, the backscattering coefficient of surfaces covered by dry snow is significantly different from the one of bare soil. For this reason, when applying the change detection method proposed by Nagler and Rott [8], it is impossible to use indistinctly or combined reference SAR images from the soil without snow and dry snow dates. In contrast, Fig. 7(c) shows the density curves and boxplots of ROI 6, which presents similar backscattering values for the bare soil and dry snow. ROIs 3 and are located at a similar altitude and in slopes oriented perpendicularly to the acquisition orbit and looking direction of the sensor. However, their positions differ: they are 18 km apart.
Finally, Fig. 7(b) shows the density curves and boxplots of ROI 5, where all the land surface cover types overlap and have different modes. For the same dates selected in the analysis of ROIs 3 and 6, in this case, it is not possible to differentiate surface cover types by analyzing the density curves.
In summary, the characteristics of backscattering for dry snow, wet snow, and bare soil period are highly different in the ROIs. The complex features of land surface types and their variability result in a suggestion to use dynamic thresholds instead of fixed thresholds.

C. Assessment of Fixed Threshold Versus First Derivatives Methodologies Applied to Snow Seasonal Evolution Time Series
When an adapted version of the algorithm for detecting melting phases proposed by Marin et al. [36] was applied to the 12 time series of the study area, five ROIs (numbers 2, 3, 9, 10, and 12) presented an intermediate phase change, corresponding to the beginning of winter season. This means that if a fixed −2 dB threshold is used in the study area, in some cases, there will be a first change of phase detection that corresponds to the transition from bare soil to snow-covered soil and, then, a second change of phase detection that corresponds to the transition from dry snow to wet snow.
What has been described previously for Table III shows that, in the area of study, the performance of methodologies for automatic snowmelt phase detection, as the one proposed by Marin et al. [36], is limited. At least, the thresholds should be adaptable to automatically detect snowmelt phase change in regions like Bariloche.
In addition, for the bitemporal approach proposed by Nagler and Rott. [8], Table III indicates that if a −2 dB threshold detects changes among two consecutive images either from bare soil or from dry snow surfaces, then either considering those surfaces indistinctly as reference images or combining them will result in inconsistent classification results.
Considering what was mentioned by Ventura et al. [41], the snow and snow-free ratio distributions often overlap. The choice of thresholds is a crucial aspect of the analysis. For this reason, we also considered different thresholds to overcome the detection of an intermediate phase change with this methodology. The thresholds were selected based on the ROIs density curves and listed by different authors for C-band SAR imagery [42], [43], [44]. However, different thresholds presented other limitations; for example, ROIs 1, 10, and 12 would not detect the melting phase with a −3 dB fixed threshold.
In order to overcome the difficulties presented by using a fixed threshold to analyze the time series in the study area, we developed a methodology for computing positive and negative anomalies with the standard deviation of the time series derivatives. We show that the time series derivatives provide more information about the underlying dynamics than the observed data, coinciding with what was presented by Buchelt et al. [25]. The time series dynamics are uncovered by their derivatives.
Another encouraging fact of applying our anomalies methodology for detecting two phase changes is that for ROI 5, the phase change between bare soil and dry snow-covered soil was not detected with Marin et al.'s [36] methodology. Therefore, for some ROIs where the algorithm for detecting wet snow was successful, the anomalies methodology also worked to detect two phase changes.
Only the beginning of wet snow is identified in ROI 6. This agrees with the results obtained from the density curves, boxplots, and existing literature.
For ROIs 3 and 5, the first and second positive anomalies correspond to the runoff period, while ROI 6 detects only one that also coincides with the runoff period. According to what was observed in Sentinel-2 snow-covered maps, on both dates identified as positive anomalies, the ROIs presented different percentages of remaining snow mixed with bare soil. As mentioned by Buchelt et al. [25], other conditions can be applied to identify events of progressive reduction of the snow cover, which may indicate more than one positive anomaly. Also, it is possible to identify events of melting followed by a freeze-up of the snowpack, as they cause a similar temporal evolution of the backscatter signal and, thus, may lead to false detection [25]. It is possible to disregard the second anomaly if two consecutive positive or negative anomalies are detected.
The results of Table IV show that with both methodologies, it is possible to detect two phase changes in seasonal snow-covered surfaces (bare soil to dry snow and dry snow to wet snow period) in the study area. However, our method contemplates each ROI's temporal and spatial variability using the derivative's standard deviation to calculate the phase changes. Our method detects events in ROIs where the fixed threshold methodology fails. The methodology presented in this study is robust and easy to implement compared with the methodology proposed by Tsai et al. [45] where interferometric, polarimetric, and backscatter features, along with elevation and land cover information, are needed. Buchelt et al. [25] proposed a method to detect the beginning of runoff using thresholds based on the seasonal minima of the snow SAR signal. In comparison to their study, our method considers less evident variations in the time series, being able to detect local minimums and, therefore, intermediate phase changes.

VI. CONCLUSION
Monitoring spatiotemporal changes in seasonal snow helps to predict and mitigate floods, avalanches, frost losses, impurities presence, and water supply, among other applications.
The combined use of field in situ measurements, meteorological data, remote sensing derived products, and albedo simulations has shown to have exceptional advantages in representing snow parameters in large and inaccessible mountainous areas.
In this work, we have seen that the seasonal backscattering characteristics of different ROI depend not only on mountainous surface features, such as slope and aspect, but also on surface roughness, soil moisture, LIA, and surface cover type, among others.
Using density curves and boxplots analysis, it was possible to determine that the statistics from bare soil, dry snow, and wet snow periods of some ROIs in Bariloche are significantly different. Such differences have not been reported in the existing bibliography in other places. However, the lack of field measurements is a limitation of this study. It is not possible to state the reasons behind the reduction of backscattering from the dry snow period with the available data.
We also showed that the mentioned surface cover differences represent limitations when using the bitemporal change detection method for wet snow mapping and also when using the automatic snowmelt phase detection algorithm based on a fixed threshold.
Instead of fixing an average threshold that would detect only the snowmelt phase, we developed a robust methodology that compensates for the heterogeneity of the study area. Our proposal can detect phase changes between the bare soil and dry snow period, as well as the dry snow to wet snow period.
This approach can identify the snow and land cover periods using positive and negative anomalies over SAR backscattering time series derivatives. It is a fast and adaptable approach for detecting phase changes since it does not consider a fixed threshold but the standard deviation of each derivative time series as a metric.
When compared with a fixed threshold methodology, the derivative method can detect two phase changes in 67% of the cases, while the fixed threshold methodology only works in 41% of them.
It was also possible to identify precipitation (rain or snow), melting, and refreezing events within the periods mentioned due to the variations of the backscattering series. These events were corroborated with meteorological satellite data, indicating that the ROIs' median backscattering values of the time series are sensitive to their overall environmental conditions and can adequately represent the land cover state at the sensor's acquisition moment.
This study advances the characterization of the Argentinean snowpack dynamics since the Patagonian Andes is one of the relatively most minor studied mountain environments globally.
Future work will compare these results against other innovative methodological approaches related to SAR and optic time series analysis and with field campaign data.