Deriving Regional Snow Line Dynamics during the Ablation Seasons 1984 – 2018 in European Mountains

Snowmelt in the mid-latitude European mountains is undergoing significant spatiotemporal changes. Regional snow line elevation (RSLE) is an appropriate indicator for assessing snow cover variations in mountain areas. To derive regional snow line dynamics during the ablation seasons 1984–2018, the present study unprecedentedly introduced a readily applicable framework. The framework constitutes four steps: atmospheric and topographic correction, snow classification, RSLE retrieval, and regional snow line retreat curve (RSLRC) derivation. The developed framework has been successfully applied to 8641 satellite images acquired by Landsat, ASTER, and Sentinel-2. The results of the intra-annual regional snow line variations show that: (1) regional snow lines in the Alpine catchments preserve the longest; (2) RSLEs are lower in the northern Pyrenees than in the southern part; (3) regional snow lines persist the shortest in the Carpathian catchments; and (4) during the end of the ablation season 2018, intermediate snowfall events in the catchments Adda, Tagliamento, and Uzh are observed. In terms of the long-term inter-annual variations, significantly accelerating snow line recession is detected in the northern Pyrenean catchment Ariege. In the Alpine catchment Alpenrhein and Drac, RSLRCs are shifting towards lower accumulated air-temperature (AT) significantly, with the magnitude of −3.77 ◦C·a−1 (Alpenrhein) and −3.99 ◦C·a−1 (Drac).


Introduction
Snowmelt during ablation seasons is essential with regards to runoff generation, winter sports/tourism, biodiversity, and natural disasters.Due to ongoing climate change, snowmelt is undergoing significant spatiotemporal variations [1][2][3][4].In the Northern Hemisphere, IPCC [5] projected with a high confidence that snow cover area (SCA) thereof was continuously decreasing during the ablation seasons, particularly in mid-latitude mountains [5,6].In Europe, mountain areas are one of the most climate-sensitive and vulnerable regions [7].They are critical habitats and natural water reservoirs and provide various ecosystem services and economic well-being.It should also be noticed that the majority of the annual runoff in European mountain areas is snowmelt-dominated [8].It is challenging to predict snow cover changes in mountain areas at a regional scale, given their different responses to climate change in temperature and precipitation [5,9].In these regards, in order to better understand regional responses and support adaptation strategies in the context of climate change, reliable information based on long-term time-series analysis of snowmelt processes during the ablation seasons is important for researchers, decision-makers, and stakeholders.
Earth Observation (EO) has a long history of application in snow monitoring.To date, different EO techniques (e.g., optical, radar, passive microwave, altimetry) have provided a timely, promising, and efficient approach to retrieve snow dynamics, including snowmelt processes [10].Previous studies characterized snowmelt processes with several spatial parameters: SCA [11,12], snow line elevation (SLE) [13][14][15], or with temporal parameters: snowmelt onset (SMO) [3,16], snow cover duration (SCD) [17,18], and snow persistence (SP) [19,20].Besides, to further capture the spatiotemporal dynamics of SCA throughout an ablation season, snow depletion curves (SDCs) are often utilized [21,22].Moreover, SDCs can also be employed in calibrating hydrological models (e.g., snow runoff model, SRM) with regards to discharge forecasting [23,24].However, it remains challenging to perform climate-change-related snow studies in complex terrains based on EO data, since long-term records (>30 years) in detailed spatial resolution (<100 m) are often required [25,26].To this end, EO is facing a typical challenge regarding the trade-off between spatial and temporal resolution, since none of the prevalently employed optical EO archives, i.e., Landsat, MODIS (moderate resolution imaging spectroradiometer), Sentinel-2, or AVHRR (advanced very high resolution radiometer), could heretofore satisfy such requirements for snow dynamics monitoring in European mountains.
To date, Landsat has the longest uninterrupted operation history, whose access became free-of-charge since 2008 [27].However, cloud obstruction and the near-two-week revisit time reduce the density of the Landsat time-series, and make it challenging to monitor the long-term snow dynamics in mountain areas solely using Landsat-based time-series.Firstly, cloud obstruction reduces the number of valid pixels in Landsat scenes.To reduce the influence of this issue, regional snow line elevation (RSLE) is introduced as an alternative to frequently-used SCA, since RSLE does not require fully cloud-free imagery.RSLE is calculated from a dynamic statistical quantile of elevation of snow/land pixels distribution [15].Secondly, the near-two-week revisit time increases the uncertainty in snowmelt observation, given that snowmelt is highly temporal dynamic and RSLE variation is not always monotonous during an ablation season due to potential intermediate snowfall events.To deal with this problem, accumulated air-temperature (AT), instead of date, is linked to snow dynamics to generate modified snow depletion curves (MSDCs) [28].It is because of the high likelihood that the relationship between AT and snow dynamic metrics (e.g., SCA) during an ablation season is monotonous.Ultimately, combining RSLE and MSDCs may provide a potential method to characterize long-term snowmelt processes during ablation seasons in mountain areas.
The overall objective of the present study is to derive 35-year regional snow line dynamics during the ablation seasons utilizing (semi-) high-resolution and free-of-charge optical dataset.For this purpose, a readily applicable framework has been developed.The framework takes advantage of the 30 m spatial resolution and more than three-decade observation history of the Landsat archive, which is also the main data source for the present study.Since the Landsat archive alone is limited by cloud obstruction, RSLE is used for deriving snow line information under partial cloud-covered scenes.To reduce the effects of near-two-week revisit time, the time domain is replaced by AT to generate regional snow line retreat curves (RSLRCs).Also, we combined Landsat observations with the images acquired by optical sensors with similar configuration (i.e., advanced spaceborne thermal emission and reflection radiometer, ASTER), to further densify the time-series.Finally, the obtained RSLRCs are then used to illustrate the timing of the ablation season and the speed of regional snow line retreat in relation to the air temperature.

