Prerpints.org logo
Preprint
Article

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

Altmetrics

Downloads

174

Views

69

Comments

0

A peer-reviewed article of this preprint also exists.

Submitted:

09 August 2023

Posted:

14 August 2023

You are already at the latest version

Alerts
Abstract
Seasonal snow cover is a fundamental component of both the global energy budget and the water cycle. Its properties such as fractional snow cover or albedo are particularly affected by climate change. Several methods based on satellite data products are available to estimate these properties, each one with its pros and cons. This work presents a novel methodology that integrates three indexes applied to MODIS satellite data (Spectral Mixture Analysis (SMA), Normalized Difference Snow Index (NDSI) and Melt Area Detection Index (MADI)), to perform a spatio-temporal reconstruction of the fractional snow cover and albedo at 250 m spatial resolution in the Brunswick Peninsula, southwest Patagonia during the cold season (April-October) for the period 2000-2020. Three main steps are included: (1) the increase of the spatial resolution of MODIS (MOD09) data to 250 m using a spectral fusion technique; (2) the snow-cloud discrimination; (3) the daily spatio-temporal reconstruction of snow extent and its albedo signature with subpixel detection using the endmembers extraction and spectral mixture analysis. The results show a 98% agreement between MODIS snow detection and ground-based snow measurement at the automatic weather station Tres Morros (53.3174 °S, 71.2790 °W), with fractional snow cover values between 20% and 50%, showing a close relationship between snow and vegetation type. The number of snow days compiled from the MODIS data indicates a good performance (Pearson correlation of 0.9) compared with the number of skiing days at Cerro Mirador ski centre near Punta Arenas. Although the number of seasonal snow days show a significant increase trend of 0.54 days/year in Brunswick Peninsula during the 2000-2020 period a significant decreasing trend of -4.64 days/year was detected during 2010-2020, and also below the 400 m a.s.l. elevation, which is the most affected area. A reconstruction of the monthly mean temperature using the ERA5 Land reanalysis product shows a significant warming trend in May (0.068 ºC/year) and October (0.098 ºC/year) in the 2000-2020 period. Under the future emission scenario RCP8.5, the regional climate model RegCM4 predicts further warming during 2020-2050 of 0.059 ºC/year in July, 0.088 ºC/year in August, and 0.019 ºC/year in October, which will further reduce snow cover.
Keywords: 
Subject: Environmental and Earth Sciences  -   Remote Sensing

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].
N D S I = ρ B 4 − ρ B 6 ρ B 4 + ρ B 6
where ρ B 4 is the green band reflectance (bandwidth 545–565 nm), and ρ 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/km2)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 v 1 l + Z v 2 l ,
where the first term Z v 1 l is the result of prediction of regression (presented in Eq. 3) and Z v 2 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 v 1 l = a l Z v l k + b l
is a linear regression model where a l and b l are the two coefficients for the l-th band and Z v l k 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 v 1 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:
γ i = β i 0 + ∑ β i k x i k + ϵ i
In this geographically-weighted linear regression of Eq. 4, each pixel (i) has specific regression coefficients ( β i 0 and β i k ) and a covariate band ( x i k 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 (R2) 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):
U Q I = 4 σ f g p ¯ g ¯ f ¯ 2 + g ¯ 2 σ f 2 + σ g 2
In the UQI equation (Eq. 5), f ¯ ( g ¯ ) and σ f ¯ 2 ( σ g ¯ 2 ) are the mean and variance of band f ( g ) and σ f g 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).
Q N R = 1 − D γ α 1 − D s β
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 γ = 1 L L − 1 ∑ l = 1 L ∑ r = 1 L Q G l ^ , G r ^ − Q G l ˜ , G r ˜ p p
where L is the number of bands to be compared (5 bands), G l ^ and 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:
D s = 1 L ∑ l = 1 L Q G l ^ , P − Q G l ˜ , P ˜ q q
where P is the high resolution image used as a reference (band 1 at 250 m resolution in our case), and 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, P i λ ^ (Eq. 9 below) is the spectral mixture in pixel i of wavelength λ, modelled as the sum of N reference of endmembers, P k λ , each weighted by the fraction f k i , and the sum of these fraction must be less or equal to 1 to each pixel. ε i λ is the residual error at λ for the fit of N endmembers, and the respective R M S E (Eq. 10) is used in the optimization process to find each weighted fraction f k i , for each endmember.
P i λ ^ = ∑ f k i ∗ P k λ + ε i λ
R M S E = ∑ ϵ λ 2 B
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 ( α ) using the equation presented below:
α = 1 − A θ 0 r B θ 0
α uses a relationship between the illumination angle (defined as solar zenith angle – terrain slope) θ 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].
MADI = ρ B 1 ρ B 7
This MADI expression uses the reflectance ratio between the red portion of the visible spectrum ( ρ B 1 , bandwith 620-670 nm, MODIS band 1) and the mid-infrared ( ρ 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].

4. Results

4.1. Spectral Fusion

The comparative study between different spectral fusion methods follows the methodology described by Wang et al. (2015) [57]. Linear regression results, following (Eq. 3), were reviewed using both snow and snowless scenes over a selected area on Brunswick Peninsula. Data from MODIS band 5 were also included, not initially considered in Wang et al. (2015) [57]. Table 4 shows a good linear fit in snowless areas between band 1 (as an explanatory variable) and bands 3, 4, 5, 6, and 7, as also indicated by Wang et al. (2015) [57]. However, it is necessary to re-evaluate this relationship in snow-covered areas since the linear fit is less convincing over snow condition.

4.1.1. The First-Term of Spectral Fusion, the Linear Relationship

As a first step, and given the different results between snow and snowless images, using the linear relationship proposed by Wang et al. (2015), coarse band i ~ band 1 (band i in range 3 to 7), we explored other relationships and found better results with the relationships detailed in Table 5.
An exploratory analysis of the different MODIS bands over a reference MODIS image of September 25, 2015, shows the need for normalization by logarithmic transformation since the reflectance values in each band have an exponential distribution [62]. After this normalization, as a first stage, an Ordinary Least Squares (OLS) adjustment is implemented using as an independent variable the MODIS band 1 and the resampled DEM at 500 meters of spatial resolution (DEM500) as indicated in Table 4. The distribution of the residuals shows a clear violation of heterogeneity, justifying the recommendation of Wang et al. (2015) [57] to use a GLS instead of an OLS model. Additionally, the existence of a spatial correlation was reviewed with a variogram analysis of its residuals (Figure 4), with a clear spatial dependence that must be considered.
Using the new relationships shown in Table 5, a GLS model was implemented for the coarse bands (GLS500), with a selection of variance structures for each band, according to the protocol as indicated by Zuur et al., (2009) [62]. For the example image (September 25, 2015), the following variance structures were selected: for Band_3500 power of the covariate variance structure; for Band_4500 exponential variance structure; for Band_5500, Band_6500 and Band_7500 constant plus power of the variance covariate function.
The GLS model performed better than the OLS model according to the AIC, but not so compared to the results obtained with a GWR500. The GWR model has the advantage of generating coefficients for each pixel. The coefficients of determination (R2) between the modelled coarse bands using GWR500 and the original coarse bands, varied between 0.97 and 0.99 as follows: Band_3500 = 0.99, Band_4500 = 0.99, Band_5500 = 0.97, Band_6500 = 0.98, and Band_7500 = 0.97, with low RMSE values in all models (Table 6).
Finally, to obtain the first term of the spectral fusion ( Z v 1 l ) of Eq. 2, we need to use these relationships (at coarse resolution) and implement them at a higher resolution (250 meters). For the cases of OLS500 and GLS500, the key issue is the estimation of the regression coefficients for the coarse bands, assumed to be universal at different spatial resolutions. Consequently, the relationship built at coarse spatial resolutions can be applied at a higher spatial resolution [58]. This assumption is the same for the GWR500 model, but since there are coefficients for each pixel (β0 and βi in Eq. 4), these are assigned to the highest resolution imagery on a pixel by pixel basis of the original data as an additional information layer.

4.1.2. The Second Term of Spectral Fusion, Kriging Interpolation

The next step is to apply the ATAK interpolation to residuals of the coarse regression (500 meters) to each model (OLS500, GLS500 and GWR500), so as to redistribute the error at a resolution of 250 m. First, we need to select the best model representing the residuals’ spatial correlation (variogram). Once this model has been estimated for each band, we proceed with implementing the ATAK interpolation to obtain the second term of Eq. 2 ( Z v 2 l ). This comparative process is detailed in Script Nº1 of Appendix A (Supplementary data and codes), and its results are shown in Table 6.
While in the coarse resolution data, the results of the GWS500 model are the best option, once implemented in the fine resolution data (250 m), the three ATAK models achieve a similar behaviour based on the QNR Index. The best results correspond to the GLS_ATAK250 model, with a QNR Index of 0.949. However, regarding the results for each band’s UQI, the GLS model is the best (Table 6). The QNR Index represents the performance of each spectral fusion model, where in our case, we consider the spatial (Ds) and spectral (Dγ) distortion with equal relevance. An experiment with the GWR250 model without ATAK interpolation was performed, obtaining a QNR index of 0.914, which improved to 0.946 after applying the ATAK interpolation (GWR_ATAK250). Obviously, the improvement is not so significant compared to the OLS and GLS models. This can be explained in part by the spatial correlation of its residuals (Eq. 2), which lacks a good model representation. On the other hand, the computation time of the GLS_ATAK250 method is almost three times faster than the GWS_ATAK250, which is preferable given that one of our goals is to have an operational methodology for a large amount of data to be processed.
All algorithms are contained in Script Nº2 (Appendix A) which includes: (i) the selection of the best variance structure of the model for each band, (ii) the selection of the best variogram model of the residuals, and (iii) controls on the possible issues in this downscaling batch process, managing missing data, outliers, image data problems (We include a threshold for detecting a lack of spectral information (for saturation or other spectral problems) of a standard deviation of the digital numbers lower than 10 considering the whole picture, based on our experiments). Figure 5 shows the residual maps of the GLS regression at 250 m and the sum of these two products (Eq. 2), compared with the corresponding coarse band.

