Snow Cover Evolution in the Gran Paradiso National Park, Italian Alps, Using the Earth Observation Data Cube

Mountainous regions are particularly vulnerable to climate change, and the impacts are already extensive and observable, the implications of which go far beyond mountain boundaries and the environmental sectors. Monitoring and understanding climate and environmental changes in mountain regions is, therefore, needed. One of the key variables to study is snow cover, since it represents an essential driver of many ecological, hydrological and socioeconomic processes in mountains. As remotely sensed data can contribute to filling the gap of sparse in-situ stations in high-altitude environments, a methodology for snow cover detection through time series analyses using Landsat satellite observations stored in an Open Data Cube is described in this paper, and applied to a case study on the Gran Paradiso National Park, in the western Italian Alps. In particular, this study presents a proof of concept of the preliminary version of the snow observation from space algorithm applied to Landsat data stored in the Swiss Data Cube. Implemented in an Earth Observation Data Cube environment, the algorithm can process a large amount of remote sensing data ready for analysis and can compile all Landsat series since 1984 into one single multi-sensor dataset. Temporal filtering methodology and multi-sensors analysis allows one to considerably reduce the uncertainty in the estimation of snow cover area using high-resolution sensors. The study highlights that, despite this methodology, the lack of available cloud-free images still represents a big issue for snow cover mapping from satellite data. Though accurate mapping of snow extent below cloud cover with optical sensors still represents a challenge, spatial and temporal filtering techniques and radar imagery for future time series analyses will likely allow one to reduce the current cloud cover issue.


Introduction
The latest scientific observations [1,2] highlight a clear warming of the global climate system in recent decades, directly affecting the atmosphere, land, oceans and the cryosphere. Mountainous environments are among the regions most sensitive to and most affected by climate change [3,4]. Several studies using both measured and modelled data show evidence that warming rates are amplified at higher elevations (e.g., [3] and references therein). Among the major effects of such warming are the shrinking of glaciers; reductions of snow cover extension, quantity and duration; and permafrost thawing [1, 2,5], bringing on consequences from water availability; e.g., [6], biodiversity and ecosystem functions and services.
To assess and understand environmental changes in mountains and to support decision-making and adaptation, regular and continuous monitoring is required. One of the Essential Climate Variables [7] deserving particular attention is snow cover. The reduction of snow cover causes an amplification of warming rates through snow-albedo feedback. Moreover, snow represents a crucial seasonal reserve of water for downstream areas and has a key role in many ecological processes.
Owing to the difficulty of monitoring high-altitude environments with in-situ station networks, satellite Earth observation (EO) data can be considered as an appropriate source to complement scattered in-situ measurements [8]. EO are well suited to map snow cover because of the good contrast of snow with most other natural surfaces, except some clouds [9,10]. Moreover, the global coverage and regular repeatability of measurements offered by satellite images allow experts to monitor the vast temporal and spatial variability of snow cover where ground measurements may be insufficient, expensive or even dangerous [10][11][12]. For more than 40 years, snow has been successfully mapped from space using a variety of sensors [9,10,[13][14][15][16][17]. Since the 1970s, the long-term data records make Landsat datasets frequently used to study snow cover around the world at medium resolution (15 to 100 m) [18][19][20][21]. With the advent of the European Space Agency's (ESA's) Sentinel constellation, it is now possible to have images every 5 days at 10-m resolution. This further enhances monitoring capabilities to provide nearly real-time information on several geophysical parameters [22].
To efficiently exploit the increasing availability of satellite EO data, Earth observations data cubes (EODC) have recently emerged [23,24]. They represent a solution to store, organize, manage and analyze large amounts of multi-sensor EO data. The ambition is to allow scientists, researchers and other possible users to harness big EO data, facilitating the access and use of analysis ready data (ARD) [25][26][27], which are consistently processed using the highest scientific standards for immediate analysis in applications and for time-series exploitation. The interest in that objective has been proven by various implementations that already exist, such as the Open Data Cube (ODC) [25], the EarthServer [28], the e-sensing platform [29], the JRC Earth Observation Data and Processing Platform (JEODPP) [30], the Copernicus Data and Information Access Services (DIAS) [31,32] and the Google Earth Engine (GEE) [33].
In this paper we show, for the first time, an application of the EO Swiss Data Cube (SDC), an EODC specifically developed for Switzerland, based on Landsat imagery. We apply a preliminary version of an algorithm for snow detection, hereafter referred to as snow observations from space (SOfS), to the SDC to extract information on snow cover area (SCA) in the period 1984-2018, focusing on the Gran Paradiso National Park (GPNP) in the western Italian Alps. Snow cover information for the Gran Paradiso National Park based on in-situ stations is very limited and not easily available. In this case remote sensing data can provide valuable additional information for park management.
The first objective of this study is to assess the characteristics of this new snow dataset, such as the frequency of available images and of the cloud-free observations per month. Owing to significant cloud obstruction in many satellite images, we then test an approach, part of the SOfS algorithm, to combine images into monthly aggregations, providing maximum snow cover area products for a given month. A similar approach to mitigate the impact of clouds has already been applied to MODIS data by combining the information provided by MODIS Terra and Aqua sensors at the weekly time scale (i.e., [34][35][36]). The SOfS algorithm also allows one to evaluate, on a seasonal (winter) scale, the number of cloud-free observations and the number of times snow was observed in the corresponding cloud-free scenes. Combining the information of the number of snow observations relative to the number of cloud-free observations, we produce "snow cover summaries" for each winter season from 1984 to 2018 and for the aggregation of all 34 winters seasons. Third, a tentative analysis of the temporal evolution of the percent area of the GPNP with high/low probability to observe snow in winter is also presented, though this has to be regarded with caution owing to the non-homogeneous frequency of the underlying observations. We finally identify and discuss the main benefits and limitations of the snow data cube and of the SOfS algorithm employed for snow discrimination. Possible strategies to further enhance the capabilities of the SOfS algorithm exploiting the EODC infrastructure are presented and discussed.
The paper is structured as follows: Section 2 describes the area of study and its climatic and environmental characteristics; Section 3 presents the Swiss Data Cube; Section 4 describes the methodology used to generate the snow cover dataset, including the procedure to identify clouds and the preliminary version of the SOfS algorithm, and describes how the snow cover monthly and seasonal summaries are generated; Section 5 shows the application of that methodology, discussing the characteristics of the snow cover data cube obtained for the GPNP area and presents the main results of the paper; Section 6 discusses the limitations and perspectives of this study and Section 7 concludes the paper.