Satellite Data and Pre-Processing
In the present study, three different (semi-) high resolution and free-of-charge optical data (see Table 1) are employed, i.e., (1) Landsat thematic mapper (TM), enhanced thematic mapper plus (ETM+), operational land imager/thermal infrared sensor (OLI/TIRS) between 1984 and 2018; (2) ASTER between 2000 and 2007 (since the ASTER/SWIR sensor was defected in late April, 2008, details see: https://asterweb.jpl.nasa.gov/swir-alert.asp);and (3) Sentinel-2A/2B multispectral instrument (MSI) between 2015 and 2018.The processing levels of employed Landsat, ASTER, and Sentinel-2 products are radiometrically calibrated and orthorectified (i.e., Landsat-L1TP, AST_L1T, and Sentinel-2-L1C).To obtain physically comparable surface reflectance datasets, atmospheric corrections were performed with ATCOR-3 [29] for each Landsat, ASTER, and Sentinel-2 scene.Topographic corrections were also integrated into ATCOR-3 by employing slope and elevation information from the ASTER Global Digital Elevation Model Version 2 (GDEM V2) [30].In this study, Landsat and ASTER data are used for analyzing long-term snow dynamics, since they are of a longer operational time span than Sentinel-2.Whereas taking advantage of shorter revisit time of combining Sentinel 2A and 2B data, the Sentinel-2 acquisitions have been used as an independent dataset for cross-validating the Landsat/ASTER-based RSLRCs.

Auxiliary Data
For the purpose of accuracy assessment, the snow cover classifications are compared with the snow depth records from the National Oceanic and Atmospheric Administration -Global Historical Climatology Network (NOAA-GHCN, https://www.ncdc.noaa.gov/ghcnd-data-access,accessed on 06 March 2019) and the European Climate Assessment & Dataset (ECA&D, https://www.ecad.eu/dailydata/predefinedseries.php,accessed on 06 March 2019).Therein, the snow depth records between 1984 and 2018 originate from 65 meteorological stations within the selected study areas.Together with the validated snow cover classification, elevation information is indispensable to determine the RSLEs.The ASTER GDEM V2 [30] is hence employed to assess the statistical quantile of elevation of snow/land pixels distribution for RSLEs derivation.To generate RSLRCs, 2 m air-temperature reanalysis data between 1984 and 2018 are obtained from ERA-Interim [31].Because glaciers influence snow line dynamics at a regional scale, it is required to calibrate RSLRCs in glaciated regions with the elevation information of the glacier(s).This elevation information is derived from the Randolph Glacier Inventory (RGI) and the ASTER GDEM.

Study Areas
Snowmelt processes are investigated in 11 montane catchments distributed across the mid-latitude European mountains (i.e., the Alps, the Carpathian Mountains, and the Pyrenees) during the ablation seasons (April to June) between 1984 and 2018.The numbers of analyzed Landsat, ASTER, and Sentinel-2 scenes for each catchment are listed in Figure 1.It should be noticed that the horizontal overlapping area of Landsat footprints is relatively large in mid-latitude regions, and the acquisition time difference between the two horizontally adjacent Landsat scenes thereof is nearly one week.If the major part of a catchment is within such an overlapping area, the data available in these catchments could be doubled.On the other hand, for large catchments that span over several horizontal overlapping footprints, it would be recommended to divide such catchments into sub-catchments.Otherwise, the mosaicked images from the Landsat scenes acquired in a one-week time difference may create artifacts in the retrieved snow dynamics.

Methods
The methods can be categorized into four main steps: (1) pre-processing, (2) snow classification, (3) RSLE retrieval, and (4) RSLRC, as described in the corresponding subsections (Figure 2).In the first step, 8641 satellite images were pre-processed, including 3903 Landsat, 726 ASTER, and 4012 Sentinel-2 images.Thereafter, the reflectance of the images was normalized.The pre-processed images were then classified into four classes: snow, land, cloud, and water/shadow (Section 3.1).The classification results were subsequently validated against the meteorological station observations.Based on the validated classification results and the DEM, RSLEs were retrieved catchment-wisely according to the elevation distribution of the snow and land pixels (Section 3.2).Afterward, the accuracy assessment was performed based on two quality indices showing the percentage of valid pixels and erroneous pixels.In the last subsection, RSLRCs were calculated by linking sigmoid-transformed Landsat-and ASTER-derived RSLEs to the contemporary AT using linear regression.The results are cross-validated against the corresponding Sentinel-2-based RSLRCs.

Snow Cover Classification and Validation
To classify snow cover, there are two main processing parts, i.e., using a decision tree based on multiple thresholds to classify snow, and using masks (i.e., cloud mask, water mask, shadow mask and/or thermal mask) to exclude misclassifications:

•
Snow classification: In this study, the algorithms developed by Klein et al. [32] and Poon and Valeo [33] is employed to classify snow.The algorithm is based on a decision tree with multiple thresholds on the normalized difference snow index (NDSI), the green band, and the near infra-red (NIR) band.To detect snow in forested areas, the NDSI-NDVI (normalized difference vegetation index) field is utilized to calibrate the snow cover classification results therein.However, such high NDSI values could also be observed in clear water bodies.Therefore, water bodies must be masked out to avoid misclassification.Because the water bodies commonly show positive normalized difference water index (NDWI) values, and the reflectance of water bodies in the green band is relatively low, the water mask is generated based on thresholding these two values.

•
Shadow mask: Shadow-cast areas are normally treated as non-valid pixels.In this study, the shadow pixels are identified following the methods from the ESA satellite snow product intercomparison and evaluation exercise (SnowPEX) Team [40].Thereafter, the shadow-cast pixels are masked out in the snow cover results.

•
Thermal mask: Both Landsat and ASTER have thermal bands.To filter out bright and warm surfaces such as warm rocks in the classification results, a thermal threshold (<288 K) introduced by Metsämäki et al. [41] is applied to Landsat-and ASTER-based snow classifications.Sentinel-2 does not have any thermal band, which could potentially commit more commission errors over bright and warm targets.
The workflow chart is displayed in Figure 2, where snow classifier and masks are applied in sequence to obtain snow classification results.Thereafter, the binary snow classification results are validated using the contemporary NOAA-GHCN and ECA&D snow depth observations.The validation is based on the snow depth threshold proposed by Parajka et al. [42], i.e., if the observed snow depth is no less than 1 cm, the corresponding pixel is regarded as snow-covered, and vice versa.

Regional Snow Line Elevation Retrieval and Accuracy Assessment
The definition of the snow line elevation may vary according to objectives of the different applications (e.g., geographical studies, meteorological studies, hydrological studies) as discussed by Krajčí et al. [15].In this study, we used the term RSLE proposed by Krajčí et al. [15] that is originally designed for remote sensing applications.RSLE is defined as the elevation where there are as few as possible snow pixels below it, and as few as possible land pixels above it.Methodologically, RSLE can be determined as the elevation where the minimum value of the sum of two cumulative histograms (i.e., cumulative histogram of snow pixels elevations and land pixels elevations) is reached (Figure 3). .Estimation of regional snow line elevation (RSLE) from the combined cumulative histograms of the snow-covered pixels (in blue) and land pixels (in orange).In the x-axis, N 1 indicates that there are N 1 land pixels above the RSLE; (N 2 -N 1 ) there are N 2 snow pixels below the RSLE, N 2 alone means the sum of land pixels above the RSLE and snow pixels below the RSLE, and TP stands for the total number of pixels.
To assess the accuracy of the retrieved RSLEs, two quality indices are introduced.The first quality index is representativeness index (RI) measuring the percentage of valid pixels (i.e., labeled as snow or land) within the spatial extent of a catchment.The calculation of RI is expressed as: where F and S are the number of snow-free and snow pixels within the catchment extent in classified images.T p represents the total number of all the pixels (i.e., snow, land, cloud, water, and shadow) within the corresponding catchment extent.The optimal RI value is 100%, which indicates that there are no invalid pixels (i.e., cloud, shadow, water, missing value) within the catchment for RSLE retrieval.The second quality index, error index (EI), measuring the percentage of erroneous pixels (i.e., the number of the snow pixels below the RSLE, noted as S b , and the number of snow-free pixels above the RSLE, noted as L a ) according to the corresponding RSLE.EI can be calculated as: By definition, the obtained RSLE should result in as few erroneous pixels as possible.EI is introduced as the ratio between the number of erroneous pixels and the total number of pixels (T p ), hence the best EI value is 0%.

Regional Snow Line Retreat Curve (RSLRC) Derivation and Validation
Once the RSLEs during each ablation season are retrieved, RSLRCs can be calculated, which describe the recession of the regional snow line over the ablation season.The RSLRC links the AT and RSLEs using the robust M-estimation [43].The typical shape of a RSLRC is sigmoid, and therefore a logistic link function is firstly applied (Equation ( 3)).The robust regression is then implemented using the 'rlm' library in R. The formulas of the RSLRC can be expressed by Equations ( 3)-( 5): where the RSLE max is the highest RSLE, which equals to the highest elevations from the DEM.In case of the glaciated areas, RSLE max is estimated by the 95th percentile of the elevations in delineated glacier outlines from the RGI.The slope (k) is always a negative value, whose absolute value represents the steepness of the RSLRC.The ratio between the intercept (b) and -k (Equation ( 4)) is a coefficient indicating the AT of the mid-ablation season (AT MA , the mid-point of the RSLRC).The generalized behavior of the RSLRC in relation to the coefficients is illustrated in Figure 4. AT i represents the AT at the ith day within the ablation season, which is the integral of the daily temperature above the base temperature (T 0 ) added by the AT from the previous month (AT 0 ).This calibrates the shift of consecutive temperatures below 0 • C at the beginning of an ablation season.In this study, AT i is approximated by the daily average of four 2 m air-temperature measures at a 6-h interval, and the base temperature is equal to the melting point 0 • C for snow/ice (Equation ( 5)).To guarantee the quality of the input RSLEs, a threshold of RI > 20% has been set to filter out RSLE results that may not be representative due to the lack of valid observations.Moreover, an end-ablation snow line correction is involved in the processing chain.It detects all input RSLEs within the maximum ± 10% of the elevation range in the catchments.If the number of end-ablation RSLEs is more than 80% of the total observation amount, then the RSLRCs would not be generated as too few observations of the regional snow line retreat are captured.If the number of end-ablation RSLEs is more than 50% of the total observation amount, only half of the RLSEs are included in the RSLRCs.This is because the breakdown point of the robust M-estimation is approximately 0.5, meaning more than 50% of end-ablation RSLEs would provide erroneous weights to the observations.With regards to the accuracy assessment, the corrected coefficient of determination (corrected-R 2 ) for robust regression [44], mean absolute error (MAE) and root mean square error (RMSE) are calculated to evaluate the retrieved RSLRCs.To date, since the Sentinel-2 data only cover the ablation seasons 2016-2018, it does not bring a significant improvement to the 35-year time-series.On the other hand, given that the Sentinel-2 archive has the best data availability within its available years, Sentinel-2-derived RSLRCs constitute a valuable independent dataset to cross-validate the Landsat/ASTER-derived RSLRCs.

Results
The present study aims to characterize regional snow line dynamics during the ablation seasons.Within this context, two aspects of the regional snow line dynamics are presented in the following subsections, i.e., intra-annual and inter-annual variations of regional snow lines.Also, a comprehensive accuracy assessment is provided in the last subsection.

Intra-Annual Variations of Regional Snow Lines during the Ablation Season 2018
To demonstrate the intra-annual variations of the regional snow lines, RSLEs during the ablation season 2018 are displayed in Figure 5.The RSLEs are derived from the Landsat OLI/TIRS and ETM+ observations, representing the general situation when dual Landsat/ASTER sensors are in orbit.Given that the acquisition dates of Landsat images vary spatially, an ablation season is divided into nine categories (i.e., early/middle/end of April/May/June) in a 10-day time interval for better demonstration.Moreover, to better visualize the spatiotemporal distribution of the RSLEs, an overview of the Alpenrhein catchment is additionally illustrated in Figure 5. Spatially, in the northern Pyrenees (i.e., Ariege) RSLEs are lower than in the southern part (i.e., Serge).Besides, the regional snow line in the northern part of Pyrenees is preserved longer compared to the southern part.In the two Carpathian catchments, particularly in Uzh, the RSLEs last much shorter than in other investigated catchments during the ablation season 2018.In contrast, the regional snow lines are still preserved at the end of the ablation season 2018 in most of the Alpine catchments.Regarding the temporal pattern, the occurrence of the reddish color indicates the snow line at the end of the ablation season (i.e., the end of June), which should be located around the highest elevation.Otherwise, once the regional snow lines in reddish color are located in lower elevation zones compared with previous ones, it indicates the occurrence of an intermediate snowfall event.Therefore, these results could help identify anomalies, such as intermediate snowfall events observed in Figure 5, where reddish/yellowish regional snow lines are of lower elevation than bluish/greenish ones.For example, the yellow-colored (end of May) regional snow line appears in Adda and Tagliamento, and the coral-colored (end of June) regional snow line is observed within Uzh.In addition, within the Carpathian region, regional snow lines are rarely observed at the end of the ablation season.

Inter-Annual Variations of Regional Snow Lines during the Ablation Seasons 1984-2018
To demonstrate the inter-annual regional snow line variations, the time-series of the AT MA and steepness from RSLRCs are shown in Figure 6 for each study area.In terms of the amount of the usable RSLRCs after the quality controls (Section 3.3), the Alpine catchments, as well as the northern Pyrenean catchments, are of the highest RSLRC quantity (29 years in average).Fewer RSLRCs are available for the catchments in the Carpathian Mountains and the southern Pyrenees (16 years on average).RSLRCs' scarcity appears mostly in the 1990s, particularly with the severest scarcity in the Carpathian Mountains.With such few RSLRC results, it is difficult to carry out statistically trustable analysis.Yet in most of the Alpine catchments and northern Pyrenean catchments, performing statistical analysis is still feasible.
The steepness of the RSLRC represents the velocity of regional snow line retreat because the AT is related to the day-of-year in the ablation season.Among the investigated catchments, only in Ariege (in the northern Pyrenees) a statistically significant (p-value = 0.036) negative trend of the RSLRCs' steepness has been observed, based on the 28-year observations.It indicates that the snow line is retreating significantly faster in this northern Pyrenean catchment.The 35-year AT time-series (see Figure 7) is indicating a tendency of increasing AT during the ablation seasons.It is thus also necessary to investigate the timing of the snow-clearance process.AT MA illustrates how much AT is required for a regional snow line to reach the middle of the RSLE.Significant negative trends of AT MA are detected in some Alpine catchments distant to the Mediterranean sea, i.e., Alpenrhein (p-value = 0.038) and Drac (p-value = 0.052).It indicates that lower than 3.77 and 3.99 Celsius per year are needed for the regional snow lines to reach the middle of the RSLRC in the Alpenrhein and Drac, respectively.As in these two catchments, no significant trend of RSLRCs' steepness is observed, it also suggests higher RSLEs in early April, and the disappearance of snow line retreat occurs earlier.

Accuracy Assessment
The uncertainties of the derived RSLRCs and the corresponding parameters are mainly caused and propagated by four sources: (1) misclassifications in the SCA maps, (2) errors of the employed DEM, (3) uncertainties of the RSLE results, and (4) imperfect model fit.Given that the accuracy of ASTER GDEM has been reported by several previous studies (e.g., [30,45,46]) and would stay systematic in the analyses, we only assess the three remaining error sources.

Accuracy Assessment of the Snow Classification Results
The overall accuracy (OA) of the SCA maps has been reported as and the Kappa coefficient is 0.72, according to 7720 ECA&D and NOAA-GHCN snow depth observations.In addition, based on the confusion matrix (Table 2), precision and recall are calculated.Precision equals to 83.25%, indicating that 83.25% of the classified snow pixels are truly snow-covered.The recall is 65.41%, meaning that 65.41% of the snow-covered pixels have been detected.Regarding the accuracies of the SCA maps derived from the acquisitions from different sensors, the OAs are reported as 94.38% (ASTER), 97.08% (Landsat), and 95.22% (Sentinel-2).The accuracy assessment is implemented by investigating two quality indices (i.e., RI and EI in Figure 8) introduced in Section 3.2.The median percentages of the valid pixels within each catchment are above 60%, and the upper and lower quantiles thereof are near 90% and 40%, respectively.From the boxplots of the EIs, in general, the obtained RSLEs are of low (<5%) erroneous pixels percentage, i.e., the snow pixels below the corresponding RSLE and the land pixels above the corresponding RSLE.

Accuracy Assessment of the Regional Snow Line Retreat Curves (RSLRCs)
To assess the performance of the robust M-estimation of the derived RSLRCs, the corrected-R 2 for robust estimation, MAE, and RMSE are calculated for each fitted RSLRC (Figure 9).The three metrics summarize the fit of the RSLRC to the RSLE inputs.Within the Alpine catchments (excluding Tagliamento), the median corrected-R 2 is around 0.90.Therein, more than 90% of the variance in the RSLEs can be explained by the variation of AT.According to the median, the RSLRCs explain approximately 10% less variation in the RSLEs in the Pyrenean catchments than the aforementioned Alpine catchments.In the Carpathian Mountains and Tagliamento, most of the RSLRCs are only able to predict 65% of the RSLE variability according to their comparably small median corrected-R 2 (near 65%).from Maira, the upper-quantile, median, and low-quantile of MAEs are generally below 30 m, 15 m, and 10 m, respectively.Most of the RMSEs are lower than 300 m.Maira is the catchment showing a comparably higher MAE and RMSEs than the other investigated catchments.To cross-validate the RSLRC steepness and AT MA obtained from the Landsat/ASTER-based RSLRCs in the past 35 years, Sentinel-2 data are employed to calculate the contemporary RSLRC steepness and AT MA following the same workflow (Figure 2) between 2016 and 2018.The performance of each RSLRC using cross-validation is displayed in Figure 10.The cross-validation shows a high agreement between the RSLRCs calculated from Landsat/ASTER imagery and Sentinel-2 imagery.Regarding the bias of the estimated RSLRC steepness, Landsat/ASTER has an approximately 2.5% overestimation than the Sentinel-2.Yet near 3.5% underestimation is made by Landsat/ASTER-derived RSLRCs with regard to AT MA .The coefficient of correlation (r) regarding the RSLRC steepness and AT MA between Landsat/ASTER-based RSLRCs and Sentinel-2-based RSLRCs are 0.75 and 0.88, respectively, indicating a good agreement.

Discussion
In this section, we first discuss the challenges in accurately deriving RSLEs and RSLRCs with regards to data data reliability, and method efficiency.Afterward, the limitations and difficulties concerning the validation data and validation scheme are summarized.Lastly, we discuss the observed regional snow line variations and the potential applications of our results.
5.1.Challenges in Accurately Deriving Regional Snow Line Elevations (RSLEs) and Regional Snow Line Retreat Curves (RSLRCs) The uncertainties of the retrieved RSLEs mainly come from the accuracies of snow classifications and the involved DEM.The accuracies of snow classifications are affected by 1) cloud contamination, 2) signal saturation, 3) thermal-band absence, and 4) snow-in-forest detection.Among these factors, cloud contamination is the most influential.Commission errors are often recognized over bright targets including snow/ice [39,47,48] when applying the cloud masks (e.g., MFmask, ACCA).The misclassified pixels near the snow lines could potentially bias the statistics/distribution of the combined cumulative histograms of the snow and land pixels (Figure 3).Signal saturation (e.g., in Landsat 4/5 TM, Landsat 7 ETM+, ASTER [49,50]) could lead to poor accuracy in snow classification because the visible bands are often used to separate snow/cloud from other bright surfaces [51].It is in line with the accuracy assessment results in Section 4.3.1.Moreover, it should be noticed that the thermal band is absent in Sentinel-2 data.It potentially leads to more commission errors when detecting snow in regions surrounded by warm bright land covers.Besides, detecting snow-in-forest using optical satellite data is challenging, since forest canopy covers up the snow beneath and casts shadow on its surroundings [32,52].Obtaining accurate DEM in mountain areas are problematic for either interferometric method or photogrammetric method based on spaceborne satellite data.Often, the correlation algorithm for photogrammetric DEM generation fails over snow-covered areas [53,54], while an interferometric DEM has the problems of the voids in the original data [55].In the present study, the errors induced by the employed DEM stay systematic, since they are largely counteracted.Given that airborne LiDAR can produce a DEM with high accuracy, further studies with regards to accuracy improvement should consider airborne LiDAR DEMs as valuable alternatives to interferometric or photogrammetric DEMs.
The accuracies of the derived RSLRCs depend on the availability and the representativeness of the input RLSE data, as well as the efficiency of the regression model.Data availability of spaceborne optical images varies spatiotemporally (especially in high mountain regions), due to cloud cover, acquisition plans (e.g., "commercialization era" [56]), sensor anomalies (e.g., Landsat 7 ETM+ SLC-off issue [57]) and footprint patterns [13].Therein, cloud obstruction is the major issue, as clouds usually persist for a long time during the winter in European mountains.The resultant missing observations lead to the failure in RSLRC derivation, as shown in Figure 6.Moreover, footprint pattern is strongly influencing the data availability at a catchment scale.Taking Landsat as an example, the acquisition date difference of two horizontally adjacent footprints is often one week, thus data availability therein can be potentially doubled.In terms of the representativeness of the input data, apart from the accuracy of input RSLEs (Section 5.1), the accuracy of RSLRCs is also associated to the intermediate snowfall events, which create anomalously low RSLEs that are often treated as the 'outliers' in the RSLRCs.Since there is a high probability of intermediate snowfall events during accumulation seasons and the beginning of the ablation seasons, the present method is more suitable for middle-to-late ablation seasons.Furthermore, the heterogeneous geographical settings (e.g., aspects, solar irradiation, slopes, climate characteristics) are influential, as RSLEs are regional products.However, due to cloud cover, splitting a catchment into several sub-zones can result in insufficient cloud-free pixels for RSLE derivation, especially around the cloud-prone windward slopes.Regarding the representativeness of the climate reanalysis data, the original spatial resolution of the applied ERA-Interim is 80 km.It can potentially induce problems for catchments with a large proportion of plain areas (e.g., Maira catchment) because the obtained AT is, in fact, more representative for the snow-free areas.Such coarse resolution also results in the omission of temperature variation in complex terrains.Further researches employing higher resolution climate reanalysis e.g., ERA-5 reanalysis dataset (32 km), COSMO-REA6 (6 km), COSMO-REA2 (2 km), are desirable.Noteworthily, these higher resolution climate reanalysis data are often only available for a limited time or spatial extent.
It must be noticed that the main objective of generating RSLRCs in the present study is not to predict intermediate RSLEs between two ATs (or two dates), but rather to characterize the regional snow line behavior during an ablation season.Therefore, the most crucial aspects of the regression models are breakdown point and goodness of model fit.Robust regression models (i.e., robust M-estimation) often have a higher breakdown point, meaning they could handle a larger amount of contaminated observations (e.g., RSLEs after an intermediate snowfall event) than conventional ordinary least squares (OLS) regression.Therefore, they can reduce the influences of high leverage outliers (e.g., intermediate snowfall events), and avoid biased-steepness from the logistic curve.In addition, to further alleviate the impacts of intermediate snowfall and temperature anomalies, the time-domain in RSLRCs is replaced by AT.To eliminate contaminated input data, RSLEs are filtered based on RI.To address the second concern (i.e., the goodness of model fit), a threshold of corrected-R 2 > 0.4 is chosen to exclude the RSLRCs that predict less than 40% of the variability of the regional snow line retreat.Section 4.3.3further confirms that the involved RSLRCs are mostly of a corrected-R 2 > 0.85 in Alpine catchments, indicating a good fit of the applied regression models.Given that it is not appropriate to evaluate the model fit solely with the coefficient of determination, we have processed the Sentinel-2 time-series in the last three years to carry out a cross-validation to check the robustness of the RSLRCs.The high agreement between Sentinel-2-based and Landsat-based RSLRCs indicates that the method is valid for catchments where RSLRCs can be obtained under the abovementioned constraints.