4.2. SMA

Before implementing the SMA, as is mentioned in the methodology, it is necessary to find the seasonal spectral signatures (endmembers) of the principal vegetation types present at Brunswick Peninsula. In an explorative analysis of the spectral signature and its representation as endmembers at Brunswick Peninsula we selected forest, grassland and peatland, in addition to the 4 categories of snow and bare ground. Table 7 shows the days selected to detect the seasonal spectral signature for each vegetation type.
Figure 6 shows how each spectral signature varies seasonally. Forest presents a rather low reflectance compared to grassland and peatland. This seasonal variability reaffirms the importance of considering seasonal changes in the vegetational endmembers for this type of spectral classification.
In Figure 7 the results of the SMA are presented. Figure 7B shows snow fraction (%) and Figure 7D snow grain size (μm), compared with the topographic elevation (Figure 7A). For broadband albedo calculation refer to Script Nº 4 in methodology and Appendix A.

4.3. Spatio-Temporal Snow Reconstruction

4.3.1. Cloud and Snow Masks

After performing the SMA, we contrasted the results of fractional snow cover, grain size, and cloud cover mask over 3 days with different atmospheric conditions: (i) with the presence of both snow and cloud (August 6, 2015), (ii) mostly cloudy (August 8, 2015), and (iii) snow with clear sky (September 25, 2015). The cloud mask confuses snow and cloud cover in all conditions, limiting its use (Figure 8). Considering that the west side of southernmost Patagonia is affected by a high number of cloudy days (77–86%) and precipitation days (81–88%) throughout the year [44], this issue must be considered especially in an area of climatic transition as occurs in the Brunswick Peninsula. The RGB composition of grain size (red), NDSI (green) and MADI (blue), shows a good discrimination of snow and cloud. A snow mask is generated by applying the thresholds indicated in the methods section. Using this snow mask to discriminate the snow pixels (value 1 as snow, value 0 as no snow, and nan as indeterminate), we improve the discrimination between clouds and snow, and apply it to reconstruct snow fraction and broadband albedo. Script N° 4 of Appendix A describes this process.

4.3.2. Temporal Interpolation

Taking advantage of the new snow mask proposal, developed in the previous section, and its implementation based on 21 years of MODIS data in the study area, we made a temporal interpolation in each pixel for fSC and albedo. Figure 9 presents how the cubic spline smoothing function, described in detail in the method section, works as a moving window of 15 days (7 days forward and 7 days backward) and can fill the missing data of fractional snow cover (fSC) and albedo based on past and present behaviour for each pixel due to cloud cover condition. We can note after the implementation of the interpolation technique that several grid points with missing values are filled now, which are particularly notable on August 6 and 8 because they have several problems with cloud interference, as we can see in the previous Figure 8. This process was also implemented for batch processing (Script Nº 5 in the Appendix).

4.4 Ground Validation with AWS Data

Figure 10 shows the annual winter variability of the snow height measured at AWS Tres Morros from July 2018 (installation date) to December 2020, also showing when the snow height exceeds 20 cm for the sampling period. Data were processed and compared with the MODIS reconstructed data corresponding to the 250 m pixel location. For the AWS albedo data at the Tres Morros a minimum snow height of 20 cm was established to consider representativeness, avoiding problems with mixed snow and vegetation signals. The valid albedo data are limited to July (only four days), August, and September of 2018 due to problems in field data, with mean (maximum) values in August of 0.82 (0.96) and September 0.72 (0.87). We only found 24 of 159 valid AWS albedo measurements (with snow height > 20 cm) that coincide on the same day with the reconstructed MODIS albedo at 250 m resolution, exhibiting a poor correlation (Pearson correlation of 0.1).
For fSC reconstruction (SF_Rec) from MODIS imagery for 2018-2020, 135 days of MODIS data at the Tres Morros pixel were compared with AWS Tres Morros snow height measurements. A Pearson correlation of 0.55 was found in the Linear relation (SF_Rec ~ aws_sh) and RMSE of 0.113. However, exploring a GAM model (SF_Rec ~ s(aws_sh)) the RMSE reduces to 0.089 and R2 increases to 0.55. AIC values are -1998 (Linear) and -252.4 (GAM), respectively, indicating that the GAM model provides the preferred approach. Nonetheless, for snow height values below 20 cm (presented in the relation GAM 20 cm) fSC values of zero are found suggesting that the threshold of fSC in the snow mask (fSC > 0.2) for this pixel leaves some fSC undetected. Table 8 shows the results of the AWS Tres Morros data vs the reconstructed values. As for the confidence of snow detection at this pixel, we find a 98% coincidence between the fSC (with a threshold of 10%) and the snow height measurement at AWS Tres Morros, using a threshold of 5 cm.

4.5 Reconstructed Snow Cover Variability at Brunswick Peninsula

Figure 11A shows a box plot of daily snow cover area during the cold season for each year for the entire study area. The highest number of outlier values coincides with the years that present a short snow season (2003, 2014, 2016 and 2017), although during those years several days showed high amounts of snow. The snow season (shown in Figure 11B,C) was defined as a minimum percentage of snow cover area of 10% (on Brunswick Peninsula), which is highly related with the number of skiing days at the ski centre located in Cerro Mirador (Pearson correlation of 0.9, RMSE of 14 days, [3]), showing an average of 68 days (Figure 11C) for the 2000-2020 period. For the last 9 years (2012-2020), 7 years are below this average. The monthly trends analysis (Figure A1 in Supplementary section) of both the 2000-2020 and 2010-2022 periods, using the MK and Sen´s slope, do not show a significant trend (at 95%). However, the monthly behaviour of the time series changed during the 2010-2020 period, with September showing a decreasing trend of 7.5 km2 per year, and after applying a Holt Winter filter it showed a significant (at 95%) decreasing trend of 4.7 km2 per year. The reconstructed snow days (snow season days) for the whole period (2000-2020) show a trend of +0.11 days/year (not significant), but applying an exponential filter it presents a significant trend of +0.54 days/year in Brunswick Peninsula (p-value < 0.005). However, for the last 10 years (2010-2020), a significant decreasing trend of -4.64 snow days/year is observed in Brunswick Peninsula (p-value < 0.001).
Figure 12A identifies the spatial distribution of the snow season duration in days (as an average for 2000-2020) for each pixel, finding 65 and 75 snow days in Cerro Mirador and Tres Morros, respectively. Figure 12B shows the significant trend of yearly snow days (snow days/year) of the whole series, with a trend of -0.5 snow days/year (decrease) and +0.3 snow days/year (increase) in Cerro Mirador and Tres Morros sites, respectively. The yearly snow days dataset (12-A and C) is correlated with terrain information (elevation and aspect) and spatial location (latitude and longitude) but these correlations are not so representative in the yearly snow days trend (12-B and D). The PCA analysis (presented as supplementary information, Figure A2) for the model of average snow days (Mod_av) identified that all terrain and spatial variables have a significant contribution, explaining 54.9% of the variance using the two principal eigenvalues. This is also indicated by its p-value (less than 0.05) of all terms in the linear regression, with the average annual snow season duration data normalized by a logarithmic function given the exponential distribution of its data. For the case of the model of snow-days trend (Mod_trend), all terrain and spatial variables have a significant contribution which explains 56% of the variance based on the two principal eigenvalues, which was also identified in the linear model with a significance (p-value) < 0.05. For Mod_av the elevation variable yielded the highest Pearson correlation (-0.63) for a linear relationship, with a Mod_trend elevation of longitude (long) and latitude (lat) of 0.11, 0.15 and 0.14, respectively. In both cases (Mod_av and Mod_trend) a GAM was implemented using a cubic spline regression as a smoothing function (Wood, 2017; Zuur et al., 2009). Figure 12A (Mod_av) and Figure 12B (Mod_trend) present the elevation component of the respective GAM which shows the highest correlation with the response variables in both cases. Elevation has the most pronounced effect on snow cover (snow days/year) in the study area with a range of 380-500 m a.s.l for yearly snow days. However, in the case of yearly trends of snow days, elevation does not have a strong correlation. The implemented GAM explain 64.7% (for snow days) and 22.7% (for yearly trend) of the total variance, respectively.

5 Discussion

5.1 Method Improvement