Study Area and Its Climatic Characteristics
This study focuses on the Gran Paradiso National Park (GPNP) protected area, located in the western Italian Alps, encompassing the Aosta Valley and Piedmont regions, as shown in Figure 1. Established in 1922, the Italian national park covers a surface of about 720 km 2 [37], with elevations ranging from about 700 m a.s.l. to about 4000 m a.s.l. More than 75% of the park area lies above 2000 m a.s.l., and it is covered by snow for several months per year. A total of 60% of the park is covered by areas with scarce or no vegetation (rocks, screes and glaciers), while 20.2% is characterized by alpine vegetation (woods and shrubs) [38]. The pastures and meadows represent 17% of the Italian park and only 0.8% of the park is dedicated to urban areas and cultivated lands. The GPNP is home to almost 15 altitude lakes and to great richness and diversity of fauna and flora [38]. high/low probability to observe snow in winter is also presented, though this has to be regarded with caution owing to the non-homogeneous frequency of the underlying observations. We finally identify and discuss the main benefits and limitations of the snow data cube and of the SOfS algorithm employed for snow discrimination. Possible strategies to further enhance the capabilities of the SOfS algorithm exploiting the EODC infrastructure are presented and discussed. The paper is structured as follows: section 2 describes the area of study and its climatic and environmental characteristics; section 3 presents the Swiss Data Cube; section 4 describes the methodology used to generate the snow cover dataset, including the procedure to identify clouds and the preliminary version of the SOfS algorithm, and describes how the snow cover monthly and seasonal summaries are generated; section 5 shows the application of that methodology, discussing the characteristics of the snow cover data cube obtained for the GPNP area and presents the main results of the paper; section 6 discusses the limitations and perspectives of this study and section 7 concludes the paper.

Study Area and Its Climatic Characteristics
This study focuses on the Gran Paradiso National Park (GPNP) protected area, located in the western Italian Alps, encompassing the Aosta Valley and Piedmont regions, as shown in Figure 1. Established in 1922, the Italian national park covers a surface of about 720 km 2 [37], with elevations ranging from about 700 m a.s.l. to about 4000 m a.s.l. More than 75% of the park area lies above 2000 m a.s.l., and it is covered by snow for several months per year. A total of 60% of the park is covered by areas with scarce or no vegetation (rocks, screes and glaciers), while 20.2% is characterized by alpine vegetation (woods and shrubs) [38]. The pastures and meadows represent 17% of the Italian park and only 0.8% of the park is dedicated to urban areas and cultivated lands. The GPNP is home to almost 15 altitude lakes and to great richness and diversity of fauna and flora [38]. Though differences in elevation, slope and aspect between valleys can give rise to diverse microclimate conditions inside the park, the average climate of the area is characterized by a relatively low mean temperature, scarce precipitation and well-defined seasons. Figure 2 shows the Though differences in elevation, slope and aspect between valleys can give rise to diverse microclimate conditions inside the park, the average climate of the area is characterized by a relatively low mean temperature, scarce precipitation and well-defined seasons. Figure 2 shows the climatological annual cycle of temperature and precipitation in the GPNP area obtained from an observation-based gridded dataset specifically developed for the Greater Alpine Region, HISTALP [39,40]. This dataset provides monthly temperature and (total and solid) precipitation at a spatial resolution of 0.08 • latitude-longitude, corresponding to about 10 km, from 1780 to 2014 (for temperature) and from 1801 to 2014 (for precipitation). Based on this dataset and considering averages over the period 1950-2014 (Figure 2), the coldest months in the GPNP are December to March, with an average temperature of about −6 • C. The warmest months are July and August, with an average temperature of about 9 • C. Precipitation has maxima in spring/early summer (April, May and June, 110 mm/month) and autumn (particularly in November with about 180 mm/month), while the driest conditions along the year are found in July-August (71 mm/month) and in December to March (76 mm/month). climatological annual cycle of temperature and precipitation in the GPNP area obtained from an observation-based gridded dataset specifically developed for the Greater Alpine Region, HISTALP [39,40]. This dataset provides monthly temperature and (total and solid) precipitation at a spatial resolution of 0.08° latitude-longitude, corresponding to about 10 km, from 1780 to 2014 (for temperature) and from 1801 to 2014 (for precipitation). Based on this dataset and considering averages over the period 1950-2014 (Figure 2), the coldest months in the GPNP are December to March, with an average temperature of about −6 °C. The warmest months are July and August, with an average temperature of about 9 °C. Precipitation has maxima in spring/early summer (April, May and June, 110 mm/month) and autumn (particularly in November with about 180 mm/month), while the driest conditions along the year are found in July-August (71 mm/month) and in December to March (76 mm/month).

Figure 2.
Climatological annual cycle of temperature (left y-axis, black) and total precipitation (right y-axis, blue) spatially averaged over the Gran Paradiso National Park area, obtained from the HISTALP observation-based dataset. The climatology was obtained after temporally averaging the data over the period 1950-2014.
The temporal extent of the HISTALP dataset allows one to study changes in temperature and precipitation that occurred over the last few decades; for example, considering the averages in the period 1951-1980 ( Figure 3, left column plots), the averages in the period 1985-2014 (middle column) and the difference between the latter and the former (right column) of temperature (panels c, d, e), total precipitation (panels f, g, h) and solid precipitation (panels i, l, m). We considered only the winter season, including the months from December to February (DJF), which are particularly  The temporal extent of the HISTALP dataset allows one to study changes in temperature and precipitation that occurred over the last few decades; for example, considering the averages in the period 1951-1980 ( Figure 3, left column plots), the averages in the period 1985-2014 (middle column) and the difference between the latter and the former (right column) of temperature (panels c, d, e), total precipitation (panels f, g, h) and solid precipitation (panels i, l, m). We considered only the winter season, including the months from December to February (DJF), which are particularly relevant for the subsequent snow analysis. The Gran Paradiso National Park area experienced positive winter temperature changes in 1985-2014 compared to 1951-1980, ranging from about 1 • C to about 1.7 • C ( Figure 3e). The spatial distribution of the observed warming is positively correlated with the elevation distribution ( Figure 3b) with a Pearson's correlation coefficient greater than 0.9. Precipitation changes, for both total ( Figure 3h) and solid ( Figure 3m) precipitation are negative everywhere in the considered domain, indicating a decrease of precipitation in 1985-2014 compared to 1951-1980. The spatial pattern of change is similar for the total and solid precipitation, showing the largest negative values in the western and north-western part of the box encompassing the GPNP; overall, snowfall decreased slightly more than total precipitation.
In addition, Table 1 summarizes the temporal trends of winter (DJF) temperature and precipitation from 1950-2014 and for the more recent sub period (1985 to 2014), for which satellite data are available. For completeness, results for the months of February and April are also reported in the table, since these are key months for understanding snow dynamics in the Alps [41]. In particular, February is representative of a cold winter snowpack while April roughly corresponds to a time of maximum snowpack in the Alps [42].  In addition, Table 1 summarizes the temporal trends of winter (DJF) temperature and precipitation from 1950-2014 and for the more recent sub period (1985 to 2014), for which satellite data are available. For completeness, results for the months of February and April are also reported in the table, since these are key months for understanding snow dynamics in the Alps [41]. In particular, February is representative of a cold winter snowpack while April roughly corresponds to a time of maximum snowpack in the Alps [42].