Challenges of Validation
Validating long-term RSLEs and RSLRCs from (semi-) high-resolution optical satellite data is extremely challenging due to the lack of data and validation techniques.Conventionally, the validation of snow-related products can be realized using field measurements or satellite data of higher resolution.For instance, Krajčí et al. [15] combined the automatic weather station (AWS) with snowfield measurements to validate the MODIS-derived RSLE results.However, the elevation interval is 100 m and the study area is regional.For a larger spatial extent and finer resolution, a considerable amount of field measurements are required, which is unpractical due to intensive labor and material requirements and low accessibility in mountain areas.In addition, meteorological station data may not be dense enough to validate RSLEs results retrieved from (semi-) high-resolution data.Also, neither historical field measurements nor meteorological station data can reach back to 1980s.Regarding the satellite data with higher spatial resolution, the drawbacks are: (1) they are usually not free-of-charge; (2) they only have visible/near-infrared (Vis/NIR) or even only visible bands, leading to great difficulties to separate snow and other bright targets; (3) The swath is much smaller than a Landsat footprint; (4) The satellite often follows an acquisition request and the revisit-time is much longer than Landsat; (5) the operational time-span rarely reaches back to 1980s.More importantly, such data are often acquired by optical sensors, which would have relatively similar cloud obstruction at a catchment extent at the same time.Therefore, the uncertainties below the cloud can be hardly revealed by such dataset.
Unmanned aerial vehicle (UAV) imagery is a potential dataset to validate the RSLE and RSLRC results independent of cloud obscurations.UAVs are flexible to deploy in the aspects of mounted sensors (e.g., multispectral, hyperspectral, thermal, microwave, LiDAR), acquisition date, and operation periods.The spatial resolution of UAV images is often higher than those acquired by satellites.Therefore, UAV imagery holds a great potential to provide highly detailed information about the snow cover, especially for cloud-covered areas.Alternatively, WebCAM data have the advantages of (1) cloud cover independence, (2) low-cost, (3) very high observation interval, and (4) relatively long history in observing some mountain regions.Apart from a relatively high demand in pre-processing, it is so far a unique potential dataset for validation.To homogenously produce such a WebCAM dataset, it requires a great among countries and institutes.Besides, spaceborne synthetic aperture radar (SAR) data can penetrate the cloud cover, thus revealing the snow cover condition under the clouds.Based on the current method [58][59][60], total snow cover areas could be derived from the SAR observation with an accuracy of up to 98.1%.However, due to the data access policy, it would be very costly to apply for these data for long-term research at a large spatial scale.
It should be noticed that the present work provides an exploratory framework to produce a long-term and (semi-) high-resolution-based snow line product with a fundamental validation process.We have noticed that the presented validation scheme is imperfect, but it is the maximum of what researchers can do using the free-accessible data.In this context, further studies in relation to the validation are urgently desirable.

