Snow Cover Variability in Brunswick Peninsula, Patagonia, Derived from a Combination of Spectral Fusion, Mixture Analysis and Temporal Interpolation of MODIS Data

1. Introduction

Snow cover is a critical component of the cryosphere interlinked with climate at local, regional, and global scales [1]. Through radiative and thermal properties, various Earth System components modulate the transfer of energy and mass at the Earth surface-atmosphere interface [2]. Snow-covered areas are sensitive to climate change [3,4] since the surface temperatures are close to the melting point and various feedback processes, thermal insulation properties and meltwater input make this sensitivity even stronger [5]. Changes in the seasonal behaviour of snow cover, such as melting patterns, soil moisture and water availability, can strongly affect both plant communities as well as the hydrological cycle of watersheds and their ecosystems [6,7,8].

Snow appears bright to the human eye because of its high and relatively flat reflectance in the visible region of the electromagnetic wavelength (400 – 700 nm), while it has stronger absorption in parts of the near and mid-infrared [9,10]. However, the pattern of the spectral response of snow depends on its density, grain size, liquid water content, and the presence of impurities such as dust, soot, and algae [9]. Several optical remote sensing data (such as AVHRR, MODIS, VIIRS, and Landsat) can represent the snow cover extent in different temporal, spectral, and spatial resolutions [9]. For instance, the Landsat suite of sensors, with a spatial resolution of 30 m in visible, near- and short-infrared bands, is ideal for snow applications in mountains. The spectral bandwidth of Landsat data allows for accurate retrievals of fractional Snow Covered area (“fSC”, as a subpixel fraction) [11,12,13]. However, the repeat-pass interval of 16 or 18 days is not adequate for short-term snow-mapping requirements [14], especially in cloudy regions. MODIS (Moderate Resolution Imaging Spectroradiometers) instead, with a spatial resolution of 250 m in 2 bands and 500 m in 5 visible and near-infrared bands, can monitor snow cover on a daily resolution since 2000, employing an automated snow cover algorithm developed at Goddard Space Flight Center (GSFC1) [14,15,16]. The standard MODIS snow map products include fSC, snow albedo and cloud-gap filling. Improvements to the algorithm are incorporated each time the entire dataset is reprocessed [14], which has currently been updated to collection version 6 (https://nsidc.org/data/MOD10A1). The basic algorithm (Eq. 1) for generating the operational snow cover product from MODIS data is the Normalized Difference Snow Index (NDSI), as described in Hall and Riggs (2011) [17], which is the most common method used for snow cover extent identification [18].

$$NDSI=\frac{\rho B_4-\rho B_6}{\rho B_4+\rho B_6}$$

(1)

where $\rho B_4$ is the green band reflectance (bandwidth 545–565 nm), and $\rho B_6$ is a mid-infrared reflectance (bandwidth 1628–1652 nm). A pixel with NDSI > 0.4 is assumed to have snow cover (binary result: snow or no snow), while a pixel with NDSI <= 0.4 is snow-free [19].

This global snow cover algorithm has had several refinements over the years and different thresholds have been incorporated into its general formula [14]. Nevertheless, it is widely known that the NDSI has problems in differentiating snow from clouds [15,19,20,21]. Consistent results have not yet been established in the actual collection (version 6) and the problems of the previous version remain [16].

The spatial resolution is something that must also be considered to assess the uncertainties of the remote sensing products in estimating snow cover. Hydrological studies and snow model validation in mountainous terrain need the greatest possible spatial detail, especially when snow interaction with vegetation is relevant [22,23]. This allows a better representation of snow season behaviour such as the snow line [24] and also more accurate discrimination of snow albedo [13,25]. Efforts to improve the spatial resolution of snow products are thus highly desirable. For snow map studies in complex terrain The Global Climate Observing System [26] recommends a spatial resolution less than 100 m, while in flat terrain 1 km is considered adequate. For the case of complex terrain snow distribution can be assessed using different techniques and sensors, for example a combination of Landsat 30 m resolution and 16 days repeat pass with MODIS at 250 m resolution and daily temporal resolution.

The development of accessible and comprehensive methods to reconstruct snow properties using satellite remote sensing data is an emerging topic, particularly regarding interaction with vegetation. For instance, Sirguey et al. (2009) [24] have developed a comprehensive method to produce regional maps for New Zealand of snow-cover extent at 250 m spatial resolution, using the MODIS Level-1B data product with corrected reflectivity. They showed that subpixel snow fractions could be retrieved with a Mean Absolute Error (MAE) of 6.8% at 250 m. Painter et al. (2009) [12] improved products from MODIS (MOD09GA) by applying a combination of atmospheric correction and a physical algorithm called MODSCAG (MODIS Snow-Covered Area and Grain size), which yields an average RMS error of 5% for snow cover at a spatial resolution of 500 meters. MODSCAG has the advantage of detecting different kinds of snow at subpixel level and thus albedo. In comparison, the MOD10 snow albedo product (version 5) has a MAE about three times higher, and no improvement has been achieved in the actual collection (version 61) [16]. The advantages of discriminating snow-cover extent, snow grain size, and albedo are very important, considering that they are critical information for studies on energy and mass balance [27,28]. However, many recent studies are still using the classical MOD10 products to both reconstruct snow cover variability [29,30] and validate snow models [31,32].

Southern Patagonia is a key location for climate change studies, particularly spatio-temporal snow cover variability, due to its unique and continuous extent in the Southern Hemisphere beyond 45ºS [33], reaching 56ºS at Cape Horn [34]. Southernmost South America is considered a real “hotspot” of paleoclimatic records, that have recorded, on a scale of hundreds to thousands of years, abrupt changes in circulation regime and atmospheric temperature associated with global climate variability, including ice cores, geomorphological features, peat bog systems, marine and lake core sediments, and tree rings [35,36,37,38,39,40].

The climate and topography interactions in southernmost South America are largely dominated by westerly winds located between the high-pressure centre of the Southeast Pacific Ocean (SPH) and the circumpolar low-pressure belt surrounding Antarctica [41]. The Southern Annular Mode (SAM) is the main mode of large-scale variability that controls the strength of the westerlies [42]. These westerlies are disturbed by both the southern Andes and the Fuegian Andes (or Darwin Cordillera), producing one of the most extreme precipitation gradients on Earth [43,44,45,46]. A recent study based on a 17-year long meteorological record of a west-east longitudinal transect has recorded an extreme precipitation gradient from 6,100 mm/year at Paso Galería station, located in Gran Campo Nevado (52º45’S; 73º01’W), to 590 mm/year in Punta Arenas (53º08’S; 70º53’W). In other words, in a west-east transect of only 150 kilometres, annual precipitation decreases by more than one order of magnitude [47].

Snow cover studies and field stations in southern Patagonia are few. A significant decrease of 19% in snow cover extent over an area of 1,637 km2 covering a large part of Brunswick Peninsula (~53ºS), Patagonia, has been detected in the last 45 years (1972-2016 [3] ), attributed to a significant increase in winter temperatures, estimated as 0.71ºC at the weather station of Punta Arenas (Figure 1).

Due to the sparse observational network and uncertainties in the remote sensing snow products, current efforts and studies do not allow to obtain a quantitative and comprehensive assessment of snow cover variability in this region. Therefore, the main objective of this work is to seek a comprehensive and robust method for spatio-temporal reconstruction of snow cover for southern Patagonia that may be replicated in other areas for comparative studies, establishing baselines for future snow and hydrologic modelling in this remote region. We address the following questions: (1) What is the best feasible method for improving the spatial resolution of the MODIS product? (2) What snow detection methods based on remote sensing are most adequate for southernmost Patagonia? (3) How can a daily spatial-temporal reconstruction of snow cover properties be implemented considering a large number of cloudy days in Patagonia? (4) How has the spatio-temporal snow cover variability evolved since 2000?

2. Study Area

Brunswick Peninsula is located in southwest Patagonia (53ºS, 71ºW). Vegetation includes deciduous forest (Nothofagus pumilio and Nothofagus antarctica), peatland (mosses of the genus Sphagnum), grassland and shrubs (i.e., Lepidophyllum cupressiforme, Empetrum rubrum and Gaultheria mucronata). Climate has a strong oceanic influence, being dominated by westerly circulation year-round. The city of Punta Arenas, located in the north-eastern sector of the Peninsula, concentrates 80% of the population of the Magallanes province, one of the areas with the lowest population density in Chile (1.1 inhabitants/km²)3. The water balance of Brunswick Peninsula is critically important regarding the water supply to the population. In addition, the presence of seasonal snow has allowed the development of the ski resort “Club Andino” since 1938, located at Cerro Mirador, at an elevation between 350-625 m, 8 km from Punta Arenas city centre. Snow cover has shown a relevant loss over the last several decades, mainly due to atmospheric warming (Aguirre et al., 2018), significantly reducing the operation of the ski centre, which has only managed to open less than 1 month during the 2016-2019 period (Table 1), whereas until the 1990s typically the skiing season lasted for over 2 months (Rodrigo Adaros, personal communication). The lack of snow at Club Andino has led the regional government to explore an alternative location at Tres Morros hill, 31 km to the southwest of Punta Arenas city, at an elevation between 550-820 m, as a possible alternative site for the development of a future skiing area.

The study area within Brunswick Peninsula is defined by 4 basins, each of which has a river runoff station at the outlet, operated by Dirección General de Aguas (DGA, General Water Directorate of Chile). The basin boundary polygons were obtained from the CAMELS-CL dataset [48] (Figure 1).

3. Data and Methods

In this section, we describe the different data and methods used to reconstruct the spatio-temporal variability of snow in southernmost Patagonia. Data include MODIS sensor Terra satellite products and meteorological records of automatic weather stations, as well as atmospheric reanalysis and regional climate model outputs. All calculations were made using the open-source programming language R v._3.6.3 (https://www.r-project.org)

3.1. Data

a)

Satellite sensor data

The basic information for snow cover reconstruction is obtained from Moderate Resolution Imaging Spectroradiometer (MODIS) sensor data at a resolution of 500 m (3 visible (RGB), 2 near-infrared, and 2 mid-infrared bands) and 250 m (1 visible (red) and 1 near-infrared bands) of MOD09GA and MOD09QA land products version 6 (including atmospherically corrected reflectance), at daily temporal resolution (Figure 2). In addition, MODIS products MOD35_L2 (daily cloud mask at 1 km resolution) and MOD44W versio 6 (water/land mask at a resolution of 250 m), available from the Land Processes Distributed Active Archive Centre (http://lpdaac.usgs.gov/) were used. MOD44W data were downloaded from https://lpdaac.usgs.gov/products/mod44wv006/ and exported to a Geotiff format with the software HEG:HDF-EOS v.2.154, adjusting the area limits for each watershed and considering the nearest neighbour resampling interpolation. Following previous analysis related to the snow season duration on Brunswick Peninsula, which showed snow cover from May to October [3], we downloaded images between April and October to reconstruct the cold-season snow variability. The topographic elevation used in this study is derived from the NASA Shuttle Radar Topographic Mission Global [49] 1 arc second (ca. 30 m) (SRTMGL1) version.0035 digital elevation model (DEM) dataset.

b)

Weather stations

In southernmost Patagonia, the scarcity of locally observed meteorological data in both spatial and temporal representation is a major difficulty for climate change studies. Here, we use data from Tres Morros automatic weather station (AWS), Brunswick Peninsula, deployed on July 27, 2018, specifically for snow cover studies including winter sports evaluation. This station, called AWS Tres Morros, is located in the top of Tres Brazos river watershed (53.3174°S, 71.2790°W) at an elevation of 554 m a.s.l. The variables measured at AWS Tres Morros are: snow height (m), atmospheric pressure (mb), air temperature (ºC), air relative humidity (%), wind direction (˚), velocity (m/s), and short wave and long wave radiation (Wm^-2), both incoming and outgoing.

3.2. Methods

3.2.1. Downscaling and Spectral Fusion

Multimodal remote-sensing image and data fusion is a relevant research field with growing development in the last few years [50]. The goal is to obtain the most complete and accurate information with an adequate spectral and spatial resolution of the study area, given that remote sensing data normally cannot yield both high spectral and high spatial resolution at the same time. Several methods have been developed to combine the advantages of spectral and spatial resolution in one image [51], and in particular regarding MODIS products including Principal Component Analysis (PCA) [52], wavelet transformation [53,54], High-Pass Filter (HPF) [55] and Kriging with External Drift [56]. Of particular relevance is the work by Wang et al. (2015) [57] on downscaling-spectral fusion which improves the MOD09GA at 500 m to higher quality in downscaled images with 250 m resolution using a combination of linear regression and area-to-point kriging interpolation (ATPRK) applied to the residuals of the regression (Eq. 2). The method proposed by Wang et al. (2015) [57] was implemented in this work since it shows that high-quality sharpened images can be obtained, as compared to other commonly-used methods mentioned previously. It is important to consider that these comparisons were made on a single image, something common among comparisons on spectral fusions, without clouds and without snow.

The downscaling-spectral fusion equation of Wang et al. (2015) [57] is composed of a sum of 2 terms, as given by Eq. 2.

$$Z_v^l=Z_{v1}^l+Z_{v2}^l,$$

(2)

where the first term $Z_{v1}^l$ is the result of prediction of regression (presented in Eq. 3) and $Z_{v2}^l$ are the residuals between the regression predictions and the real values for each pixel, using MODIS coarse bands at 500 m resolution (l = [band 3, band 7]).

The prediction of regression

$$Z_{v1}^l=a_l{Z}_v^{lk}+b_l$$

(3)

is a linear regression model where $a_l$ and $b_l$ are the two coefficients for the l-th band and $Z_v^{lk}$ is the covariate band used (band 1 and/or band 2).

Both terms are downscaled to 250 meters. For the Eq. 3, the assumption that this relationship is universal at different spatial resolutions is followed, therefore the relationship obtained for coarse spatial resolutions can be applied to the highest spatial resolution [58]. Therefore, we can use $a_l$ and $b_l$ to predict $Z_{v1}^l$ at 250 m resolution. For the second term of Eq. 2, we use a geostatistical approach to interpolate the residuals from a resolution of 500 m to 250 m.

The geostatistics-based spatial scaling in area-to-point regression (ATPRK) and area-to-area regression kriging (ATARK) have been widely used in remote sensing applications, including downscaling process [57,59,60]. The most important difference between area-to-point and area-to-area kriging is the covariance calculation between samples where area-to-area kriging has the advantage of yielding areal values directly [59]. A recent work of Jin et al. (2018) [61] proposes to modify the linear regression with a geographically-weighted regression (GWR) as a more suitable model that considers the spatial heterogeneity and non-stationarity (Eq. 4), changing the first term in the interpolation approach as follows:

$${\gamma }_i={\beta }_{i0}+\sum \left({\beta }_{ik}{x}_{ik}+{ϵ}_i\right)$$

(4)

In this geographically-weighted linear regression of Eq. 4, each pixel (i) has specific regression coefficients (${\beta }_{i0}$ and ${\beta }_{ik}$) and a covariate band ($x_{ik}$ and errors (${ϵ}_i$).

In this work we will compare both downscaling methods, the GWR (with and without kriging residuals interpolation) and ATARK, exploring Ordinary Least Squares (OLS) and also incorporating a Generalized Least Squares (GLS) method to deal with the problem of violation of homogeneity (also called heteroscedasticy) in the case it is detected [62]. To preserve the spatial correlation (variogram or semivariogram) we fit a spatial correlation model (variogram model). The Akaike criterion (AIC) is used for model evaluation along with the coefficient of the determination (R²) and Root Mean Square Error (RMSE). In contrast to the work of Wang et al. (2015) [57], here we compered the different approaches on an image with snow and other without snow since this study focuses on improving the representation of snow cover and its spectral interaction with the other land covers.

To implement these downscaling schemes, we use the following R packages: raster [63] to manipulate all raster variables and transformation; nlme [64] to implement and analyse linear and nonlinear mixed models; gstat [65] to implement spatial interpolation and analysis; atakrig [59] to implement ATARK; and GWmodel [66] to perform the GWR model. All values out of valid range ([-100, 16000]) were filtered and a water mask at 250-meter spatial resolution was used (MOD44W).

We use different indexes for the quality assessment and selection of the downscaling methodology. To check the quality of the spectral properties for each downscaled band, we use the Universal Quality Index (UQI), a simple and intuitive measure of spectral fidelity between two bands [51,67], with each image having the same dimension (Eq. 5):

$$UQI=4\frac{{\sigma }_{fg}\overline{p}\overline{g}}{\left({\overline{f}}^2+{\overline{g}}^2\right)\left({\sigma }_f^2+{\sigma }_g^2\right)}$$

(5)

In the UQI equation (Eq. 5), $\overline{f}$($\overline{g}$) and ${\sigma }_{\overline{f}}^2$ (${\sigma }_{\overline{g}}^2$) are the mean and variance of band $f$ ($g$) and ${\sigma }_{fg}$ is the covariance between the two bands f and g. This index must be implemented in a sliding window and the average index is then calculated over all these windows. For its application, we assume that the downscaled image at high resolution is incorporated and added to the original low-resolution image [57]. For this purpose, we apply a low pass filter and compare the downscaled image to the original low-resolution image.

Although the above scheme allows for assessing the quality of each downscaled image, additionally we add the ‘Quality Index without Reference’ (QNR, Eq. 6), an index that does not require a reference image [68]. QNR is composed of UQI, Spectral Distortion Index (Dγ, Eq. 7) and Spatial Distortion Index (Ds, Eq. 8).

$$QNR={\left(1-D_{\gamma }\right)}^{\alpha }{\left(1-D_s\right)}^{\beta }$$

(6)

where α and β exponents attribute the relevance of spectral and spatial distortions to the overall quality. The two exponents jointly determine the non-linear response and have values between 0 and 1. The Spectral Distortion Index (Dγ, Eq. 7) is defined as:

$$D_{\gamma }=\sqrt[p]{\frac{1}{L\left(L-1\right)}{{\displaystyle \sum }}_{l=1}^L{{\displaystyle \sum }}_{r=1}^L{\left|Q\left(\widehat{G_l},\widehat{G_r}\right)-Q\left(\widetilde{G_l},\widetilde{G_r}\right)\right|}^p}$$

(7)

where L is the number of bands to be compared (5 bands), $\widehat{G_l}$ and $\widetilde{G_l}$ are the fused and low-resolution version of band l, respectively. Exponent p is a positive integer chosen to emphasize large spectral differences. For p = 1 (as in this study), all differences are equally weighted.

The Spatial Distortion Index ($D_s$) is defined as:

$$Ds=\sqrt[q]{\frac{1}{L}{{\displaystyle \sum }}_{l=1}^L{\left|Q\left(\widehat{G_l},P\right)-Q\left(\widetilde{G_l},\tilde{P}\right)\right|}^q}$$

(8)

where $P$ is the high resolution image used as a reference (band 1 at 250 m resolution in our case), and $\tilde{P}$ a spatially degraded version of the $P$ image obtained by filtering with a lowpass filter. Analogously, $D_s$ is proportional to the q-norm of the difference vector, where q may be chosen to emphasize larger differences. Following the criteria used in Eq. 7 (p=1), q=1 is used to keep all differences equally weighted.

3.2.2 Spectral Mixture Analysis

The spectral mixture analysis (SMA), also known as spectral unmixing, follows the assumption that the spectral intensity in an image (o better said a pixel’s spectral response) is caused by mixtures of a limited number of surface materials, or land classes in our case. Moreover, those weighted “endmembers” fractions (pure spectral signature of a specific land cover measured in the laboratory, in the field, or from the image itself) can be found for each pixel [9,69]. For the SMA implementation, we use the Multiple Endmember Spectral Mixture Analysis (MESMA), an advanced method of linear SMA (Eq. 9) that calculates the component fractions within a pixel allowing endmembers to vary in type and number on a per-pixel basis [70]. For its implementation, the RStoolbox package [71] is used.

In the linear SMA, $\widehat{P_{i\lambda }}$ (Eq. 9 below) is the spectral mixture in pixel i of wavelength λ, modelled as the sum of N reference of endmembers, $P_{k\lambda }$, each weighted by the fraction $f_{ki}$, and the sum of these fraction must be less or equal to 1 to each pixel. ${\epsilon }_{i\lambda }$ is the residual error at λ for the fit of N endmembers, and the respective $RMSE$ (Eq. 10) is used in the optimization process to find each weighted fraction $f_{ki}$, for each endmember.

$$\widehat{P_{i\lambda }}=\sum \left(f_{ki}\ast P_{k\lambda }\right)+{\epsilon }_{i\lambda }$$

(9)

$$RMSE=\sqrt{\frac{\sum {\left({ϵ}_{\lambda }\right)}^2}{B}}$$

(10)

with B the number of bands used.

For the endmember selection of a multispectral image, this methodology allows to relate the number of endmembers to the number of image bands plus 1 (number of bands + 1), thus allowing that the set of equations can be solved simultaneously to produce an exact solution without any error term [9]. Given that we use the first 7 MODIS bands (previously downscaled) in this study, selecting a maximum of 8 endmembers for the algorithm is required. For this, we follow the strategy of endmember selection proposed by Painter et al. (2009) [12], where four endmembers were used for different kinds of snow (based on its grain size), rock, forest, grassland, and peatland. The ability to recognize different types of snow is given by the sensitivity of the spectral reflectance of snow (Figure 2) to its grain size [72], being also possible to estimate its broadband albedo ($\alpha$) using the equation presented below:

$$\alpha =1-A\left({\theta }_0\right)r^{B\left({\theta }_0\right)}$$

(11)

$\alpha$ uses a relationship between the illumination angle (defined as solar zenith angle – terrain slope) ${\theta }_0$ and the grain size of the snow detected, $r$ in the SMA for each pixel. The constants A and B depend on the wavelength of bands and the illumination angle [12]. We use the coefficients of all solar wavelengths, since in the SMA method all MODIS bands are considered: A: 0.0765, 0.0648, and B: 0.2205, 0.2258 for 30º and 60º illumination angles, respectively for short-wavelength frequencies.

For the vegetation endmember representation (forest, grassland, and peatland), we use a seasonal endmember of vegetation classes for each study site. This avoids possible misclassification due to seasonal changes in spectral properties which depend on plant phenology [73]. For this, we use different MODIS images as representatives of each season. We select endmembers through a supervised classification using the Semi-Automatic Classification Plugin (SCP) developed by Congedo (2016) [74] and supported in QGIS, v 3.18 [75].

3.2.3. Spatio-Temporal Snow Reconstruction

While a significant improvement of spatial resolution is necessary for spatial-temporal studies of snow variability in complex terrain, several confounding factors (i.e., the presence of clouds, sensor noise and miscellaneous artifacts) make it difficult to use the snow-cover algorithms directly. This problem of misclassification can be improved implementing a space-time data cube approach as a second step over the first results (SMA) that enables a better estimate of snow detection [76,77]. For this we follow the steps presented below.

a)

Cloud and snow masks

The problems of differentiating snow from clouds using global MODIS products are widely known [15,19,20,24]. Nonetheless, these problems still persist in the current Collection v6 [16]. Since the methodology for improving the cloud mask proposed by Dozier et al. (2008) [76] was not effective in the study region, we explore an alternative, called snow mask, based on the fSC > 0.2 detected by the previous SMA. The snow mask was validated through a combination of NDSI (Eq. 1) and Melt Area Detection Index (MADI, Eq. 12). The underlying principle is based on identifying liquid water layers that coat snow and ice in order to discriminate between frozen and liquid water [78]. The threshold proposed is MADI ≥ 6, which also has been used in studies of snow detection [3,79].

$$\mathrm{MADI}=\frac{\rho B_1}{\rho B_7}$$

(12)

This MADI expression uses the reflectance ratio between the red portion of the visible spectrum ($\rho B_1$, bandwith 620-670 nm, MODIS band 1) and the mid-infrared ($\rho B_7$, bandwith 2105-2155 nm, MODIS band 7).

b)