The Swiss Data Cube
The Swiss Data Cube (SDC) is the second operational national Earth Observation Data Cube (EODC) worldwide, after Australia's [25,43,44]. It is an initiative supported by the Swiss Federal Office for the Environment and jointly developed, implemented and operated by the United Nations Environment Programme (UNEP)/GRID-Geneva and the University of Geneva. The SDC aims at supporting the Swiss government and institutions for their environmental monitoring and reporting mandates, and supports national scientific institutions to benefit from EO data for research and education. The SDC is built upon the ODC, a storage and analytical open source framework supported by Geoscience Australia, the Commonwealth Scientific and Industrial Research Organization, the United States Geological Survey (USGS), the National Aeronautics and Space Administration (NASA) and the Committee on Earth Observations Satellites (CEOS) [26,45,46]. The SDC currently holds 35 years of Landsat 5,7 and 8 data from the NASA Landsat program [47], covering the period from 1984 to 2019, updated daily, along with the entire Sentinel-2 archives [48,49] and a part of the Sentinel-1 archives [50,51].
The data available in the SDC are in an analysis ready form (after LEDAPS/LaSRC processing by USGS), meaning that all radiometric, geometric, solar and atmospheric corrections have already been applied and the data are spatially segmented into 30 × 30 m resolution grids that cover the area considered. All Landsat data used have undergone the same transformations to become analysis ready data (ARD) in order to get a consistent and harmonized multi-sensor dataset. ARD ensures that EO measurements are radiometrically comparable and geometrically aligned. Indeed, as mentioned by the CEOS [52], ARD is expected to limit as far as possible, barriers to interoperability both through time and with data from different sensors. Therefore, this multi-sensor dataset can be considered to be derived from a single sensor [53]. The unique ARD archive of Switzerland accounts for a total data volume of approximately 6 TB and more than 150 billion observations over the entire country that can be analyzed both in the spatial and temporal dimensions.
The Swiss Data Cube was originally developed employing eight Landsat tiles following the tiling notation system for Landsat data (the Worldwide Reference System-2 paths and rows) to cover the whole of Switzerland, also including parts of northern Italy, such as the Gran Paradiso National Park. Sentinel data, incorporated subsequently into the data cube, also cover the same domain. Despite the availability of Sentinel data in the archive, this proof of concept study focuses on Landsat sensors only. Indeed, for reasons of consistency and homogeneity, it is preferable for an initial test study to use data from the same satellite series with the same spatial resolution and having undergone the same ARD process. This study used the Collection 1 Level 1 data from Landsat 5 Thematic Mapper (TM) available from March 1984 to November 2011, Landsat 7 Enhanced Thematic Mapper Plus (ETM+) available since June 1999 and Landsat 8 Operational Land Imager (OLI)/Thermal Infrared Sensor (TIRS), available since March 2013. These datasets are characterized by a spatial resolution of 30 m and a revisit time of 16 days. Spectral bands of Landsat series satellites useful for snow detection and used in this study are the green and the shortwave Infrared (SWIR1) bands processed to orthorectified surface reflectance. Moreover, we employed the Collection 1 Level-1 Quality Assessment (QA) 16-bit Band, which provides information on the surface characteristics, cloud cover and sensor conditions ( Table 2) that can be used to generate masks according to the application and user needs. In our case, we used information from the QA product to create a cloud mask. The main characteristics of the dataset used in this study are summarized in Table 3.
For the generation of snow cover products over the GPNP, we processed more than 480 Landsat images for the months of December to April in the period 1984-2018.

Snow Observation from Space (SOfS) Algorithm
This study presents a preliminary version of the snow observation from space (SOfS) algorithm, which exploits consolidated techniques to detect snow cover from satellite images in the SDC and provides further tools to combine snow cover information over different temporal scales, for example, monthly and seasonal. A detailed description of the SOfS procedure is provided in the following sub-sections.

Cloud and Water Masks
Satellite observations are affected by many factors which can lead to poor observational quality, such as instrument failure or cloud interference [55]. During the winter season and in mountain regions especially, the presence of clouds reduces the capability of snow cover mapping through optical remote sensing [56,57]. Additionally, since clouds can have spectral signatures similar to snow, it is sometimes difficult to discriminate them from spaceborne multispectral sensors [16,21,58,59]. In order to exclude non-clear sky pixels, a cloud mask can be applied before evaluating the presence of snow. In this paper, to discriminate cloud covered pixels, we considered the cloud cover information included in the Landsat Collection 1 Level-1 Quality Assessment Band, derived using the C Function of Mask (CFMask) algorithm [60]. This algorithm provides information on the presence of clouds, along with a measure of the level of confidence that a pixel is actually cloud covered. The CFMask is very conservative and may have issues over bright targets, such as snow and ice surfaces, leading to possible overestimation of cloud cover, thus a loss of useful information for snow cover mapping. Therefore, we reprocessed the pixels and we considered as "cloud-covered," those previously identified as cloudy with a medium (34%-66%) or high (67%-100%) level of confidence, and pixels with cirrus clouds (at any level of confidence). More precisely, we considered as cloud covered, all pixels having the following attributes (see Table 2): "fill" (no data value), "low confidence cirrus," "high confidence cirrus," "medium confidence cloud" and "high confidence cloud".
Next, water pixels were filtered using a water mask from the Global Surface Water Explorer [61].