Observed Regional Snow Lines Dynamics
In term of the intra-annual variations, the RLSEs between April and June illustrate the snow clearance in the studied catchments during the ablation season 2018.However, the number of retrieved RSLEs strongly depends on the elevation range and the number of valid pixels.In some catchments, the number of RSLE observations is limited, particularly in the Carpathian Mountains.Taking the Carpathian catchment Uzh as an example, its maximum elevation is around 1500 m while the elevations of the rest catchments range from 2200 m to 4100 m.The regional snow lines are seldom persevered after May in Uzh, due to such relatively low maximum elevation.Although it would be logical to switch the observation period towards mid-winter, the chance of intermediate snowfall events is high during that time, leading to very high uncertainty.Also, the high cloud persistence in mid-winter [13] would further decrease the data usability.Thus, for such catchment of lower elevation, the method is not well-suited.Meanwhile, as presented in Section 4.1, we have identified intermediate snowfall events in three catchments: Adda (end-of-May), Tagliamento (end-of-May), and Uzh (end-of June).To confirm these three observed snowfall events, we have checked the SCA conditions from global snow pack (GSP), a MODIS-based daily gap-filled SCA product [56], at a catchment scale.From GSP results, three snowfall events have been confirmed.
Inter-annual variations results indicate that there is a significant negative trend of the regional snow line retreat in the northern Pyrenean catchment Ariege.In the context of the increasing AT, it would result in a shortened ablation season, which is in line with the results from Buisan et al. [57].On the other hand, temporal shifts of RSLRCs have been observed.Based on our linear trend analysis on AT MA , there are significant negative trends in the Alpine catchment Alpenrhein and Drac.It indicates a tendency of a shortened ablation season, and higher RSLEs in April.The results are in line with previous studies [16,61,62].Meanwhile, we discovered an obvious RSLRC scarcity in the Carpathian catchments compared to the other investigated catchments.It indicates that our method is not well-suited for deriving RSLRCs for the Carpathian Mountains, as the robust M-estimation failed to predict the majority of the regional snow line recession represented by the low corrected-R 2 values.There are two main reasons, i.e., insufficient valid RSLE inputs and intermediate snowfall events [13].Besides, snow can be completely melted in the early months due to the low maximum elevation.The factors above make the RSLEs noisy and difficult to fit by a simple curve.
The presented results can be further used in applications with regards to snow cover phenology, climate change, winter tourism/sport, as well as flora and fauna.Snow cover phenology is a vital perspective of snow cover dynamics, which includes snow accumulation onset (SAO), snowmelt onset (SMO), snowmelt end (SME), snow cover duration, to name a few.However, given the revisit time and intermediate snowfall event occurrence, (semi-) high-resolution optical EO sensors like Landsat is rather for monitoring snow dynamics in middle-to-late ablation season.Since RSLEs indicate the spatial snow cover in the altitudinal direction, it could better illustrate the dynamics of snow clearance rather than the whole snowmelt process.This is because of the lower snowfall frequency and cloud coverage in the late ablation season.Furthermore, quantitate analyses linking RSLE dynamics with other ECVs is beneficial to enhance the comprehension of regional response climate change.Thereby, better adaptation strategies be made accordingly.Snowmelt run-off is vital to water management and disaster warning in snowmelt-dominant catchments.Run-off can be predicted by calibrating hydrological models using spatial snow coverage derived from the RSLRCs.It would be intriguing to compare the performance using (semi-) high-resolution-based snow information and medium-resolution-based information [63] in mountain areas.Winter tourism/sport management requires detailed spatial snow cover data for hiking routine planning, artificial snow planning, as well as the operation period length.Flora and fauna in high mid-latitude mountain regions are influenced by the snowline dynamics.The spatiotemporal snow dynamics affect the plant phenology [64][65][66], and animal activity [67][68][69].