Temporal interpolation

Considering the need to analyze data in both space and time [80], we apply a gap-filling technique after resolving uncertainties between cloud and snow. This technique consists of a time interpolation on each pixel based on the methodology described by Redpath et al. (2019) [77]. For this, we use a time window of 15 days (i.e., the results are saved on the 7th day) with a maximum of consecutive values of 5 days to be filled. For temporal interpolation in each pixel, we use a smoothing cubic spline interpolation [81] that also is used in the work of Dozier et al. (2008) [76]. The fractional snow cover and broadband albedo interpolations are implemented separately using the snow mask previously defined. Interpolated values out of range were changed by its limits (0 and 1). The gap-filled time series of snow-properties provides the basis for further analysis of snow-cover climatology and variability. For the implementation, we use both raster [63] and space-time [82] R packages.

The scripts (5) and MODIS products referred to in the previous sections are attached as part of the supplementary material.

3.2.4. Ground Validation

To validate the spatial-temporal reconstruction of snow properties, we use daily values of albedo and snow height measurements at AWS Tres Morros averaged over three hours at 11:00, 12:00 and 13:00 local time, considering that MODIS Terra passes over the study area at around 11:00 h local time.

As for the quality assessment of meteorological data, there are conflicting criteria for its implementation. On the one hand, some authors indicate that "erroneous" data can alter the analysis of the time series. On the other hand, other authors believe that the implementation of quality control methods can exclude valid values [83], often dismissing extreme values necessary in climate analysis [84]. In this work, we only review the existence of anomalous values, but no additional statistical tests were implemented for data validation or confidence ranges regarding the exclusion of extreme values.