Snow Cover Detection and the Generation of Monthly Products
After filtering the dataset with a cloud and water mask, in order to detect the presence of snow in each pixel, we employed the well-known normalized difference snow index (NDSI) test [62][63][64] which has proven to be effective for the binary (i.e., for distinguishing between "snow" or "not snow") monitoring of snow cover [9,15,16,58,[65][66][67]. NDSI is defined as the reflectance difference between visible (green) and shortwave infrared (SWIR) wavelengths: NDSI = (r green − r SWIR )/(r green + r SWIR ) (1) where r green and r SWIR are the surface reflectances in the green and in the SWIR bands, respectively (see Table 3). The index ranges from −1 to +1. A pure snow pixel is characterized by a high NDSI, while a pixel with mixed elements (e.g., snow, water, vegetation, bare ground, etc.) is characterized by lower NDSI values [9]. A threshold test is then required to classify a pixel as snow covered or snow free. Based on a previous study performed using Landsat TM data in Sierra Nevada [64], pixels with at least 50% snow cover are found to have a NDSI value greater than 0.4. Though this value is commonly used [12], it may vary according to landscape conditions [15,20,[68][69][70]. For the Gran Paradiso National Park, we decided to apply the same threshold as the one adopted for a study over Switzerland (in preparation), i.e., 0.45, given the similarity of land cover and altitude characteristics between the two areas.
The low temporal resolution of Landsat satellites combined with the presence of clouds in mountainous environments makes it difficult to get cloud-free observations. Therefore, numerous techniques, such as multi-date composite and a combination of different sensors, have been developed to aggregate multiple observations over time [15,21,71,72]. In this paper, we combine the single NDSI maps into monthly aggregations (for each month from December to April) using a mosaicking method. As shown in Figure 4, this method takes the maximum NDSI value per pixel from any available image independently from the sensor during the considered month. The NDSI threshold test is finally applied to the maximum NDSI value. For each month (December to April) we then obtained one snow cover area (SCA) map, showing the maximum snow cover extent for that month. Figure 5 shows an example of the SCA map for February 2000, for which cloud-free and water-free pixels were classified as snow-covered (NDSI > 0.45) or snow-free (NDSI ≤ 0.45).
This technique minimizes, though without completely removing, the number of cloud-covered pixels in the monthly composite. Indeed, the percentage of scenes with less than 30% cloud cover increases from 51% to 63% after considering the maximum NDSI value by month. This multi-temporal and multi-sensor approach allows one to reduce limitations due to Landsat temporal resolution and to increase the number of pixels tested for the presence of snow [67,73].
A flow diagram of the SOfS algorithm which we set up is shown in Figure 6. we then obtained one snow cover area (SCA) map, showing the maximum snow cover extent for that month. Figure 5 shows an example of the SCA map for February 2000, for which cloud-free and water-free pixels were classified as snow-covered (NDSI > 0.45) or snow-free (NDSI ≤ 0.45). This technique minimizes, though without completely removing, the number of cloud-covered pixels in the monthly composite. Indeed, the percentage of scenes with less than 30% cloud cover increases from 51% to 63% after considering the maximum NDSI value by month. This multi-temporal and multi-sensor approach allows one to reduce limitations due to Landsat temporal resolution and to increase the number of pixels tested for the presence of snow [67,73]. difference snow index (NDSI) for a given month. From the available images (three in this example), NDSI was derived and the highest NDSI value was retained for the NDSI threshold test. This technique minimizes, though without completely removing, the number of cloud-covered pixels in the monthly composite. Indeed, the percentage of scenes with less than 30% cloud cover increases from 51% to 63% after considering the maximum NDSI value by month. This multi-temporal and multi-sensor approach allows one to reduce limitations due to Landsat temporal resolution and to increase the number of pixels tested for the presence of snow [67,73]. A flow diagram of the SOfS algorithm which we set up is shown in Figure 6.

Seasonal Snow Cover Summaries
Besides monthly aggregations, we also performed a seasonal-scale analysis, focusing on the winter season (DJF). For each winter season at the pixel level, we calculated the seasonal snow cover summary, Si, representing the relative frequency of snow cover observations in the corresponding monthly products, with respect to the total number of available cloud-free observations:

Seasonal Snow Cover Summaries
Besides monthly aggregations, we also performed a seasonal-scale analysis, focusing on the winter season (DJF). For each winter season at the pixel level, we calculated the seasonal snow cover summary, Si, representing the relative frequency of snow cover observations in the corresponding monthly products, with respect to the total number of available cloud-free observations: where nc (from 0 to 3) is the number of cloud-free observations in the three-monthly snow cover products and ns (0 to 3) is the number of snow occurrences in the three snow cover products. Figure 7 shows two examples of seasonal snow cover summaries for winters 1988/89 (upper panels) and 2001/02 (lower panels), the former characterized by substantial cloud cover in each of the monthly snow cover products, and the latter characterized by mostly cloud-free conditions. The number of cloud-free observations (nc) are displayed in Figure 7a Figure 7b,e: that corresponds to the number of times one pixel was classified as snow during the three winter months, with possible values of 0 (indicating absence of snow or no data either because of clouds or sensor technical issues), 1 (snow in only one out of three months), 2 (snow in two out of three months) or 3 (snow in all three months). Figure 7c,f shows the snow summaries (Si), expressing the percentage of snow cover observations with respect to the total number of cloud-free observations. In that case, each pixel can take one of the following values: 0%, (no snow cover in none of the three-monthly products), 33%, 50%, 66% (both snow cover and snow-free conditions in the three-monthly products) or 100% (snow cover in each of the cloud-free monthly products). The white pixels in Figure 7c,f are those with no data, or water or cloud covered in all three-monthly products.
The winter snow cover summary maps were employed to evaluate each season individually and explore the temporal evolution of the winter snow cover in the period 1984-2018. In addition to those time series, a snow cover summary considering all winter months in the period of study  was calculated. In the following, snow cover summary maps were interpreted in terms of the probability of observing snow cover in a given period of time in the GPNP.
To summarize, the preliminary version of the snow observation from space algorithm based on the NDSI approach, widely used to map snow cover in mountain regions, was applied to derive (i) monthly (December, January, February and April) SCA maps, defined as the maximum snow cover extent in the study area for that month; (ii) the winter (DJF) SCA analysis representing the percentage of snow observations with respect to the total number of cloud-free observations in all December, January and February in the period 1984-2018; (iii) winter (DJF) SCA time series, defined as the percentages of snow observations with respect to the total cloud-free observations during each winter season in the period 1984-2018. Implemented in an Earth Observation Data Cube environment, the algorithm could process a large amount of remote sensing data ready for analysis and benefited from all Landsat series since 1984 being put into one single multi-sensor dataset. Inspired by a methodology often used for MODIS sensors data, a multi-date composite of one month was applied to the snow dataset to reduce the impact of cloud cover. This proof-of-concept study tested, for the first time, the potential of the data cube infrastructure to monitor snow cover evolution in a small region of the Alps. In this context, Landsat images, with their long historical records and finer spatial resolution, become very attractive for multi-decade time series analyses, despite a 16-days repeat cycle. Applying the data cube methodology to EO datasets attempts to address new computation capabilities and to shift paradigm from scene-based analysis to pixel-based time series analysis. Time series analysis, with the aim to assess and monitor the variability and trends in snow cover extent in the GPNP.
observations (ns) is displayed in Figure 7b,e: that corresponds to the number of times one pixel was classified as snow during the three winter months, with possible values of 0 (indicating absence of snow or no data either because of clouds or sensor technical issues), 1 (snow in only one out of three months), 2 (snow in two out of three months) or 3 (snow in all three months). Figure 7c,f shows the snow summaries (Si), expressing the percentage of snow cover observations with respect to the total number of cloud-free observations. In that case, each pixel can take one of the following values: 0%, (no snow cover in none of the three-monthly products), 33%, 50%, 66% (both snow cover and snow-free conditions in the three-monthly products) or 100% (snow cover in each of the cloud-free monthly products). The white pixels in Figure 7c,f are those with no data, or water or cloud covered in all three-monthly products. The winter snow cover summary maps were employed to evaluate each season individually and explore the temporal evolution of the winter snow cover in the period 1984-2018. In addition to those time series, a snow cover summary considering all winter months in the period of study  was calculated. In the following, snow cover summary maps were interpreted in terms of the probability of observing snow cover in a given period of time in the GPNP.
To summarize, the preliminary version of the snow observation from space algorithm based on the NDSI approach, widely used to map snow cover in mountain regions, was applied to derive (i) monthly (December, January, February and April) SCA maps, defined as the maximum snow cover

Snow Cover and Climatic Data Correlation
As temperature and precipitation are the two principal drivers of snow [74], we found it interesting to evaluate their trends over the past three decades in the Gran Paradiso National Park and analyze them jointly with the snow cover information resulting from the application of the SOfS algorithm. We used the HISTALP climate data already described in Section 2. We calculated the regression and correlation coefficients between climatic and snow cover data by least-squares linear fitting over the time period 1984-2014, for the winter season and for the months of February and April separately. Correlations were calculated excluding months with a cloud cover exceeding 30% (15%) of the entire park's surface in the monthly (seasonal) analyses. Figure 8 shows in green colors, for every December, January, February, March and April in the period 1984-2018, the availability of Landsat scenes with less than 30% of cloud cover, while a cross indicates months with cloud cover exceeding 30% of the total GPNP area. The month of December is particularly under-represented with only 41% of months containing one or more scenes with less than 30% of cloud cover. Spring months (March and April) were also frequently affected by clouds. More than half of the cloudy scenes refer to the years before 2000 when only Landsat-5 (dark green in Figure 8) data were available. Since June 1999, the availability of two sensors (Landsat-5 and Landsat-7 for the period 2000-2012, and Landsat-7 and Landsat-8 for the period 2013-2018) considerably increased the number of scenes with less than 30% of cloud cover. The multi-temporal and multi-sensor approach employed in this study, i.e., the mosaicking method, allowed us to increase the number of valid monthly products with cloud cover of <30%. The last two columns of Figure 8 shows the percentage of monthly SCA products with less than 30% of cloud cover before and after the application of the mosaicking method. The percentage values referring to the monthly aggregations are always greater than the corresponding ones referring to the single satellite scenes. As an example, for April, the percentage of valid monthly SCA products (cloud cover <30% of the total area) increased from 47% to 56%, gaining three additional months with less than 30% of cloud cover out of the 34 years studied. For completeness, Figure 9 shows the percentage of cloud cover which still remains in the monthly SCA products after application of the mosaicking procedure. Despite the improvements brought by the applied methodology, on average, 58% of all considered months were valid using a cloud cover threshold of 30% (see Figure 9). The percentage of cloud cover varied widely from year to year and from month to month. However, since 2000, and particularly for the months of January, February and April, we observed a decrease in the percentage of cloud cover, owing to the availability of two different Landsat sensors. For example, for the month of April (February) the percentage of cloudy months (cloud cover >30% of the total area) decreased from 62% (50%) for the period 1985-2000 to 28% (22%) for the period 2001-2018.

Figure 8. Green cells show the availability of Landsat (TM, ETM+ and OLI) scenes for the Gran Paradiso
National Park with less than 30% cloud cover (after application of the CFMask) for December, January, February, March and April since 1984. Cells filled with more than one color indicate that valid scenes are available for more than one sensor. The last two columns represent the percentages of months over the full period (1984-2018) with at least one valid scene per month before and after the application of the mosaicking procedure respectively.
For completeness, Figure 9 shows the percentage of cloud cover which still remains in the monthly SCA products after application of the mosaicking procedure. Despite the improvements brought by the applied methodology, on average, 58% of all considered months were valid using a cloud cover threshold of 30% (see Figure 9). The percentage of cloud cover varied widely from year to year and from month to month. However, since 2000, and particularly for the months of January, February and April, we observed a decrease in the percentage of cloud cover, owing to the availability of two different Landsat sensors. For example, for the month of April (February) the percentage of cloudy months (cloud cover >30% of the total area) decreased from 62% (50%) for the period 1985-2000 to 28% (22%) for the period 2001-2018.

Winter Snow Cover Probability Map
Individual monthly maximum SCA products for December, January and February were employed to evaluate the "snow cover summary" over 34 years in the period 1984-2018 in the Gran Paradiso National Park area. Figure 10a shows the percentage of cloud-free observations per pixel in all December, January and February products in the period 1984-2018 (102 images considered in total; i.e., three winter months times 34 years). The percentage of cloud-free observations with respect to the total number of images available is at best equal to 88%, and rarely above 80%, as shown in Figure 10a in dark green colors. Most of the area (72%) had between 60% and 80% cloud-free observations, 27% of the area had between 40% and 60% cloud-free observations, while a small part of the park (0.63%) had less than 40% of cloud-free observations. This last class includes high-altitude areas (>3000 m a.s.l.) in the center of the park, that are particularly prone to cloud cover and showed few cloud-free observations available (20% to 40%) during the considered period. For those areas, highlighted in light green colors in Figure 10a, the uncertainty of snow observations was high and did not allow for a robust analysis. The percentage of occurrence in each of the five considered classes (0%-20%, 20%-40%, 40%-60%, 60%-80%, 80%-100%) is summarized in the pie chart inset into Figure 10a.
The percentage of snow observations in all December, January and February products in the period 1984-2018, with respect to the total number of cloud-free observations, is shown in Figure 10b. More than 80% of the park surface (about 512 km 2 SCA) is usually covered by snow during winter (the probability to observe snow is greater than 80%), while only 1.8% of the park surface (about 11 km 2 ) is mainly snow-free in winter (probability to observe snow < 20%). Active areas, where ephemeral snow covers the ground (probability to observe snow between 20% and 80%), can be considered transitional zones between snow-covered and snow-free areas, covering about 17% of the Gran Paradiso park (corresponding to 109 km 2 ). The probability to observe snow in winter months was evaluated as well, for different 500 m-thick elevation bins between 500 m and 4500 m, as shown in Figure 10c. For each elevation bin, the frequency of the five different probability classes (0%-20%, 20%-40%, 40%-60%, 60%-80% and 80%-100%) were expressed in terms of percent area in that class with respect to the total area of that elevation bin. Areas above 2500 meters, corresponding to about 58% of the total area of the GPNP, are characterized by a high probability (80%-100%) of observing snow cover during winter. Between 1500 m and 2000 m about 50% of the area still has a high (80%-100%) probability to be snow-covered, while below 1500 m the area with high snow probability rapidly decreases. At elevations higher than 1500 m, the area where snow cover is rarely observed is negligible. Figure 8. Green cells show the availability of Landsat (TM, ETM+ and OLI) scenes for the Gran Paradiso National Park with less than 30% cloud cover (after application of the CFMask) for December, January, February, March and April since 1984. Cells filled with more than one color indicate that valid scenes are available for more than one sensor. The last two columns represent the percentages of months over the full period  with at least one valid scene per month before and after the application of the mosaicking procedure respectively.
For completeness, Figure 9 shows the percentage of cloud cover which still remains in the monthly SCA products after application of the mosaicking procedure. Despite the improvements brought by the applied methodology, on average, 58% of all considered months were valid using a cloud cover threshold of 30% (see Figure 9). The percentage of cloud cover varied widely from year to year and from month to month. However, since 2000, and particularly for the months of January, February and April, we observed a decrease in the percentage of cloud cover, owing to the availability of two different Landsat sensors. For example, for the month of April (February) the percentage of cloudy months (cloud cover >30% of the total area) decreased from 62% (50%) for the period 1985-2000 to 28% (22%) for the period 2001-2018.

Winter Snow Cover Evolution
The snow cover summary referring to individual winter seasons in the period 1984-2018 (two examples in Figure 7) provides information on which pixels were classified as snow covered (i) in all the cloud-free images in that winter (high probability to observe snow cover); (ii) in none of the cloud-free images in that winter (low probability to observe snow cover). For each season we calculated the number of observations in those two classes (for which the probability to observe snow cover was high and low, respectively) and the corresponding area, expressing it as a percentage of the total GPNP area. The temporal evolution of the percent snow cover area, along with the cloud cover area, is shown in Figure 11c. The representativeness of those measures is displayed in Figure 11d, showing, for each winter season, the percent of pixels with three, two, one or no cloud-free observations available in the corresponding snow cover summary. For the first 15 winter seasons, cloud-free observations were generally not available for the three winter months considered (dark green bars), but rather on two (medium green bars) or one (light green bars) winter months only. Since 2000/01, the availability of cloud-free observations for the three winter months increased considerably, with almost half of the winter seasons having more than 50% of their observations being cloud-free for the three months. The uneven consistency of the data hampers a proper long-term trend analysis, for which more advanced techniques to reconstruct missing data are required. The information currently available seems to suggest, although with a high level of uncertainty, relatively large interannual variations in the percentage of the area with a high probability of observing snow (blue bars), varying from 58% to 98%, with the largest value being in winter 2005/06 and the smallest in 2015/16. There was a small negative, though not statistically significant, trend (decrease of about 0.5 km 2 /season, y = −0.0695x + 86.157 and R = 0.0044). Please note that for the trend analysis, only the seasons with less than 15% cloud cover were retained (seasons 1987/88, 1988/89, 1999/91 and 1995/96 exceeding this threshold have been excluded). The percentage of the area with a low probability of observing snow (red bars) was generally very small, varying between 0% and 18%, the latter being observed in 1991/92 and 2011/12 winters.
For completeness, winter mean temperature and (total and solid) precipitation time series from the HISTALP dataset are displayed in the upper panels of Figure 11a,b, showing no statistically-significant trend in the period 1985-2014 (see Table 1). Winter mean temperature varied from about 8 • C (during the 2009/10 season) to about −3 • C (during the 2006/07 season). Figure 11b highlights some relatively dry winters in 1991/92, 1992/93, 1994/95 and 1999/00 with less than 22 mm/month. Inversely, the 2003/04 winter season was characterized by large quantities of precipitation, exceeding 150 mm/month, most being snow. We jointly analyzed the snow cover information presented in the previous section with the climatic data from the HISTALP dataset, by quantifying the time correlation between the meteorological variables and the percentages of SCA ( Figure 12). The percentage of SCA with a high probability of observing snow was positively and significantly correlated with solid precipitation (R = 0.46), while the percentages of SCA with low probability of observing snow and solid precipitation were significantly anticorrelated (R = −0.54). The role of winter temperature on the extent snow is more complex. The percentage of SCA with a high probability to observe snow cover shows a weak, not statistically-significant, negative correlation with temperature; and the percentage of SCA with low probability to observe snow cover shows a weak, not statistically-significant, positive correlation with temperature. The weak correlation could be attributed to the fact that surface air temperature is often below the freezing point at such altitudes during winter, with limited effects on the precipitation phase. Table 4 provides a general summary of the statistics (minimum and maximum values, means and standard deviations) of the various variables, for DJF and for the months of February and April. For the minimum and maximum values also, the year of occurrence is reported in parentheses.
The same analysis performed for the winter season and shown in Figures 11 and 12 was repeated for the months of February and April, since they are key months to understand the dynamics and the changes of snow cover in mountainous areas. The results are displayed and commented on in Annex 5 ( Figures S4 and S5) and are briefly summarized in the following text. February's percentage of SCA exhibited a large interannual variability, with some years (1994, 2002 and 2006) in which almost the entire extent of the Gran Paradiso National Park appeared snow-covered, and others (1993 and 2000) in which the park appeared poorly snow-covered, according to the available satellite observations in that month. April's percentages of SCA and precipitation varied widely over the period 1984-2014 (Annex 5, Figure S5). April 1987 and 2006 were particularly snow-covered with more than 90% of the Gran Paradiso territory covered by snow. For both February and April, the analysis of the temporal percentage of SCA's evolution and trend was seriously hampered by the high number of years excluded from the analysis, owing to the presence of clouds. considered (dark green bars), but rather on two (medium green bars) or one (light green bars) winter months only. Since 2000/01, the availability of cloud-free observations for the three winter months increased considerably, with almost half of the winter seasons having more than 50% of their observations being cloud-free for the three months. The uneven consistency of the data hampers a proper long-term trend analysis, for which more advanced techniques to reconstruct missing data are required. The information currently available seems to suggest, although with a high level of uncertainty, relatively large interannual variations in the percentage of the area with a high probability of observing snow (blue bars), varying from 58% to 98%, with the largest value being in winter 2005/06 and the smallest in 2015/16. There was a small negative, though not statistically significant, trend (decrease of about 0.5 km 2 /season, y = −0.0695x + 86.157 and R = 0.0044). Please note that for the trend analysis, only the seasons with less than 15% cloud cover were retained (seasons 1987/88, 1988/89, 1999/91 and 1995/96 exceeding this threshold have been excluded). The percentage of the area with a low probability of observing snow (red bars) was generally very small, varying between 0% and 18%, the latter being observed in 1991/92 and 2011/12 winters. high probability to observe snow cover shows a weak, not statistically-significant, negative correlation with temperature; and the percentage of SCA with low probability to observe snow cover shows a weak, not statistically-significant, positive correlation with temperature. The weak correlation could be attributed to the fact that surface air temperature is often below the freezing point at such altitudes during winter, with limited effects on the precipitation phase. Figure 12. Correlation matrix between seasonal (DJF) snow extent and climate variables (temperature and solid precipitation). The shape, color and the orientation of the ellipse visually describe the correlation coefficients, whose value are also reported. The (*) indicate statistically significant correlations at a 95% confidence level (p-value < 0.05).

Figure 12.
Correlation matrix between seasonal (DJF) snow extent and climate variables (temperature and solid precipitation). The shape, color and the orientation of the ellipse visually describe the correlation coefficients, whose value are also reported. The (*) indicate statistically significant correlations at a 95% confidence level (p-value < 0.05).

Discussion
In this paper, we presented a methodology to exploit Landsat series satellite data stored in the Swiss Data Cube (SDC) for snow cover detection, applied to the Gran Paradiso National Park (GPNP). The frequency of available images and the cloud cover analysis over the analyzed 34 years showed that the availability of cloud-free observations remains a major limitation for mapping snow with optical remote sensing data. Indeed, almost half of the Landsat scenes available in the SDC for the period 1984-2018 in the GPNP are cloud-covered by more than 30% of the total area. Built in an EODC environment, the preliminary version of the SOfS algorithm implemented in the SDC can take advantage from the entire Landsat satellite series from 1984 to 2018 to combine different satellite scenes and generate maximum snow cover area products on a monthly time scale. Temporal filtering is the most common approach employed, particularly with Moderate Resolution Imaging Spectroradiometer (MODIS) satellites, to mitigate cloud contamination (i.e., [34,36,59]). In our study, the mosaicking method used to generate the monthly SCA products allowed to reduce cloud obstruction by approximately 7% and to gain about two additional monthly observations with cloud cover less than 30% of the total area. Besides the multi-day combination, we created a multi-sensor snow dataset combining observations from three different Landsat missions into a single dataset. This multi-sensor data fusion allowed us to go back further in time with the Landsat-5 mission, but the 2000s benefited from two Landsat missions at the same time (doubling the number of observations per month). By comparing the period 1984-2000, for which only observations from one satellite are available (two observations per month), with the period 2001-2018, for which the data are available from two satellites (four observations per month, Landsat-7 and Landsat-8 for the period 2013 to 2018), we observed that monthly cloud contamination was reduced by half after the 2000s. Using the potential of an EODC infrastructure to overcome the spatiotemporal constraints associated with the conventional single-sensor satellite mission, the combination of different sensors data into a unique dataset seemed to be an effective approach to increase the number of available/usable data. This new way of handling satellites data storage and analysis allows users to access the same ARD for different purposes and facilitates data analyses using a multitude of observations.
Monthly maximum snow cover products for December to February were aggregated to create "snow cover summaries" over the full 34 year-period and over each winter season. Snow cover summaries represented the percentages of snow cover observations with respect to the number of the total cloud-free observations in the considered period. Given that, we have to keep in mind that winter snow summaries are generally not representative of the three winter months for 34 years, particularly before the year 2000 when only the Landsat-5 sensor was active. However, the map showing the number of cloud-free observations available per pixel in the park over the full 34 year-period illustrates that cloud-free observations are relatively homogeneous, except in higher altitude areas where less cloud-free observations are found. As expected, the winter snow cover probability map shows that regions lying above 2000 m a.s.l. have high probability to be snow-covered in winter months. The probability to observe snow cover in winter is lowest in the south-eastern part of the park, where snowfall is also less abundant (from about 40 to 60 mm/month, see Figure 3l) compared to the rest of the park. Those are the areas lying at relatively low elevations. At higher elevations, lower air temperature and more abundant solid precipitation create favorable conditions for the formation and maintenance of snow cover (as shown in Supplementary Material, Annex 4). Using historical snow cover observations since 1984, our analysis highlights where snow cover is usually present and where it is seldom observed. Information on snow cover distribution and its evolution is essential for studying Alpine ecosystems, as well as highlighting the opportunities and risks for winter tourism's development. For example, according to a study carried out in Switzerland, it has been estimated that under warmer conditions, only ski areas at high elevations (above 2000 m a.s.l.) and with installations for producing artificial snow, will be able to open before Christmas time [76]. With more than one million people visiting the GPNP every year [38], information on the snow cover distribution during winter is essential to manage the local economy in a sustainable way.
The seasonal "snow cover summaries" produced in this study are based mostly on one or two monthly products per season, thus providing only partial information on the presence of snow during that period. As mentioned before, the month of December is more frequently covered by clouds than January and February, and thus it was less represented in our winter snow cover summary analysis. Despite this, winter SCA products showed some consistency with the climatic variables (particularly with solid precipitation) derived from the HISTALP dataset specifically developed for the alpine region. These preliminary results encourage further investigations using the proposed methodology.

Limitations
One of the major problems of snow detection using optical sensors, such as Landsat, is their inability to deliver surface information under cloudy conditions [77]. Our analysis has shown that, despite a monthly mosaicking methodology and a multi-sensors analysis, the cloud issue persists and limits the number of usable observations for snow cover mapping. On the monthly scale, abundant cloud cover strongly limits snow observations and the possibility to evaluate long-term trends and significant correlations between snow cover and climatic data.
The multi-day composite methodology tested in this study is subject to several other limitations. Using data from Landsat, which has a revisit time of 16-days, forced us to create a one-month composite with only two dates available, at least for the first 15 years (1984-2000). Previous studies have confirmed that the overall accuracy of multi-day combined snow cover products improves as the number of composition-days increase. In other words, the higher the number of images in the composite, the higher the overall accuracy. Similarly, the number of composition days was not homogeneous over the 34 years. Indeed, for the period 1984-2000 where only Landsat-5 was active, we used two images to calculate the maximum NDSI per month, while for the period 2001-2018, we use four images in a month to create the maximum NDSI products (with the combination of Landsat-5 and 7 and Landsat-7 and 8). Moreover, as mentioned by [78], the overall accuracy of a mosaicking method decreases when increasing the compositing time period window (one month here). Based on the assumption that snow will persist on the surface for a certain period of time [31], a one-month compositing period based on two or four images might be insufficient to assess snow cover change. It would thus be beneficial to exploit all satellites currently available in the SDC (Sentinel-1 and Sentinel-2, besides the Landsat series) and MODIS data to increase the number of remote sensing observations per month, and hence reduce the time window used in the mosaicking procedure. A previous study using MODIS observations over Austria has demonstrated that seven days of temporal filtering can remove more than 95% of clouds while maintaining an accuracy exceeding 92% [79].

Perspectives
In order to mitigate the influence of cloud cover, several methods have been proposed, such as spatial and temporal filtering or interpolation techniques [66,78,80,81]. Spatial filtering is designed to replace cloud-contaminated pixels by using information on surrounding adjacent pixels that are not obscured by clouds on the same date. According to the spatial position of the pixel, the estimation of the regional snow line (defined as the elevation above which all pixels are covered by snow) can help when reclassifying cloudy pixels as snow-covered or snow-free [66]. The study undertaken by Parajka et al. [79] using the daily MODIS snow cover products shows that the snow line method provides robust snow cover mapping, even if cloud cover is as large as 90%. In our case, with only Landsat-5 data available in the Swiss Data Cube from 1984 to 2000, that spatial filtering could be used to effectively reduce the impact of clouds. Similarly, as suggested by Dietz et al. [82], digital elevation model (DEM) information can be used to infer information on cloudy pixels. For example, Qobilov et al. [83] suggested that if a pixel is cloud covered, the nearest cloud-free pixel with the same elevation, azimuth and slope angle can be used instead. Another possibility to get information on cloud covered areas is to combine different optical datasets (e.g Landsat, Sentinel-2 and MODIS) and radar-based (Sentinel-1) sensors [84]. In fact, the synthetic aperture radar (SAR) from the first Sentinel satellite in the European Copernicus program, Sentinel-1, demonstrates an effective ability to map snow in all-weather and day-and-night conditions [10,15]. This will allow us to combine potentialities and limit their individual drawbacks. EODC seems to be an appropriate environment to conduct such multi-sensor data fusion and analysis and to improve the SOfS algorithm.

Conclusions
In this paper, we presented a first study on the snow cover evolution in the Gran Paradiso National Park, north-western Italian Alps, from 1984 to 2018, using the preliminary version of the snow observation from space (SOfS) algorithm benefiting from the Swiss Data Cube technology. The SOfS algorithm was applied to all Landsat series images (Landsat-5, Landsat-7 and Landsat-8) available in the period 1984-2018 in an analysis-ready form from the SDC, which corresponds to more than 480 Landsat scenes, to study the snow cover evolution in the GPNP.
The Earth observation data cube environment allowed to generate a monthly snow cover area time series in a totally automated-way for the Gran Paradiso National Park. The SOfS algorithm was employed for the first time in a multi-date (monthly composite) and a multi-sensor (combination of Landsat-5, 7 and 8 in the same dataset) mosaicking approach for increasing the number of cloud-free observations. Indeed, by using three different Landsat sensors incorporated into the same dataset, cloud cover contamination in the monthly SCA products after the 2000s was reduced by more than half.
However, despite the improvements deriving from our approach, the availability of cloud-free scenes remains a major limitation for mapping snow cover with optical remote sensing. There is a need to further investigate the detection of snow under cloud conditions using alternative approaches, such as spatial and temporal filtering techniques, and a combination of additional satellite data, such as Sentinel-1 and Sentinel-2. The integration and combination of consistent and comparable data provides new avenues and opportunities, increasing information density to minimize cloud contamination and to look further back in the past [52].
The preliminary version of the snow detection tool contained in the SDC and described in this paper needs to be further enhanced, in order to increase the quality, availability and accessibility of snow cover information across the Alps. With those improvements, the data cube technology might support a full analysis of the snow cover evolution, compared to a more traditional, scene-based sampling approach that might limit our ability to detect changes [26], allowing us better understand them and their consequences. As mentioned by Guo et al. [85] "multi-temporal Earth observation data may reveal large-scale processes and features that are not observable via traditional methods." Climate change must be taken as an opportunity to gain a more comprehensive understanding of the past and present evolution of our environment, for rethinking our mountainous regions from both environmental and economic perspectives through the establishment of adaptation measures to ensure the sustainable use of alpine resources [86].