Summary and Conclusions
This paper presents a readily applicable framework for retrieving regional snow line elevations (RSLEs) and their dynamics based on free-of-charge optical remote sensing and climate reanalysis datasets.Snow cover areas (SCAs) are firstly classified using a multi-threshold decision tree and corrected with cloud/shadow/water/thermal masks.The overall accuracy of the SCA results are around 96.71% according to 7720 meteorological snow-depth observations.The RSLEs are subsequently determined by the spatial distribution of the classified snow-covered and snow-free areas.The majority of the RSLEs are based on satellite images with more than 60% valid pixels (cloud-free and non-shadow) and are representative (>95% of the snow/land pixels) of the spatial distribution of snow in the altitudinal direction.Then the RSLEs are combined with accumulated air-temperatures (AT) to derive regional snow line retreat curves (RSLRCs) in a sigmoid shape.Based on the robust M-estimation, the RSLRC steepness and AT MA (AT of mid-ablation season) are derived to characterize the regional snow line dynamics during the ablation seasons 1984-2018.The applied robust M-estimation models predict approximately 90% (Alps), 80% (Pyrenees), and 65% (Carpathian Mountains) of the regional snow line retreat variations during the ablation seasons 1984-2018, according to the corrected coefficients of determination.Finally, the Landsat-derived results are cross-validated against the Sentinel-2-derived RSLRCs (2016-2018).The cross-validation shows that the Landsat-derived RSLRC steepness and AT MA have only 2.5% overestimation and 3.5% underestimation, respectively, indicating good agreement with those derived from Sentinel-2 (r = 0.75 for the RSLRC steepness, r = 0.88 for AT MA ).
This framework has been applied to 11 mid-latitude catchments in Europe.The results with regards to the intra-annual variation illustrate the recession of the regional snow line during the ablation season.The highest RSLRC quantity (29 years in average) is obtained in the investigated Alpine and northern Pyrenean catchments.The missing RSLRCs occur mostly in the 1990s, particularly in the Carpathian Mountains.Regarding the intra-annual variations, the results show that: (1) The Alpine catchments preserve regional snow lines longer than the other investigated catchments; The applicability of the presented framework is mainly constrained by data availability and intermediate snowfall events.The framework performs well in the Alpine catchments but problematically in the Carpathian catchments, where historical cloud-free Landsat data are scarce, and intermediate snowfall events occur frequently.The sparse observations and weak model performance also lead to difficulties in interpreting the statistical results.To address these problems, the next step would be to increase the number of EO data, and to involve climate reanalysis data with improved spatial resolution.The introduced RSLRCs are able to ingest snow classification results from other optical sensors, e.g., CBERS (China-Brazil Earth Resources Satellite program), SPOT