The methods presented here include tests of recent downscaling methodologies [57,61] regarding both accuracy and speed of implementation (computing cost) to improve the spatial resolution of wintertime acquisitions, since these procedures have mainly been performed under summer conditions [57]. The method proposed here is evaluated for each band (UQI), as well as using a global analysis (QNR Index) for the preservation of spectral ( D γ ) and spatial ( D s ) characteristics. An improvement is achieved for all indices, and also regarding computational implementation. It is important to consider that our methos was compared on two images, one image with snow and other without snow, in contrast to the work of Wang et al. (2015) [57]. where these comparisons were made on a single and clear image (without clouds and without snow). It would be interesting for future applications of this methodology create a long-time series to provide a more extensive assessment for all the images processed but for our study this is out of the main scope since it is the first time this technique is implemented over this region.
Of high importance are cloud and snow discrimination. After testing different approaches between cloud and snow discrimination (i.e., [76,77,94]), we propose here a new method. The proposed method combines SMA, NDSI and MADI indexes, automatically processing 21 years of MODIS images, obtaining daily snow fSC and albedo values at a spatial resolution of 250 m.
The proposed methodological validation includes: (1) the spectral fusion for downscaling process, (2) the spectral mixture analysis, (3) the snow/cloud discrimination using a snow mask, and (4) the temporal interpolation. These steps were compared with ground measurements at AWS Tres Morros. Our reconstruction shows a coincidence of 98% between the fSC (with a threshold fSC ≥ 10%) and the measurements of snow height (aws_sh) at the location of AWS Tres Morros (with a threshold Snow_h ≥ 5 cm). We also find a relationship (Generalized Additive Model) between SF_Rec ~ s(aws_sh) that explains variance of 55 %. However, this relationship depends on the height of vegetation, so it cannot be extrapolated to other pixels, but may be extrapolated to the past for the pixel tested. Of a total of 159 days of valid field AWS albedo measurements, only a period of 35 days overlap with the satellite reconstructed data, exhibiting a poor correlation (Pearson correlation of 0.1). AWS Tres Morros is located below the tree line, so the pixel used for this validation has a relevant vegetation component, with a fSC of 20% to 50%, and a maximum snow height of 80 cm. It would be necessary to measure albedo in the field at several locations within the 250 x 250 m pixel to obtain an average representation of albedo to obtain a better correlation between field and satellite data.
Regarding snow cover variability, we limit the representation to the behaviour of the entire study area to daily, monthly and yearly resolution referred to the cold season (April-October). The snow days (season duration) have a good performance (Pearson correlation of 0.9) compared with the number of skiing days at Cerro Mirador presented in Aguirre et al. (2018) [3], with a significant negative snow days trend of 0.5 days/year at Cerro Mirador for 2000-2020. For the whole of Brunswick Peninsula the snow days show a significant positive trend of 0.54 snow day/year for the 2000-2020 period, while a strong significant decreasing trend of -4.64 snow day/year was detected for the last 10 years (2010-2020). In a long-term analysis, the work of Aguirre et al., (2018) shows a significant decrease trend of 19% of snow cover in Brunswick Peninsula for the last 45 years (2000-2016) which can be attributed to a statistically significant long-term warming of 0.71°C at Punta Arenas during the winter.

5.2 Climate Forcing

The spatio-temporal model shown in Figure 12, which represents both the yearly snow days trend and snow days as annual average (2000-2020), explains 22.7% and 64.7% of the data deviance, respectively. The most affected areas regarding snow cover reduction are those located within the elevation band below 400 m a.s.l. In an earlier study [3], the primary forcing responsible for this snow-cover reduction in Patagonia was attributed to atmospheric warming. Here we expand and update this climate analysis.
The Southern Annular Mode (SAM), the main atmospheric mode of the Southern Hemisphere, controls the position and strength of the southern westerly winds [33,95], and may be described as the difference between the mean sea level pressure at mid-latitudes (~40°S) and Antarctica (~65°S). Over the last 50 years, the SAM has trended towards its positive phase due to an increased pressure gradient between midlatitudes and Antarctica. This has resulted in a southern migration and enhancement of the westerlies, which may continue according to the projected future climate change scenarios [96,97,98].
Depending on the latitudinal shift of the westerlies compared to the position of the study area in Patagonia, this could imply changes in the precipitation totals. However, a more evident effect of a positive SAM is reinforcing the longitudinal rainfall gradient between the western (higher rainfall) and eastern (lower rainfall) margins of the Andes at Patagonian latitudes due to enhanced westerly air flow (i.e., [99]).
An unequivocal forcing of the snow cover reduction is due to the regional temperature increase in southern Patagonia. Significant warming of 0.72 ºC (0.12 ºC/decade) has been detected in 1958-2016 during the winter season [3], based on data of the Punta Arenas Airport weather station. The work by Srur et al., (2018) [100] in Patagonia (49.5 ºS / 73 ºW) points to a similar observation, where a displacement of the tree line (Nothofagus pumilio) to higher elevations was detected, also pointing to atmospheric warming.
To study the climate in remote environments with poor availability of weather stations [101,102] atmospheric model outputs data can be used, such as the ERA5 reanalysis dataset. This ECMWF reanalysis product covers the period from 1979 to present, at a coarse resolution of 30 km. A higher resolution (ERA5 Land) of 9 km [103], which has recently been used in southern Patagonia (i.e., [35,104,105]) is recently available. Here we analyze ERA5 Land hourly temperature and total precipitation data to extend and compare climate variables with the MODIS snow reconstruction at Brunswick Peninsula. Based on ERA5 data, we build mean, minimum and maximum monthly temperature and monthly degree-hours (defined as the number of monthly hours above 0ºC). For solid and liquid precipitation, we use a threshold relation with air temperature to estimate the solid/liquid fraction in total precipitation as indicated in the work of Schaefer et al. (2013) [106], with an upper-temperature threshold of 1.6 ºC, and a lower threshold of 0.9 ºC, as proposed by Aguirre et al. (2018) [3]. Above 1.6 ºC only liquid precipitation occurs, while below 0.9 ºC only solid precipitation is observed, with a linear change between these two thresholds. Table 7 shows the cross Pearson correlations between monthly snow cover area (km2), reconstructed from MODIS data, and the ERA5 Land climate variables for the cold season (April to October).
Table 7. Cross Pearson correlation between snow cover area from MODIS data at the AWS Tres Morros and ERA5 land monthly variables from ERA5 data at the same location for the period 2000 to 2020, between April to October. ERA 5 variables are liquid and solid precipitation; mean, minimum and maximum temperature; and monthly degree-hours.
Table 7. Cross Pearson correlation between snow cover area from MODIS data at the AWS Tres Morros and ERA5 land monthly variables from ERA5 data at the same location for the period 2000 to 2020, between April to October. ERA 5 variables are liquid and solid precipitation; mean, minimum and maximum temperature; and monthly degree-hours.
ERA5 climate variables Liquid precipitation Solid precipitation Mean temperature Maximum temperature Minimum temperature Degree-hours
Cross Pearson correlation with MODIS snow cover area -0.49 0.48 -0.70 -0.65 -0.60 -0.70
This result confirms that the main climatic driver for the snow-cover reduction is related to the mean atmospheric temperature increase, or degree-hours, both of which show a cross Pearson correlation of -0.70. A linear regression of monthly mean atmospheric temperature and MODIS snow cover area shows a determination coefficient R2 of 0.48 (p<0.001). We then adjusted a linear model and a non-linear GAM model to explain the MODIS snow area based on mean atmospheric temperature, liquid precipitation and solid precipitation. In those two models, only the mean temperature proved to be significant, while the other two variables (liquid and solid precipitation) did not result in a significant relation with snow cover extent.
The regional climate model RegCM4 at a resolution of 10 km has recently been run for all of Chile: hindcast 1979 to 2015; historical 1975 to 2005; modelled near future 2020 to 2050. Two different climate change scenarios have been considered: RCP 2.6 and 8.5 [107]. Figure 13A shows the daily temperature correlation (R2=0.86) between ERA5 Land product and AWS Tres Morros data for 2018-2020. This correlation increases significantly to R2=0.98 on a monthly scale (Figure 13B). With this validation, we use ERA5 land temperature data to check the performance of RegCM4 at AWS Tres Morros (2000-2005), obtaining R2 = 0.53, which a GAM model adequately represents due to its non-linear behaviour (Figure 13C).
The trend in ERA5 Land monthly mean temperature extracted for the location at AWS Tres Morros shows significant warming in May (0.068 ºC/year) and October (0.098 ºC/year) for the 2000-2020 period, which coincide with the start and end of the cold season. For the future (2020-2050) scenario RCP 8.5 at Brunswick Peninsula, we use the regional climate model RegCM4 corrected by ERA5 Land, considering the relationship found in Figure 13C. This scenario presents a significant warming trend during July (0.059 ºC/year), August (0.088 ºC/year) and October (0.019 ºC/year), which predicts that snow cover reduction will continue in Patagonia.

6. Conclusions