To validate the proposed methodology, we first compared the snow height (aws_sh) and albedo (aws_al) measurements from AWS Tres Morros and the snow properties reconstruction, based on the steps described above, for the respective pixel and for the days where both coincide. The values considered for comparison are the percentage (%) of snow detection (more than 5% in fSC) based on AWS measurements, considering days with snow height more than 5 centimetres based on the vegetation height below the snow sensor as in the description realized by Rodell and Houser (2004) [85].

We use Pearson correlation and RMSE, exploring linear and non linear relationships between the variables, For the non linear case we use the Generalized Additive Models (GAM) approach using cubic regression spline as a smoothing function [86]. GAM allow for non-linear relationships between the response variable and multiple explanatory variables, also called additive models, fitting a smoothing curve through the data [62,86]. The AIC and deviance explained (as a measure of the quality of fit GAM models and for its Gaussian family, as in our case, is equal to the variance explained) are also used for model evaluation.

3.2.5. Snow Cover Variability

A minimum threshold of fSC > 20% is used to consider a pixel as snow-covered, which is different from the 50% proposed by Redpath et al., (2019) [77] due to the high interaction between the snow and vegetation in the study area. In the pixel corresponding to AWS Tres Morros we have a combination of peatland, Nothofagus forest, shrubs and grasses. Daily and monthly temporal variability of the snow-covered area and the percentage of the study area (referred to the 4 watersheds) covered by snow is used to determine the seasonal snow duration in the cold season (April-October). The number of snow days are also considered a spatial indicator of snow persistence at Brunswick Peninsula.

To study snow-cover variability and trend analysis, we use Mann Kendall [87] and Sen’s slope [88] non-parametric methods, widely applied in climate trends analysis [89]. Before applying these trend analyses, we implement an exponential filter [90], specifically in our case the Holt-Winter filter, which considers both seasonal and non-seasonal behaviour that depends on the time series analysed, and also minimizing the square prediction error [91]. We also use the time series decomposition additive approach [92]. To complement the temporal analysis that considers the overall trends over the entire study area, we calculate trends of each pixel considering the variation of the number of snow days in the time series of annual snow days. With both products, covering spatial and temporal variability, we explore the relationships that could explain the spatial distribution of the trends detected. For this, we use Principal Component Analysis (PCA) and further linear and non-linear relationships with spatial variables as latitude, longitude, exposure orientation (aspect) and height. The R packages used are wql [93], stats [91] and mgcv [86].