Figure 1 .
Figure 1.Overview of the study areas and number of available Landsat scenes from TM, ETM+, OLI/TIRS, ASTER and Sentinel-2 multispectral instrument (MSI) between April and June (1984-2018).

Figure 3
Figure3.Estimation of regional snow line elevation (RSLE) from the combined cumulative histograms of the snow-covered pixels (in blue) and land pixels (in orange).In the x-axis, N 1 indicates that there are N 1 land pixels above the RSLE; (N 2 -N 1 ) there are N 2 snow pixels below the RSLE, N 2 alone means the sum of land pixels above the RSLE and snow pixels below the RSLE, and TP stands for the total number of pixels.

Figure 4 .
Figure 4. Theoretical regional snow line retreat curve (RSLRC) based on the accumulated air-temperature (AT) and RSLE.AT EA is the AT at the end of the ablation season, whose RSLE is the RSLE max .AT min and AT max are the minimum/maximum AT observed within the ablation season.AT MA is the mid-point of the RSLRC, which is the AT of mid-ablation season.The solid line represents the standard situation, and the dash lines are the RSLRCs with different regression coefficients.The dash lines and dot-dash lines represent the behaviors of the RSLRCs in different slope and intercept values.

Figure 10 .
Figure 10.Comparison of estimated steepness of the RSLRCs and the AT MA obtained from Landsat/ASTER imagery and Sentinel-2 imagery.
(2) RSLEs are lower in the northern Pyrenees than in the southern part; (3) RSLEs last the shortest in the Carpathian catchments; and (4) intermediate snowfall events occurred in Adda and Tagliamento during the end of May 2018, and in Uzh during the end of June 2018.Meanwhile, significantly (p-value = 0.036) fastened snow line retreat is detected in the northern Pyrenean catchment Ariege.In the Alpine catchment Alpenrhein and Drac, there are statistically significant trends of the RSLRC shift towards lower AT, in the magnitude of −3.77 • C•a −1 (Alpenrhein) and −3.99 • C•a −1 (Drac).

Table 1 .
Configurations of utilized bands of the selected earth observation (EO) sensors in this study: Landsat thematic mapper (TM), enhanced thematic mapper plus (ETM+), operational land imager/thermal infrared sensor (OLI/TIRS), advanced spaceborne thermal emission and reflection radiometer (ASTER), and Sentinel-2 (S2).CW stands for central wavelength in µm, and SR stands for spatial resolution in meters.

Table 2 .
Confusion matrix relating satellite-derived snow classifications and ground snow-depth observations.