This work presents a novel and robust methodology for a high-resolution spatio-temporal snow cover reconstruction in southern Patagonia. Results include snow extent and albedo at subpixel level (250 m), including an improvement of snow-cloud discrimination as well at these latitudes (~53ºS), validated at Brunswick Peninsula,
The results show a significant decreasing trend of the snow season of 4.64 days/year within the period 2010-2020. The most affected elevation is below 400 meters, with a non-significant negative trend about -0.5 snow days per year, mainly due to warming at the beginning and end of the snow season (May and October). The projection under the future RCP8.5 scenario, using results from the RegCM4 regional climate model at a resolution of 10 km, shows significant warming of 0.059 ºC/year during July, 0.088 ºC/year during August, and 0.019 ºC/year during October in the 2020-2050 period, in contrast to the study period (2000-2020) where ERA5 Land shows significant warming in May (0.068 ºC/year) and October (0.098 ºC/year) on Brunswick Peninsula.
Given the results of this work, it is important to improve the quantification of snow cover (i.e., fSC, snow depth, snow water equivalent), especially in vegetated areas, such as in the low-mid elevations in Patagonia, since relevant snow underestimations can occur. Furthermore, due to the relevance of snow cover for human development, water resources and the environment, its evolution under future climate change scenarios should be studied in detail, deploying a robust ground-based monitoring data collection, over a broader territorial representation.

Acknowledgments

This study is co-funded by the FNDR GORE Magallanes “Programa de Transferencia Científico Tecnológico Modelamiento Climático Planificación, XII Región, Código BIP No 30462410”, known as the MoCliM Programme, executed by the Chilean Antarctic Institute, with the collaboration of University of Magallanes. We also appreciate the support and information provided by Corporación Nacional Forestal (CONAF) of Magallanes and Club Andino of Punta Arenas. Collaboration by Inti González and Rodrigo Adaros during field work is appreciated. Intense computational processes, such as the downscaling algorithms, were performed at the National Laboratory of High-Performance Computing (NLHPC) of Chile. We are also grateful for the contributions made by Cape Horn International Centre (CHIC), ANID/BASAL FB210018, and Gaia Antarctica Research Centre (CIGA) of Universidad de Magallanes.

Appendix A. Supplementary data and codes

MODIS Data dawload
To optimise the download process, we used the script “order_MWS.pl”, available from:
https://ladsweb.modaps.eosdis.nasa.gov/tools- and-services/.
Script Nº1: Comparative spectral fusion methods.
Script Nº2: Quality Index for comparative spectral fusion.
Script Nº3: Bach process with spectral fusion for an operative implementation.
Script Nº4: Batch process for snow detection, cloud correction and albedo.
Script Nº5: Spatio-temporal snow cover reconstruction.
Preprints 82050 g0a1
Figure A1. Monthly trend of snow cover area at Brunswick Peninsula.
Figure A1. Monthly trend of snow cover area at Brunswick Peninsula.
Preprints 82050 g0a1
Preprints 82050 g0a2
Figure A2. PCA of Snow days average and Trend snow days.
Figure A2. PCA of Snow days average and Trend snow days.
Preprints 82050 g0a2
Preprints 82050 g0a3
Figure A3. Figure A3. Generalized Additive Models (GAM) to Mod_av.
Figure A3. Figure A3. Generalized Additive Models (GAM) to Mod_av.
Preprints 82050 g0a3
Preprints 82050 g0a4
Figure A4. Figure A4. Generalized Additive Models (GAM) to Mod_trend.
Figure A4. Figure A4. Generalized Additive Models (GAM) to Mod_trend.
Preprints 82050 g0a4

References

  1. Chen, X.; Long, D.; Liang, S.; He, L.; Zeng, C.; Hao, X.; Hong, Y. Developing a Composite Daily Snow Cover Extent Record over the Tibetan Plateau from 1981 to 2016 Using Multisource Data. Remote Sens Environ 2018, 215, 284–299. [Google Scholar] [CrossRef]
  2. Frei, A.; Tedesco, M.; Lee, S.; Foster, J.; Hall, D.K.; Kelly, R.; Robinson, D.A. A Review of Global Satellite-Derived Snow Products. Advances in Space Research 2012, 50, 1007–1029. [Google Scholar] [CrossRef]
  3. Aguirre, F.; Carrasco, J.; Sauter, T.; Schneider, C.; Gaete, K.; Garín, E.; Adaros, R.; Butorovic, N.; Jaña, R.; Casassa, G. Snow Cover Change as a Climate Indicator in Brunswick Peninsula, Patagonia. Front Earth Sci (Lausanne) 2018, 6, 130. [Google Scholar] [CrossRef]
  4. Hughes, L. Climate Change and Australia: Trends, Projections and Impacts. Austral Ecol 2003, 28, 423–443. [Google Scholar] [CrossRef]
  5. Beniston, M.; Farinotti, D.; Stoffel, M.; Andreassen, L.M.; Coppola, E.; Eckert, N.; Fantini, A.; Giacona, F.; Hauck, C.; Huss, M.; et al. The European Mountain Cryosphere: A Review of Its Current State, Trends, and Future Challenges. Cryosphere 2018, 12, 759–794. [Google Scholar] [CrossRef]
  6. IPCC Summary for Policymakers. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, 2013; 33. [CrossRef]
  7. Mernild, S.H.; Liston, G.E.; Hiemstra, C.A.; Malmros, J.K.; Yde, J.C.; McPhee, J. The Andes Cordillera. Part I: Snow Distribution, Properties, and Trends (1979-2014). International Journal of Climatology 2017, 37, 1680–1698. [Google Scholar] [CrossRef]
  8. Walker, D.A.; Halfpenny, J.C.; Walker, M.D.; Wessman, C.A. Long-Term Studies of Snow-Vegetation Interactions. Bioscience 1993, 43, 287–301. [Google Scholar] [CrossRef]
  9. Lillesand, T.M.; Kiefer, R.W.; Chipman, J.W. Remote Sensing and Image Interpretation; Seventh.; John Wiley & Sons, 2015. ISBN 9781118343289.
  10. Tedesco, M. Electromagnetic Properties of Components of the Cryosphere. In Remote sensing of the cryosphere; Tedesco, M., Knight, P., Eds.; Wiley-Blackwell: Oxford, UK, 2015; pp. 17–29. ISBN 9781118368862. [Google Scholar]
  11. Cortés, G.; Girotto, M.; Margulis, S.A. Analysis of Sub-Pixel Snow and Ice Extent over the Extratropical Andes Using Spectral Unmixing of Historical Landsat Imagery. Remote Sens Environ 2014, 141, 64–78. [Google Scholar] [CrossRef]
  12. Painter, T.H.; Rittger, K.; McKenzie, C.; Slaughter, P.; Davis, R.E.; Dozier, J. Retrieval of Subpixel Snow Covered Area, Grain Size, and Albedo from MODIS. Remote Sens Environ 2009, 113, 868–879. [Google Scholar] [CrossRef]
  13. Rittger, K.; Painter, T.H.; Dozier, J. Assessment of Methods for Mapping Snow Cover from MODIS. Adv Water Resour 2013, 51, 367–380. [Google Scholar] [CrossRef]
  14. Hall, D.K.; Frei, A.; Dery, S. Remote Sensing of Snow Extent. In Remote sensing of the cryosphere; Tedesco, M., Knight, P., Eds.; Wiley-Blackwell: Oxford, UK, 2014; pp. 31–47. ISBN 9781118368862. [Google Scholar]
  15. Hall, D.K.; Riggs, G.A. Accuracy Assessment of the MODIS Snow Products. Hydrol Process 2007, 21, 1534–1547. [Google Scholar] [CrossRef]
  16. Riggs, G.; Hall, D. MODIS Snow Products User Guide to Collection 6. 2017.
  17. Hall, D.K.; Riggs, G.A. Normalized-Difference Snow Index (NDSI). In Encyclopedia of snow, ice and glaciers; Singh, V.P. (Vijay P.), Singh, P. (Hydrologist), Haritashya, U.K., Eds.; Springer, Dordrecht: Zurich, Switzerland, 2011; pp. 779–780. [Google Scholar]
  18. Frei, A.; Lee, S. A Comparison of Optical-Band Based Snow Extent Products during Spring over North America. Remote Sens Environ 2010, 114, 1940–1948. [Google Scholar] [CrossRef]
  19. Hall, D.K.; Riggs, G.A.; Salomonson, V. V; DiGirolamo, N.E.; Bayr, K.J. MODIS Snow-Cover Products. Remote Sens Environ 2002, 83, 181–194. [Google Scholar] [CrossRef]
  20. Lopez, P.; Sirguey, P.; Arnaud, Y.; Pouyaud, B.; Chevallier, P. Snow Cover Monitoring in the Northern Patagonia Icefield Using MODIS Satellite Images (2000–2006). Glob Planet Change 2008, 61, 103–116. [Google Scholar] [CrossRef]
  21. Sirguey, P.; Mathieu, R.; Arnaud, Y. Subpixel Monitoring of the Seasonal Snow Cover with MODIS at 250 m Spatial Resolution in the Southern Alps of New Zealand: Methodology and Accuracy Assessment. Remote Sens Environ 2009, 113, 160–181. [Google Scholar] [CrossRef]
  22. Musselman, K.N.; Molotch, N.P.; Brooks, P.D. Effects of Vegetation on Snow Accumulation and Ablation in a Mid-Latitude Sub-Alpine Forest. Hydrol Process 2008, 22, 2767–2776. [Google Scholar] [CrossRef]
  23. Varhola, A.; Coops, N.C.; Weiler, M.; Moore, R.D. Forest Canopy Effects on Snow Accumulation and Ablation: An Integrative Review of Empirical Results. J Hydrol (Amst) 2010, 392, 219–233. [Google Scholar] [CrossRef]
  24. Sirguey, P.; Mathieu, R.; Arnaud, Y. Subpixel Monitoring of the Seasonal Snow Cover with MODIS at 250 m Spatial Resolution in the Southern Alps of New Zealand: Methodology and Accuracy Assessment. Remote Sens Environ 2009, 113, 160–181. [Google Scholar] [CrossRef]
  25. Dumont, M.; Gardelle, J.; Sirguey, P.; Guillot, A.; Six, D.; Rabatel, A.; Arnaud, Y. Linking Glacier Annual Mass Balance and Glacier Albedo Retrieved from MODIS Data. Cryosphere 2012, 6, 1527–1539. [Google Scholar] [CrossRef]
  26. GCOS, G. Systematic Observation Requirements for Satellite-Based Data Products for Climate 2011 Update: Supplemental Details to the Satellite-Based Component of the “Implementation Plan for the Global Observing System for Climate in Support of the UNFCCC (2010 Upd. In Tech. rep; World Meteorological Organisation (WMO) 7 bis, 2011.
  27. Dozier, J.; Painter, T.H. MULTISPECTRAL AND HYPERSPECTRAL REMOTE SENSING OF ALPINE SNOW PROPERTIES. Annu Rev Earth Planet Sci 2004, 32, 465–494. [Google Scholar] [CrossRef]
  28. Painter, T.H.; Rittger, K.; McKenzie, C.; Slaughter, P.; Davis, R.E.; Dozier, J. Retrieval of Subpixel Snow Covered Area, Grain Size, and Albedo from MODIS. Remote Sens Environ 2009, 113, 868–879. [Google Scholar] [CrossRef]
  29. Pérez, T.; Mattar, C.; Fuster, R. Decrease in Snow Cover over the Aysén River Catchment in Patagonia, Chile. Water (Basel) 2018, 10, 619. [Google Scholar] [CrossRef]
  30. Stehr, A.; Aguayo, M. Snow Cover Dynamics in Andean Watersheds of Chile (32.0–39.5° S) during the Years 2000–2016. Hydrol Earth Syst Sci 2017, 21, 5111–5126. [Google Scholar] [CrossRef]
  31. Malmros, J.K.; Mernild, S.H.; Wilson, R.; Tagesson, T.; Fensholt, R. Snow Cover and Snow Albedo Changes in the Central Andes of Chile and Argentina from Daily MODIS Observations (2000–2016). Remote Sens Environ 2018, 209, 240–252. [Google Scholar] [CrossRef]
  32. Mernild, S.H.; Beckerman, A.P.; Yde, J.C.; Hanna, E.; Malmros, J.K.; Wilson, R.; Zemp, M. Mass Loss and Imbalance of Glaciers along the Andes Cordillera to the Sub-Antarctic Islands. Glob Planet Change 2015, 133, 109–119. [Google Scholar] [CrossRef]
  33. Garreaud, R.; Lopez, P.; Minvielle, M.; Rojas, M.; Garreaud, R.; Lopez, P.; Minvielle, M.; Rojas, M. Large-Scale Control on the Patagonian Climate. J Clim 2013, 26, 215–230. [Google Scholar] [CrossRef]
  34. Rozzi, R.; Massardo, F.; Gallardo, M.R.; Plana, J. La Reserva de Biosfera Cabo de Hornos: Un Desafío Para La Conservación de La Biodiversidad e Implementación Del Desarrollo Sustentable En El Extremo Austral de América. Environ Res 2007, 35, 55–70. [Google Scholar]
  35. Aguirre, F.; Squeo, F.A.; López, D.; Grego, R.D.; Buma, B.; Carvajal, D.; Jaña, R.; Casassa, G.; Rozzi, R. Gradientes Climáticos y Su Alta Influencia En Los Ecosistemas Terrestres de La Reserva de La Biosfera Cabo de Hornos, Chile. In Proceedings of the Anales del Instituto de la Patagonia; 2021; Vol. 49. [Google Scholar]
  36. Aravena, J.C.; Luckman, B.H. Spatio-Temporal Rainfall Patterns in Southern South America. International Journal of Climatology 2009, 29, 2106–2120. [Google Scholar] [CrossRef]
  37. Davies, B.J.; Darvill, C.M.; Lovell, H.; Bendle, J.M.; Dowdeswell, J.A.; Fabel, D.; García, J.L.; Geiger, A.; Glasser, N.F.; Gheorghiu, D.M.; et al. The Evolution of the Patagonian Ice Sheet from 35 Ka to the Present Day (PATICE). Earth Sci Rev 2020, 204, 103152. [Google Scholar] [CrossRef]
  38. Kilian, R.; Lamy, F. A Review of Glacial and Holocene Paleoclimate Records from Southernmost Patagonia (49–55°S). Quat Sci Rev 2012, 53, 1–23. [Google Scholar] [CrossRef]
  39. Markgraf, V.; Huber, U.M. Late and Postglacial Vegetation and Fire History in Southern Patagonia and Tierra Del Fuego. Palaeogeogr Palaeoclimatol Palaeoecol 2010, 297, 351–366. [Google Scholar] [CrossRef]
  40. Villalba, R.; Grosjean, M.; Kiefer, T. Long-Term Multi-Proxy Climate Reconstructions and Dynamics in South America (LOTRED-SA): State of the Art and Perspectives. Palaeogeogr Palaeoclimatol Palaeoecol 2009, 281, 175–179. [Google Scholar] [CrossRef]
  41. Mayewski, P.A.; Meredith, M.P.; Summerhayes, C.P.; Turner, J.; Worby, A.; Barrett, P.J.; Casassa, G.; Bertler, N.A.N.; Bracegirdle, T.; Naveira Garabato, A.C.; et al. State of the Antarctic and Southern Ocean Climate System. Reviews of Geophysics 2009, 47, RG1003. [Google Scholar] [CrossRef]
  42. Marshall, G.J. Trends in the Southern Annular Mode from Observations and Reanalyses. J Clim 2003, 16, 4134–4143. [Google Scholar] [CrossRef]
  43. Sauter, T. Revisiting Extreme Precipitation Amounts over Southern South America and Implications for the Patagonian Icefields. Hydrol Earth Syst Sci 2020, 24, 2003–2016. [Google Scholar] [CrossRef]
  44. Carrasco, J.F.; Casassa, G.; Rivera, A. Meteorological and Climatological Aspects of the Southern Patagonia Icefield. In; Springer, Boston, MA, 2002; pp. 29–41.
  45. Schneider, C.; Glaser, M.; Kilian, R.; Santana, A.; Butorovic, N.; Casassa, G. Weather Observations Across the Southern Andes at 53°S. Phys Geogr 2003, 24, 97–119. [Google Scholar] [CrossRef]
  46. Garreaud, R.; Lopez, P.; Minvielle, M.; Rojas, M. Large-Scale Control on the Patagonian Climate. J Clim 2013, 26, 215–230. [Google Scholar] [CrossRef]
  47. Weidemann, Stephanie. ; Sauter, T.; Kilian, R.; Steger, D.; Butorovic, N.; Schneider, C. A 17-Year Record of Meteorological Observations Across the Gran Campo Nevado Ice Cap in Southern Patagonia, Chile, Related to Synoptic Weather Types and Climate Modes. Front Earth Sci (Lausanne) 2018, 6, 53. [Google Scholar] [CrossRef]
  48. Alvarez-Garreton, C.; Mendoza, P.A.; Boisier, J.P.; Addor, N.; Galleguillos, M.; Zambrano-Bigiarini, M.; Lara, A.; Cortes, G.; Garreaud, R.; McPhee, J. The CAMELS-CL Dataset: Catchment Attributes and Meteorology for Large Sample Studies-Chile Dataset. Hydrol Earth Syst Sci 2018, 22, 5817–5846. [Google Scholar] [CrossRef]
  49. Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L.; et al. The Shuttle Radar Topography Mission. Reviews of Geophysics 2007, 45, RG2004. [Google Scholar] [CrossRef]
  50. Pohl, C.; Van Genderen, J. Remote Sensing Image Fusion: A Practical Guide; Crc Press, 2017. ISBN 1498730035.
  51. Canty, M.J. Image Analysis, Classification and Change Detection in Remote Sensing: With Algorithms for Python; Fourth.; CRC Press, 2019. ISBN 0429875347.
  52. Shettigara, V.K. A Generalized Component Substitution Technique for Spatial Enhancement of Multispectral Images Using a Higher Resolution Data Set. Photogramm Eng Remote Sensing 1992, 58, 561–567. [Google Scholar]
  53. Nunez, J.; Otazu, X.; Fors, O.; Prades, A.; Pala, V.; Arbiol, R. Multiresolution-Based Image Fusion with Additive Wavelet Decomposition. IEEE Transactions on Geoscience and Remote Sensing 1999, 37, 1204–1211. [Google Scholar] [CrossRef]
  54. Ranchin, T.; Aiazzi, B.; Alparone, L.; Baronti, S.; Wald, L. Image Fusion—the ARSIS Concept and Some Successful Implementation Schemes. ISPRS Journal of Photogrammetry and Remote Sensing 2003, 58, 4–18. [Google Scholar] [CrossRef]
  55. Aiazzi, B.; Alparone, L.; Baronti, S.; Garzelli, A. Context-Driven Fusion of High Spatial and Spectral Resolution Images Based on Oversampled Multiresolution Analysis. IEEE Transactions on Geoscience and Remote Sensing 2002, 40, 2300–2312. [Google Scholar] [CrossRef]
  56. Sales, M.H.; Souza, C.M.; Kyriakidis, P.C. Fusion of MODIS Images Using Kriging With External Drift. IEEE Transactions on Geoscience and Remote Sensing 2013, 51, 2250–2259. [Google Scholar] [CrossRef]
  57. Wang, Q.; Shi, W.; Atkinson, P.M.; Zhao, Y. Downscaling MODIS Images with Area-to-Point Regression Kriging. Remote Sens Environ 2015, 166, 191–204. [Google Scholar] [CrossRef]
  58. Gao, F.; Kustas, W.; Anderson, M. A Data Mining Approach for Sharpening Thermal Satellite Imagery over Land. Remote Sens (Basel) 2012, 4, 3287–3319. [Google Scholar] [CrossRef]
  59. Hu, M.; Huang, Y. Atakrig: An R Package for Multivariate Area-to-Area and Area-to-Point Kriging Predictions. Comput Geosci 2020, 139, 104471. [Google Scholar] [CrossRef]
  60. Pardo-Iguzquiza, E.; Atkinson, P.M.; Chica-Olmo, M. DSCOKRI: A Library of Computer Programs for Downscaling Cokriging in Support of Remote Sensing Applications. Comput Geosci 2010, 36, 881–894. [Google Scholar] [CrossRef]
  61. Jin, Y.; Ge, Y.; Wang, J.; Heuvelink, G.; Wang, L. Geographically Weighted Area-to-Point Regression Kriging for Spatial Downscaling in Remote Sensing. Remote Sens (Basel) 2018, 10, 579. [Google Scholar] [CrossRef]
  62. Zuur, A.; Ieno, E.N.; Walker, N.; Saveliev, A.A.; Smith, G.M. Dealing with Heterogeneity. In Mixed effects models and extensions in ecology with R; Springer Science & Business Media, 2009; pp. 71–101. ISBN 0387874585.
  63. Hijmans, R.J.; van Etten, J. Raster: Geographic Data Analysis and Modeling. R package version 2016, 2. [Google Scholar]
  64. Pinheiro, J.; Bates, D.; DebRoy, S.; Sarkar, D.; Heisterkamp, S.; Van Willigen, B.; Maintainer, R. Package ‘Nlme. ’ Linear and nonlinear mixed effects models, version 2017, 3. [Google Scholar]
  65. Pebesma, E.J. Multivariable Geostatistics in S: The Gstat Package. Comput Geosci 2004, 30, 683–691. [Google Scholar] [CrossRef]
  66. Gollini, I.; Lu, B.; Charlton, M.; Brunsdon, C.; Harris, P. GWmodel: An R Package for Exploring Spatial Heterogeneity Using Geographically Weighted Models. J Stat Softw 2013, 63, 1–50. [Google Scholar] [CrossRef]
  67. Wang, Z.; Bovik, A.C. A Universal Image Quality Index. IEEE Signal Process Lett 2002, 9, 81–84. [Google Scholar] [CrossRef]
  68. Alparone, L.; Aiazzi, B.; Baronti, S.; Garzelli, A.; Nencini, F.; Selva, M. Multispectral and Panchromatic Data Fusion Assessment without Reference. Photogramm Eng Remote Sensing 2008, 74, 193–200. [Google Scholar] [CrossRef]
  69. Somers, B.; Asner, G.P.; Tits, L.; Coppin, P. Endmember Variability in Spectral Mixture Analysis: A Review. Remote Sens Environ 2011, 115, 1603–1616. [Google Scholar] [CrossRef]
  70. Roberts, D.A.; Gardner, M.; Church, R.; Ustin, S.; Scheer, G.; Green, R.O. Mapping Chaparral in the Santa Monica Mountains Using Multiple Endmember Spectral Mixture Models. Remote Sens Environ 1998, 65, 267–279. [Google Scholar] [CrossRef]
  71. Leutner, B.; Horning, N.; Schwalb-Willmann, J.; Hijmans, R.J. RStoolbox: Tools for Remote Sensing Data Analysis. R package version 0.2.6 2019. [Google Scholar]
  72. Painter, T.H.; Roberts, D.A.; Green, R.O.; Dozier, J. The Effect of Grain Size on Spectral Mixture Analysis of Snow-Covered Area from AVIRIS Data. Remote Sens Environ 1998, 65, 320–332. [Google Scholar] [CrossRef]
  73. Dudley, K.L.; Roth, K.L.; Coates, A.R. A Multi-Temporal Spectral Library Approach for Mapping Vegetation Species across Spatial and Temporal Phenological Gradients. Remote Sens Environ 2015, 167, 121–134. [Google Scholar] [CrossRef]
  74. Congedo, L. Semi-Automatic Classification Plugin Documentation. Release 2016, 4, 29. [Google Scholar]
  75. QGIS Development-Teams QGIS Geographic Information System. Open Source Geospatial Foundation Project 2020.
  76. Dozier, J.; Painter, T.H.; Rittger, K.; Frew, J.E. Time-Space Continuity of Daily Maps of Fractional Snow Cover and Albedo from MODIS. Adv Water Resour 2008, 31, 1515–1526. [Google Scholar] [CrossRef]
  77. Redpath, T.A.N.; Sirguey, P.; Cullen, N.J. Characterising Spatio-Temporal Variability in Seasonal Snow Cover at a Regional Scale from MODIS Data: The Clutha Catchment, New Zealand. Hydrol Earth Syst Sci 2019, 23, 3189–3217. [Google Scholar] [CrossRef]
  78. Chylek, P.; McCabe, M.; Dubey, M.K.; Dozier, J. Remote Sensing of Greenland Ice Sheet Using Multispectral Near-Infrared and Visible Radiances. J Geophys Res 2007, 112, D24S20. [Google Scholar] [CrossRef]
  79. Bormann, K.J.; McCabe, M.F.; Evans, J.P. Satellite Based Observations for Seasonal Snow Cover Detection and Characterisation in Australia. Remote Sens Environ 2012, 123, 57–71. [Google Scholar] [CrossRef]
  80. Medeiros, E.S. de; de Lima, R.R.; Olinda, R.A. de; Dantas, L.G.; Santos, C.A.C. dos Space–Time Kriging of Precipitation: Modeling the Large-Scale Variation with Model GAMLSS. Water (Basel) 2019, 11, 2368. [Google Scholar] [CrossRef]
  81. Wang, Y. Smoothing Splines: Methods and Applications; CRC press, 2011. ISBN 1420077562.
  82. Pebesma, E. Spacetime: Spatio-Temporal Data in R. J Stat Softw 2012, 51, 1–30. [Google Scholar] [CrossRef]
  83. Easterling, D.R.; Kunkel, K.E.; Wehner, M.F.; Sun, L. Detection and Attribution of Climate Extremes in the Observed Record. Weather Clim Extrem 2016, 11, 17–27. [Google Scholar] [CrossRef]
  84. Mudelsee, M. Trend Analysis of Climate Time Series: A Review of Methods. Earth Sci Rev 2019, 190, 310–322. [Google Scholar] [CrossRef]
  85. Rodell, M.; Houser, P.R. Updating a Land Surface Model with MODIS-Derived Snow Cover. J Hydrometeorol 2004, 5, 1064–1075. [Google Scholar] [CrossRef]
  86. Wood, S.N. Generalized Additive Models: An Introduction with R 2 An Introduction with R. Generalized Additive Models 2017, 10, 9781315370279. [Google Scholar]
  87. Mann, B.H. Non-Parametric Test Against Trend. Econometrica 1945, 13, 245–259. [Google Scholar] [CrossRef]
  88. Sen, P.K. Estimates of the Regression Coefficient Based on Kendall’s Tau. J Am Stat Assoc 1968, 63, 1379–1389. [Google Scholar] [CrossRef]
  89. Helsel, D.R.; Hirsch, R.M. Echniques of Water Resources Investigations, Book 4, Chapter A3; U.S. Geological Survey., 2002. ISBN 0444814639/9780444814630.
  90. Rosenbluth, B.N.; Fuenzalida, H.A.; Aceituno, P. RECENT TEMPERATURE VARIATIONS IN SOUTHERN SOUTH AMERICA. INTERNATIONAL JOURNAL OF CLIMATOLOGY Int. J. Climatol 1997, 17, 67–85. [Google Scholar] [CrossRef]
  91. Hyndman, R.; Koehler, A.B.; Ord, J.K.; Snyder, R.D. Forecasting with Exponential Smoothing: The State Space Approach; Springer Science & Business Media, 2008. ISBN 3540719180.
  92. Chatfield, C. The Analysis of Time Series : An Introduction, Sixth Edition. The Analysis of Time Series 2003. [Google Scholar] [CrossRef]
  93. Jassby, A.D.; Cloern, J.E. Wq: Exploring Water Quality Monitoring Data 2017.
  94. Frey, R.A.; Ackerman, S.A.; Liu, Y.; Strabala, K.I.; Zhang, H.; Key, J.R.; Wang, X. Cloud Detection with MODIS. Part I: Improvements in the MODIS Cloud Mask for Collection 5. J Atmos Ocean Technol 2008, 25, 1057–1072. [Google Scholar] [CrossRef]
  95. Wachter, P.; Beck, C.; Philipp, A.; Höppner, K.; Jacobeit, J. Spatiotemporal Variability of the Southern Annular Mode and Its Influence on Antarctic Surface Temperatures. Journal of Geophysical Research: Atmospheres 2020, 125, e2020JD033818. [Google Scholar] [CrossRef]
  96. Clem, K.R.; Renwick, J.A.; McGregor, J.; Fogt, R.L. The Relative Influence of ENSO and SAM on Antarctic Peninsula Climate. Journal of Geophysical Research: Atmospheres 2016, 121, 9324–9341. [Google Scholar] [CrossRef]
  97. Garreaud, R. Record-Breaking Climate Anomalies Lead to Severe Drought and Environmental Disruption in Western Patagonia in 2016. Clim Res 2018, 74, 217–229. [Google Scholar] [CrossRef]
  98. Turner, J.; Lu, H.; White, I.; King, J.C.; Phillips, T.; Hosking, J.S.; Bracegirdle, T.J.; Marshall, G.J.; Mulvaney, R.; Deb, P. Absence of 21st Century Warming on Antarctic Peninsula Consistent with Natural Variability. Nature 2016, 535, 411–415. [Google Scholar] [CrossRef] [PubMed]
  99. Schneider, C.; Gies, D. Effects of El Niño–Southern Oscillation on Southernmost South America Precipitation at 53°S Revealed from NCEP–NCAR Reanalyses and Weather Station Data. International Journal of Climatology 2004, 24, 1057–1076. [Google Scholar] [CrossRef]
  100. Srur, A.M.; Villalba, R.; Rodríguez-Catón, M.; Amoroso, M.M.; Marcotti, E. Climate and Nothofagus Pumilio Establishment at Upper Treelines in the Patagonian Andes. Front Earth Sci (Lausanne) 2018, 6, 57. [Google Scholar] [CrossRef]
  101. Sato, K.; Inoue, J. Comparison of Arctic Sea Ice Thickness and Snow Depth Estimates from CFSR with in Situ Observations. Clim Dyn 2018, 50, 289–301. [Google Scholar] [CrossRef]
  102. Wang, C.; Graham, R.M.; Wang, K.; Gerland, S.; Granskog, M.A. Comparison of ERA5 and ERA-Interim near-Surface Air Temperature, Snowfall and Precipitation over Arctic Sea Ice: Effects on Sea Ice Thermodynamics and Evolution. Cryosphere 2019, 13, 1661–1679. [Google Scholar] [CrossRef]
  103. Muñoz-Sabater, J.; Dutra, E.; Agustí-Panareda, A.; Albergel, C.; Arduini, G.; Balsamo, G.; Boussetta, S.; Choulga, M.; Harrigan, S.; Hersbach, H.; et al. ERA5-Land: A State-of-the-Art Global Reanalysis Dataset for Land Applications. Earth Syst Sci Data 2021, 13, 4349–4383. [Google Scholar] [CrossRef]
  104. Bravo, C.; Ross, A.N.; Quincey, D.J.; Cisternas, S.; Rivera, A. Surface Ablation and Its Drivers along a West–East Transect of the Southern Patagonia Icefield. Journal of Glaciology 2021, 1–14. [Google Scholar] [CrossRef]
  105. Quilodrán, C.S.; Sandvig, E.M.; Aguirre, F.; de Aguilar, J.R.; Barroso, O.; Vásquez, R.A.; Rozzi, R. The Extreme Rainfall Gradient of the Cape Horn Biosphere Reserve and Its Impact on Forest Bird Richness. Biodivers Conserv 2022, 1–15. [Google Scholar] [CrossRef]
  106. Schaefer, M.; Machguth, H.; Falvey, M.; Casassa, G. Modeling Past and Future Surface Mass Balance of the Northern Patagonia Icefield. J Geophys Res Earth Surf 2013, 118, 571–588. [Google Scholar] [CrossRef]
  107. Bozkurt, D.; Rojas, M.; Boisier, J.P.; Rondanelli, R.; Garreaud, R.; Gallardo, L. Dynamical Downscaling over the Complex Terrain of Southwest South America: Present Climate Conditions and Added Value Analysis. Clim Dyn 2019, 53, 6745–6767. [Google Scholar] [CrossRef]
Preprints 82050 g001
Figure 1. Study area within Brunswick Peninsula, showing the 4 selected basins. White circles are runoff stations and yellow circles are referential sites for ski activities (Cerro Mirador and Tres Morros). Contour line interval (in orange) is 100 m. Figure designed using open source QGIS v_3.162.
Figure 1. Study area within Brunswick Peninsula, showing the 4 selected basins. White circles are runoff stations and yellow circles are referential sites for ski activities (Cerro Mirador and Tres Morros). Contour line interval (in orange) is 100 m. Figure designed using open source QGIS v_3.162.
Preprints 82050 g001
Preprints 82050 g002
Figure 2. Snow spectral signatures of different grain sizes and the MODIS bands. Image modified from Painter et al. (1998) [72].
Figure 2. Snow spectral signatures of different grain sizes and the MODIS bands. Image modified from Painter et al. (1998) [72].
Preprints 82050 g002
Preprints 82050 g003
Figure 3. Reflectance values of MODIS bands in snow and snowless imagery. The reflectance values ([0,1]) are stretched to match the MODIS reflectance scale ([-100, 16000]).
Figure 3. Reflectance values of MODIS bands in snow and snowless imagery. The reflectance values ([0,1]) are stretched to match the MODIS reflectance scale ([-100, 16000]).
Preprints 82050 g003
Preprints 82050 g004
Figure 4. Variograms of linear regression residuals for each coarse band at Brunswick Peninsula watersheds. The colours on the upper left image represent the reflectance values (orange high, black low values) of band 1 of an image from September 25, 2015.
Figure 4. Variograms of linear regression residuals for each coarse band at Brunswick Peninsula watersheds. The colours on the upper left image represent the reflectance values (orange high, black low values) of band 1 of an image from September 25, 2015.
Preprints 82050 g004
Preprints 82050 g005
Figure 5. Spectral fusion results using the GLS-ATAK models. From left to right: Residual of ATAK at 250 m, GLS at 250 m, GLA-ATAK at 250 m and coarse band at 500 m. The 16 bit reflectance units for each band are expressed at a scale of [-100, 16000], and the coordinates of the map are in UTM. MODIS image of September 25 of 2015.
Figure 5. Spectral fusion results using the GLS-ATAK models. From left to right: Residual of ATAK at 250 m, GLS at 250 m, GLA-ATAK at 250 m and coarse band at 500 m. The 16 bit reflectance units for each band are expressed at a scale of [-100, 16000], and the coordinates of the map are in UTM. MODIS image of September 25 of 2015.
Preprints 82050 g005
Preprints 82050 g006
Figure 6. Endmember signatures and their seasonal variation at Brunswick Peninsula. Reflectances are in % and represent spectral albedo for each MODIS band.
Figure 6. Endmember signatures and their seasonal variation at Brunswick Peninsula. Reflectances are in % and represent spectral albedo for each MODIS band.
Preprints 82050 g006
Preprints 82050 g007
Figure 7. Spectral mixture analysis (SMA): (A) elevation model; (B) fractional snow cover for each pixel; (C) broadband albedo; and (D) average grain size of snow for each pixel (image for (B), (C) and (D) dating from September 25, 2015).
Figure 7. Spectral mixture analysis (SMA): (A) elevation model; (B) fractional snow cover for each pixel; (C) broadband albedo; and (D) average grain size of snow for each pixel (image for (B), (C) and (D) dating from September 25, 2015).
Preprints 82050 g007
Preprints 82050 g008
Figure 8. Snow and cloud cover, and new snow classifications. The top row shows RGB true colour images; the second row, cloud mask product of MODIS (MOD35_L2); third row, first snow classification based on NDSI (snow threshold ≥ 0.4); fourth row, second snow classification based on MADI (snow threshold ≥ 6); fifth row, third snow classification based on fractional snow cover fSC (snow threshold ≥ 0.2); bottom row, the snow mask generated based on the condition that all three snow classification indexes must agree.
Figure 8. Snow and cloud cover, and new snow classifications. The top row shows RGB true colour images; the second row, cloud mask product of MODIS (MOD35_L2); third row, first snow classification based on NDSI (snow threshold ≥ 0.4); fourth row, second snow classification based on MADI (snow threshold ≥ 6); fifth row, third snow classification based on fractional snow cover fSC (snow threshold ≥ 0.2); bottom row, the snow mask generated based on the condition that all three snow classification indexes must agree.
Preprints 82050 g008
Preprints 82050 g009
Figure 9. Snow fraction (fSC) spatial-temporal reconstruction, before (top) and after the interpolation (bottom).
Figure 9. Snow fraction (fSC) spatial-temporal reconstruction, before (top) and after the interpolation (bottom).
Preprints 82050 g009
Preprints 82050 g010
Figure 10. Snow height measurements at AWS Tres Morros. Each black line represents the time when the snow height decreases below 20 cm.
Figure 10. Snow height measurements at AWS Tres Morros. Each black line represents the time when the snow height decreases below 20 cm.
Preprints 82050 g010
Preprints 82050 g011
Figure 11. Snow cover variability in Brunswick Peninsula from MODIS reconstruction for the entire study area. (A) Daily snow cover area by year as boxplots, with empty circles representing outliers. (B) Monthly snow cover area as both km2 and percentage. The horizontal pink line represents 10% of the total study area as a definition of the snow season and the numbers in blue above represent the number of months during the snow season. (C) The blue line represents the number of days per year that the snow cover exceeds 10% of the total study area. The orange dashed line represents the exponential filter of the number of days per year for the whole series, that shows a significant trend of +0.54 days/year. The red dashed line represents the exponential filter of the number of days per year for 2010-2020 period, showing a significant decreasing trend of -4.64 day/year. The green line presents the average of the snow season days for the period 2000-2020.
Figure 11. Snow cover variability in Brunswick Peninsula from MODIS reconstruction for the entire study area. (A) Daily snow cover area by year as boxplots, with empty circles representing outliers. (B) Monthly snow cover area as both km2 and percentage. The horizontal pink line represents 10% of the total study area as a definition of the snow season and the numbers in blue above represent the number of months during the snow season. (C) The blue line represents the number of days per year that the snow cover exceeds 10% of the total study area. The orange dashed line represents the exponential filter of the number of days per year for the whole series, that shows a significant trend of +0.54 days/year. The red dashed line represents the exponential filter of the number of days per year for 2010-2020 period, showing a significant decreasing trend of -4.64 day/year. The green line presents the average of the snow season days for the period 2000-2020.
Preprints 82050 g011
Preprints 82050 g012
Figure 12. Spatial and temporal variability of snow cover at Brunswick Peninsula from MODIS reconstruction. (A) Spatial variability of snow days per year as an average for the period 2000-2020; (B) yearly significant trend of reduction of snow days for each pixel (in days/year). (C,D) show the elevation component of the Generalized Additive Model (GAM). The red line is a spline fit, with confidence interval (at 95%) shown as blue dashed lines. (C) Model showing the altitudinal distribution of snow days as an average for the period 2000-2020; and (D) Model showing the altitudinal distribution of yearly trend of snow days.
Figure 12. Spatial and temporal variability of snow cover at Brunswick Peninsula from MODIS reconstruction. (A) Spatial variability of snow days per year as an average for the period 2000-2020; (B) yearly significant trend of reduction of snow days for each pixel (in days/year). (C,D) show the elevation component of the Generalized Additive Model (GAM). The red line is a spline fit, with confidence interval (at 95%) shown as blue dashed lines. (C) Model showing the altitudinal distribution of snow days as an average for the period 2000-2020; and (D) Model showing the altitudinal distribution of yearly trend of snow days.
Preprints 82050 g012
Preprints 82050 g013
Figure 13. Reanalysis and regional climate model validation. (A) scatter plot of daily mean temperature between ERA5 Land and AWS Tres Morros (2018-2020); (B) scatter plot of monthly mean temperature between ERA5 Land and AWS Tres Morros (2018-2020); and (C) scatter plot of monthly mean temperature between ERA5 Land and RegCM4 (2000-2005).
Figure 13. Reanalysis and regional climate model validation. (A) scatter plot of daily mean temperature between ERA5 Land and AWS Tres Morros (2018-2020); (B) scatter plot of monthly mean temperature between ERA5 Land and AWS Tres Morros (2018-2020); and (C) scatter plot of monthly mean temperature between ERA5 Land and RegCM4 (2000-2005).
Preprints 82050 g013
Table 1. Club Andino’s operational days 2010-2020.
Table 1. Club Andino’s operational days 2010-2020.
year 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020
season length 78 85 65 76 57 95 10 0 30 30 Closed due to pandemic
Table 4. The results of Ordinary Linear Squares (OLS) regression proposed by Wang et al. (2015) at Brunswick Peninsula.
Table 4. The results of Ordinary Linear Squares (OLS) regression proposed by Wang et al. (2015) at Brunswick Peninsula.
Snow Image Snowless Image
Bandi Bandi ~ Band1 Bandi ~ Band1
Band3 R2 = 0.96
RMSE = 0.229		 R2 = 0.91
RMSE = 0.066
Band4 R2 = 0.99
RMSE = 0.094		 R2 = 0.84
RMSE = 0.074
Band5 R2 = 0.39
RMSE = 0.264		 R2 = 0.22
RMSE = 0.165
Band6 R2 = 0.54
RMSE = 0.229		 R2 = 0.73
RMSE = 0.133
Band7 R2 = 0.41
RMSE = 0.409		 R2 = 0.93
RMSE = 0.075
Table 5. Downscaling using spectral fusion, showing the first term (Eq. 2) as a new regression relationship based on coarse MODIS bands 3 to 7 at 500 m, using as independent variables the coarse Band 1 and the resampled DEM at 500 m resolution ( D E M 500 ).
Table 5. Downscaling using spectral fusion, showing the first term (Eq. 2) as a new regression relationship based on coarse MODIS bands 3 to 7 at 500 m, using as independent variables the coarse Band 1 and the resampled DEM at 500 m resolution ( D E M 500 ).
OLS Regression
Band 3 B a n d 3 ~ B a n d 1 R2 = 0.96
RMSE = 0.230
Band 4 B a n d 4 ~ B a n d 1 R2 = 0.99
RMSE = 0.094
Band 5 B a n d 5 ~ B a n d 1 + D E M 500 R2 = 0.52
RMSE = 0.234
Band 6 B a n d 6 ~ B a n d 1 + D E M 500 R2 = 0.68
RMSE = 0.319
Band 7 B a n d 7 ~ B a n d 1 + D E M 500 R2 = 0.61
RMSE = 0.333
Table 6. Models and downscaling results for subset selection in the snow-covered image area at Brunswick Peninsula. All indexes are dimensionless: Spatial Distortion Index (Ds), Spectral Distortion Index (Dγ), Quality Index without Reference (QNR Index), and Universal Quality Index (UQI).
Table 6. Models and downscaling results for subset selection in the snow-covered image area at Brunswick Peninsula. All indexes are dimensionless: Spatial Distortion Index (Ds), Spectral Distortion Index (Dγ), Quality Index without Reference (QNR Index), and Universal Quality Index (UQI).
OLS500 GLS500 GWR500 OLS_ATAK250 GLS_ATAK250 GWR_ATAK250
Band_3500 R2 = 0.97
RMSE = 0.230
AIC = -175		 R2 = 0.97
RMSE = 0.372
AIC = -2877		 R2 = 0.99
RMSE = 0.110
AIC = -2500		 UQI= 0.92 UQI= 0.94 UQI= 0.93
Band_4500 R2 = 0.993
RMSE = 0.085
AIC = -3713		 R2 = 0.993
RMSE = 0.101
AIC = -6053		 R2 = 0.99
RMSE = 0.038
AIC = -6296		 UQI= 0.93 UQI= 0.93 UQI= 0.93
Band_5500 R2 = 0.64
RMSE = 0.152
AIC = -1654		 R2 = 0.63
RMSE = 0.153
AIC =-1727		 R2 = 0.97
RMSE = 0.047
AIC = -5404		 UQI= 0.89 UQI= 0.98 UQI= 0.89
Band_6500 R2 = 0.69
RMSE = 0.304
AIC = 822		 R2 = 0.65
RMSE = 0.392
AIC = 222		 R2 = 0.98
RMSE = 0.074
AIC = -3838		 UQI= 0.96 UQI= 0.99 UQI= 0.91
Band_7500 R2 = 0.58
RMSE = 0.302
AIC = 794		 R2 = 0.53
RMSE = 0.359
AIC = 498		 R2 = 0.97
RMSE = 0.080
AIC = -3545		 UQI= 0.97 UQI= 0.99 UQI= 0.89
Ds 0.042 0.043 0.039
Dγ 0.015 0.009 0.016
QNR Index 0.943 0.949 0.946
Table 7. Image selection for seasonal endmember extraction at Brunswick Peninsula.
Table 7. Image selection for seasonal endmember extraction at Brunswick Peninsula.
Season Image date
Summer January 17 2015
January 16 2016
Autumn April 29 2016
May 05 2016
Winter September 03 &11 2016
Spring October 16 & 18 2016
Table 8. AWS Tres Morros snow height measurements (Snow_h) vs fSC reconstruction from MODIS imagery (SF_Rec) for 2018-2020.
Table 8. AWS Tres Morros snow height measurements (Snow_h) vs fSC reconstruction from MODIS imagery (SF_Rec) for 2018-2020.
Relations evaluated
Linear
SF_Rec ~ aws_sh		 R2 = 0.30
RMSE = 0.113
AIC = -199.8
GAM
SF_Rec ~ s(aws_sh)		 R2 = 0.55
RMSE = 0.089
AIC = -252.4
GAM 20 cm
SF_Rec ~ s(aws_sh)		 R2 = 0.05
RMSE = 0.066
AIC = - 251.6
1
http://modis-snow-ice.gsfc.nasa.gov
2
QGIS Geographic Information System. QGIS Team (2017). QGIS Geographic Information System. Open Source Geospatial Foundation Project. Available Online at: https://qgis.org
3
http://www.subdere.gov.cl/división-administrativa-de-chile/gobierno-regional-de-magallanes-y-antártica-chilena
4
Available in https://wiki.earthdata.nasa.gov/display/DAS/HEG%3A++HDF-EOS+to+GeoTIFF+Conversion+Tool
5
https://search.earthdata.nasa.gov/search?q=SRTMGL1
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

© 2024 MDPI (Basel, Switzerland) unless otherwise stated