Snow Cover and Snow Persistence Changes in the Mocho-Choshuenco Volcano (Southern Chile) Derived From 35 Years of Landsat Satellite Images

Mountain regions have experienced above-average warming in the 20th century and this trend is likely to continue. These accelerated temperature changes in alpine areas are causing reduced snowfall and changes in the timing of snowfall and melt. Snow is a critical component of alpine areas - it drives hibernation of animals, determines the length of the growing season for plants and the soil microbial composition. Thus, changes in snow patterns in mountain areas can have serious ecological consequences. Here we use 35 years of Landsat satellite images to study snow changes in the Mocho-Choshuenco Volcano in the Southern Andes of Chile. Landsat images have 30 m pixel resolution and a revisit period of 16 days. We calculated the total snow area in cloud-free Landsat scenes and the snow frequency per pixel, here called “snow persistence” for different periods and seasons. Permanent snow cover in summer was stable over a period of 30 years and decreased below 20 km2 from 2014 onward at middle elevations (1,530–2,000 m a.s.l.). This is confirmed by negative changes in snow persistence detected at the pixel level, concentrated in this altitudinal belt in summer and also in autumn. In winter and spring, negative changes in snow persistence are concentrated at lower elevations (1,200–1,530 m a.s.l.). Considering the snow persistence of the 1984–1990 period as a reference, the last period (2015–2019) is experiencing a −5.75 km2 reduction of permanent snow area (snow persistence > 95%) in summer, −8.75 km2 in autumn, −42.40 km2 in winter, and −18.23 km2 in spring. While permanent snow at the high elevational belt (>2,000 m a.s.l.) has not changed through the years, snow that used to be permanent in the middle elevational belt has become seasonal. In this study, we use a probabilistic snow persistence approach for identifying areas of snow reduction and potential changes in alpine vegetation. This approach permits a more efficient use of remote sensing data, increasing by three times the amount of usable scenes by including images with spatial gaps. Furthermore, we explore some ecological questions regarding alpine ecosystems that this method may help address in a global warming scenario.

Mountain regions have experienced above-average warming in the 20th century and this trend is likely to continue. These accelerated temperature changes in alpine areas are causing reduced snowfall and changes in the timing of snowfall and melt. Snow is a critical component of alpine areas -it drives hibernation of animals, determines the length of the growing season for plants and the soil microbial composition. Thus, changes in snow patterns in mountain areas can have serious ecological consequences. Here we use 35 years of Landsat satellite images to study snow changes in the Mocho-Choshuenco Volcano in the Southern Andes of Chile. Landsat images have 30 m pixel resolution and a revisit period of 16 days. We calculated the total snow area in cloud-free Landsat scenes and the snow frequency per pixel, here called "snow persistence" for different periods and seasons. Permanent snow cover in summer was stable over a period of 30 years and decreased below 20 km 2 from 2014 onward at middle elevations (1,530-2,000 m a.s.l.). This is confirmed by negative changes in snow persistence detected at the pixel level, concentrated in this altitudinal belt in summer and also in autumn. In winter and spring, negative changes in snow persistence are concentrated at lower elevations (1,200-1,530 m a.s.l.). Considering the snow persistence of the 1984-1990 period as a reference, the last period (2015-2019) is experiencing a −5.75 km 2 reduction of permanent snow area (snow persistence > 95%) in summer, −8.75 km 2 in autumn, −42.40 km 2 in winter, and −18.23 km 2 in spring. While permanent snow at the high elevational belt (>2,000 m a.s.l.) has not changed through the years, snow that used to be permanent in the middle elevational belt has become seasonal. In this study, we use a probabilistic snow persistence approach for identifying areas of snow reduction and potential changes in alpine vegetation. This

INTRODUCTION
In addition to increasing average temperatures globally, climate change is also disrupting and altering seasonal patterns that drive ecosystem function (IPCC, 2021). There is growing evidence that the rate of warming is amplified at higher elevations, such that high-mountain environments experience more rapid changes in temperature than environments at lower elevations (Beniston et al., 1997). Indeed, mountain regions have been identified as having experienced above-average warming in the 20th century, a trend likely to continue (Beniston et al., 1997;Liu and Chen, 2000). These accelerated changes in alpine areas means that ecosystems are being impacted by both earlier snowmelt and reductions in snowfall. Seasonal snow makes up the largest portion of the cryosphere with up to 31% of the earth's surface being covered by seasonal snow (Vaughan et al., 2013). However, globally snow cover has decreased by 10% compared with snow coverage in the late 1960's (Walther et al., 2002) and since 1972 snow cover duration in the Northern hemisphere has been declining at a rate of 5 days/decade, attributable to earlier spring snowmelt (Choi et al., 2010). A recent global assessment of mountain regions reported that around 78% of the world's mountain areas have experienced a snow cover decline between the 2000 and 2018 study period (Notarnicola, 2020).
Snow is a critical component of alpine areas. It drives hibernation of animals and determines the length of the growing season for plants. In addition, snow provides an insulating layer for plants and animals, protecting them from temperature extremes during winter (Callaghan et al., 2011;Pauli et al., 2013). Warming in alpine areas will not only reduce snowfall but also advance the timing of snowmelt. Ecologically, this timing is particularly important for the plants and animals that rely on the insulation provided by the snow to cope with freezing temperatures (Wipf et al., 2009;Callaghan et al., 2011;Pauli et al., 2013;Wheeler et al., 2014). Early snowmelt has been found to dramatically change soil microbial composition (Zinger et al., 2009), expose plants to freezing events they cannot tolerate (Bokhorst et al., 2009;Gerdol et al., 2013) and advance plant reproduction and growth (Lambert et al., 2010;Livensperger et al., 2016). Hence, a change in the time of the snowmelt, as well as a reduction in snowfall, can have significant consequences for high elevation ecosystems. Furthermore, snow cover and duration also influences lowland ecosystems, as the decrease in snow fall can decrease freshwater for downstream ecosystems. Snow water equivalent, a parameter that takes into account snow depth and density, has a direct influence on the hydrologic cycle and water supply in snowmelt-dominated regions (Barnett et al., 2005;Gan et al., 2013). Shrinkage of spring snow water equivalents over the last few decades has been reported for the Arctic region (Brown et al., 2010), East Europe (Bednorz, 2004), and the Northern Hemisphere (Brown, 2000;Brown and Robinson, 2011).
In spite of the crucial role played by snow in alpine areas, there are gaps in the knowledge as to how alpine snow is responding to a changing climate. It is still unclear whether all alpine areas will respond in a similar fashion or if differences exist between hemispheres, climate types, and environments. A search on the "Web of Science" with the keywords: snow, snow duration, mountains, climate change, and warming shows 60 studies ranging from 1999 till 2020 that specifically investigated how snow regimes were changing in alpine areas. Of these 60 studies, 39 were done in the Northern hemisphere, including 21 studies undertaken in Europe. In the Southern hemisphere, 21 studies were undertaken and 15 of these were based in South America (two of these studies also included North American sites), with 11 studies along the Andes Cordillera, particularly the Northern and Central Andes. Fourteen of these 21 Southern hemisphere studies used satellite images to quantify snow cover. Most of the South American studies describe glacier dynamics (Bodin et al., 2010;Durán-Alarcón et al., 2015), measure snow effects on albedo (Lhermitte et al., 2014) or validate methodologies to describe snow patterns (Soncini-Sessa and Volta, 2005). Surprisingly, only four studies focused particularly on the effect of climate change on snow accumulation (Masiokas et al., 2012;Migliavacca et al., 2015;López-Moreno et al., 2017;Notarnicola, 2020). The Andes Cordillera maintains both permanent and seasonal snowpack (Masiokas et al., 2020), with a seasonal snowpack coverage of approximately 61,000 km 2 Saavedra et al., 2017). The seasonal snowmelt from the Andes provides the majority of the water required for human consumption, agriculture, industrial use, electricity generation, and aquifer recharge (Masiokas et al., 2006(Masiokas et al., , 2012(Masiokas et al., , 2020. Therefore, understanding how the snowpack duration in these mountains will be impacted by a changing climate is extremely important. Given this, studies of the effects of climate change on snow cover and duration in the Andes are starting to emerge. Existing studies have mainly used data on snowfall and accumulation extracted from the closest weather stations. Thus, this kind of data does not consider spatial variation in snowfall or duration of snow, which is known to vary over a short spatial scale in the mountains. To overcome the spatial and temporal constraints imposed by the use of weather station data, remote sensing appears as a cost-efficient approach to measure and monitor snow cover systematically over large areas. The main limitations are the availability of cloud free images, especially when we want to study seasonal patterns. Besides measuring snow cover in cloudfree images, in this study, we used a per-pixel approach to take advantage of incomplete images due to cloud cover and other causes of data loss (e.g., the SLC instrument failure of the Landsat 7 ETM+ sensor). A similar approach was used by Macander et al. (2015) to study snow persistence in Alaska, although the classification algorithm they proposed returns maps of snowfree dates rather than snow persistence. Here we use Landsat satellite images from 1984 till 2019, including cloud free and incomplete scenes, to study the spatio-temporal changes in both snow cover and snow persistence in the Mocho-Choshuenco Volcano in the Southern Andes of Chile. Besides snow cover per date, traditionally studied from available cloud-free remote sensing data, here we propose to calculate a pixel based snow frequency measure per season (e.g., the number of times that a pixel is covered by snow considering the available satellite records), from now on called "snow persistence." This provides robust snow persistence maps per season and per different time periods. Pixels showing snow persistence of 100% during the four seasons will show the areas with permanent snow cover. Changes in snow persistence, for example from one decade to the next one, can show a reduction of the area with permanent snow (less pixels with 100% snow persistence) or highlight the areas where the permanent snow changed to seasonal snow (e.g., 50 or 75% snow persistence) or ephemeral snow (e.g., <25%). In the present study, we derived two remote sensing metrics (snow cover area and snow persistence) to address the following questions regarding temporal variations in snow cover: Is the total area covered by snow changing through the years? What is the seasonality of the total snow cover in the volcano? and questions looking at the spatial and temporal variation in snow persistence: Are there temporal differences in snow seasonality at different elevational belts? What are the spatial patterns of snow persistence in the different seasons? Is there a reduction of snow persistence in the different seasons through the years? And finally, which elevational belts are changing the most in regard to snow persistence? We expect the total area covered by snow in this volcano to be reduced through the years as seen in other alpine areas around the world and we expect mid and low elevation belts to have a more pronounced reduction of snow persistence particularly in Spring and Autumn, indicating early snowmelt and late snowfall as has been predicted for alpine areas.

Study Area
We studied snow cover and its persistence through time in the Mocho-Choshuenco volcano complex, located in Los Ríos Region in the South of Chile (Figure 1). The Mocho-Choshuenco is an active stratovolcano compound located in Los Ríos Region, inside the Southern Chile Volcanic Zone (Rawson et al., 2015(Rawson et al., , 2016. The volcanic compound covers an area of 250 km 2 , an approximate volume of 100 km 3 with a summit of 2,422 m a.s.l. (Moreno and Lara, 2007). Currently, the compound is 7th in the specific risk ranking developed by National Service of Geology and Mining of Chile. Most of the study area belongs to the Mocho-Choshuenco National Reserve under the administration of the Chilean National Forest Service CONAF. We selected the area above the treeline (above 1,200 m.a.s.l.) equivalent to the alpine area to analyse the snow duration. The alpine area in this volcano is dominated by shrubs such as Adesmia longipes, Escallonia alpina, Gaultheria poepigii, Ovidea andina, and Senecio chionophilus and grasses such as Anthoxanthun juncifolium, Festuca thermarum, and Poa obvallata. Toward higher elevations small perennial herbs like Nassauvia revoluta and Nassauvia lagascae sparsely inhabit the scoria (Teillier et al., 2012).

Landsat Images Availability and Definition of Periods of Analysis
All available Landsat images, including Landsat 5, 7, and 8, were accessed and analyzed in the cloud using the Google Earth Engine capabilities (Gorelick et al., 2017). The specific code of this and other sections of the manuscript are available in this link. 1 Specifically, we used the atmospherically corrected product Landsat Surface Reflectance Collection 1 Tier 1 corresponding to the highest radiometric and precision terrain quality. All Landsat scenes have 30 m pixel resolution. A total of 678 scenes from Path 232/233 and Row 88 were used, from which 23.5% of the scenes were completely cloud free and 76.5% were incomplete (with data gaps). Incomplete images are the result of partial cloud cover or sensor data loss, e.g., failure of the SLC instrument on Landsat 7 (Supplementary Material and Supplementary Figure 1). Cloud-free images were used for calculating total snow cover area per date of the study area and both cloud-free and incomplete images were used for snow persistence calculations which are applied on a pixel basis. Before 1998, Landsat data was scarce with a temporal gap of 7 years between 1990 and 1997 (Supplementary Material and Supplementary Figure 1A). From 1998 on, the dataset provides images consistently over time with 25-50 images per year, considering the revisit time of 16 days that each satellite has, but also that the different satellites overlapped in time, e. g., Landsat 5 and 7 between 1999 and 2013 and Landsat 7 and 8 between 2013 and 2020. Following the meteorological conditions of the South Hemisphere, there is a marked seasonality in data availability with a minimum of data during winter, particularly in June, due to increased cloud cover. In order to have sufficient data points for snow frequency calculations (snow persistence maps), we grouped the Landsat data by season: summer (Dec-Jan-Feb), autumn (Mar-Apr-May), winter (Jun-Jul-Aug), and spring (Sep-Oct-Nov) and by five periods of time: P1 (1984-1990), P2 (1998-2003) Total Snow Area Estimation and Validation Using the NDSI Snow detection was automatically performed using the Normalized Difference Snow Index (NDSI) proposed by Dozier (1984) and the spectral bands of Landsat images, specifically the short-wave infrared (SWIR) and green (G) bands (Figure 2). This snow index has been widely used for snow detection in different regions of the world (Burns and Nolin, 2014;Zhang et al., 2019; 1 https://github.com/JoseLastra/Frontiers_paper Masiokas et al., 2020). The formula of the NDSI is given as follows: From all original Landsat images (Figures 2A-D), where the snow is already recognizable by simple visual inspection in different seasons (A: summer, D: winter), the NDSI images were retrieved (Figures 2B,E). To classify pixels as "snow" or "not-snow" several NDSI thresholds have been proposed (Figures 2C,F). In the literature, thresholds of 0.4 (or similar) have been used or recommended (Hall and Riggs, 2007;Sankey et al., 2015;Masiokas et al., 2020) but this threshold has shown large standard deviations across different areas (Härer et al., 2018). In this study, we carried out a calibration exercise to define a NDSI threshold for the Mocho-Choshuenco volcano and later a validation of the snow mapping exercise. To achieve this, we set two transects (Figure 2) in a gradient from permanent snow (top of the volcano) to ephemeral snow (downward the volcano). In the west-east transect we defined four calibration points (C1-C4) and in the east-west transect four validation points (V1-V4  Figure 1), we visually inspected whether the corresponding pixel was "snow" or "not-snow" using a false color composition (SWIR-NIR-Red) of the original Landsat images (Figures 2B,E) supported by available Google Earth high spatial resolution images, especially in the transition between snow and not snow. With this database we were able to study the dynamic range of the NDSI in a gradient of permanent snow to ephemeral snow and across the seasons and years. By constructing boxplots of NDSI for pixels with "snow" and "not snow" of the four calibration points, we identified the optimal threshold to detect snow. Finally, we constructed a confusion matrix and calculated the overall accuracy of the defined NDSI threshold using the validation dataset. By using this NDSI threshold, binary snow and not-snow maps were constructed from which the total snow area per date was calculated (Figures 2C,F). Additionally, we constructed a confusion matrix of the snow flag included with the Landsat Surface Reflectance product level 2.0 to compare the results obtained from NDSI thresholding and the default Landsat snow flag. Visualizations of the annual cycle of total snow area for the entire volcano and also at three different elevational belts (1,200-1,530, 1,530-2,000, and >2,000 m a.s.l.) were constructed using kernel density estimations of the snow area -time space. For this purpose, we used the PhenKplot function from the "npphen" R function (Chávez et al., 2017).

Snow Persistence Modeling and Maps
As indicated in previous sections, cloud free satellite images are scarce at high latitudes and in mountainous areas, limiting remote sensing based assessments of snow cover dynamics. For example, in our study site, only 23% of the 678 Landsat images were cloud-free. To make use of incomplete satellite images (e.g., the 77% of our Landsat dataset) we use a complementary approach to total snow cover area calculation (Section "Total Snow Area Estimation and Validation Using the NDSI") which we call "snow persistence". Snow persistence is the probability of the presence of snow for a given pixel at a given day of the year (DOY) or through a time interval. Frequencies are calculated per pixel and per DOY based on all available Landsat data (cloud free and incomplete images) and can also be calculated for different periods of time (e.g., the five periods defined in this study, see Supplementary Material and Supplementary Figure 1). Based on the "snow"/"notsnow" binary images (Figures 2C,F) we extract for each pixel the time series of snow presence per date. In order to have more binary observations for the fitting process, the data can be aggregated by month or season. Then we fit a binomial generalized linear model to the "snow"/"not-snow" data per month or season. Different from the study of Macander et al. (2015), where the classification algorithm returns the snow-free date per pixel, in this study, the response variable corresponds to the probability that a given pixel is covered by snow at a given date or time interval. Here we calculated snow persistence per season for the whole Landsat time series and also for the five periods of analysis. Pixels with snow persistence of 100% in summer, winter, autumn, and spring can be considered as permanent snow. Changes in snow persistence through time can demonstrate the potential impacts of climate change: e.g., areas where snow persistence in summer was 100% and now is 80% indicates permanent snow is becoming seasonal or ephemeral. Since this is calculated at a pixel basis, there is no need for complete cloudfree images and thus all valid pixels of all available satellite images can be used. We developed a function in R Core Team. (2020) to calculate snow persistence maps for the entire study area and also at three different elevational belts (1,200-1,530, 1,530-2,000, and >2,000 m a.s.l.) using the capabilities of the "raster" package (Hijmans, 2019) to handle large Landsat rasterstacks and also R capabilities to take advantage of multi-core processing. To evaluate the quality of the binomial GLM models, we calculate the Area Under the ROC Curve (AUC-ROC) per pixel. AUC values close to 0.5 means the classifier is no better than a random classifier, and values close to 1 means the classifier has a high quality performance. We finally constructed AUC-ROC maps per period, showing the spatial distribution of the model performances in terms of AUC-ROC grouped in four categories suggested by Hosmer et al. (2013) which are <0.7 (poor), 0.7-0.8 (acceptable), 0.8-0.9 (excellent), and 0.9-1.0 (outstanding).

Total Snow Cover Area and Seasonality
The boxplots obtained from the calibration exercise showed that a threshold of NDSI of 0.5 maximizes classification accuracy (Supplementary Material and Supplementary Figure 2). Using this threshold to classify snow and not-snow for the validation dataset, we achieved an overall accuracy of 98.27% (Supplementary Material and Supplementary Table 1). The NDSI > 0.5 threshold showed higher overall accuracy than the default Landsat snow flag (98.27% against 90.88%) and lower omission error for the snow class (2.8% against 15%) (Supplementary Material and Supplementary Tables 1, 2). Hereafter, all calculations were performed using this threshold.
The total snow cover area of our study site during the 1986-2019 time frame is presented in Figure 3. This is the longest time series available in remote sensing records at this spatial scale (30 m). As expected, there is a strong seasonality in total snow cover with a total area reaching between 80 and 105 km 2 in winter and a minimum and consistent area of about 18-20 km 2 in summer. In Figure 3B, we show a kernel density estimation of the total snow area -DOY space using all annual cycles of snow area. This way, the long-term snow seasonality of the Mocho-Choshuenco volcano is deployed. The dark-red line shows the most likely snow area at different DOY's along to its frequency distribution (colored area). From this figure, we observed that around DOY 150 (end of May) the snow starts to fall and there is a sharp increase in the area covered by snow, reaching its maximum by the end of the winter (August, about DOY 240). A wide range of total snow is observed at the peak of the winter season demonstrated by the broad frequency distribution around the most probable value of about 100 km 2 . In spring, the melting process begins and therefore this period has the most variation in snow cover through time. By the DOY 300 (end of October) the melting process ends and the total snow area goes down to about 40-30 km 2 , and keeps going down to the basal summer snow cover of 20 km 2 . This minimum area covered by snow in summer was relatively stable over 25 years until 2014 when it decreased and remained below 20 km 2 until 2019 (basal red line in Figure 3A). Besides this decrease in the summer cover area, there is no other clear evidence of total snow area changes through time.
When we decomposed the total snow area time series into three different elevational belts (Figure 4), we observed clear seasonal variation in the area covered by snow at the low and middle elevations, but not at high elevation (>2,000 m a.s.l.). In the middle elevation, the snow area in winter remains more stable (about 40 km 2 ) and for longer (about 170 days) compared with low elevation (40-55 km 2 lasting about 100 days). From this Figure, we can see that the decrease in total snow area in the Mocho-Choshuenco volcano, starting in 2014, takes place mostly at middle elevations ( Figure 4C).

Snow Persistence Spatio-Temporal Patterns
The AUC-ROC maps calculated for the five different periods (Supplementary Material and Supplementary Figure 3), show that the pixel based models used for snow persistence mapping present high efficiency for practically the entire study area. Supplementary Table 3 in the Supplementary Material, summarizes the area at different levels of model efficiency, showing that for the five periods less than 0.5% of the study area presented poor efficiency (AUC-ROC < 0.7). Furthermore, less than 5% of the study area presented "acceptable" efficiency for all periods while the categories "excellent" and "outstanding" were the most represented. Overall, these results suggest a good performance of the snow persistence modeling and mapping for all periods.
The snow persistence maps calculated with 29 years of Landsat data (678 satellite images) show distinctive spatial features for the four different seasons (Supplementary Material and  Supplementary Figure 4). The area with highly persistent snow (i.e., snow persistence > 75%) is much larger in winter and spring than in summer and autumn. Consequently, the low elevational belt is snow-free in summer and autumn, while in the middle elevational belt, some areas remain with highly persistent snow, Seasonality of total snow cover area was calculated using Gaussian kernel density estimations of the snow cover along the year based on all historical records. Red color in panel (B) indicates high frequency of a given snow cover at different days of the year and yellow color lower frequencies of snow cover areas. Summer is very stable and the area covered by snow is about 20 km 2 . The end of winter and the start of spring are the periods with larger snow cover, approximately 100 km 2 . Then the snow cover decreases gradually in spring till the basal 20 km 2 of snow cover in summer. Note that from 2014 onward the lowest snow cover area is lower than the basal 20 km 2 (red dotted line). especially along the East-South slope of the volcano. During winter and spring this middle elevational belt is permanently covered by snow. Finally, the high elevational ring is covered by snow practically all year long, besides some specific areas with high slope near the crater and on the northern slope, which are snow-free during summer and autumn.
The total area with different snow persistence ranges observed for the five different periods at three elevational belts is presented in Figure 5 while the changes per period considering all elevational belts together is summarized in the Supplementary Material and Supplementary Figure 5. Pixels with snow persistence > 95% in summer can be considered as permanent snow and when setting the first period (with 20 km 2 ) as a baseline, our results show that permanent snow area has decreased by 17% in period 2, then recovered 4% in period 3 with an accumulated decrease of 13%, then decreased again by 20% in period 4 to finally reach its minimum extent in the last period with 5.75 km 2 (28% reduction). Temporal changes in autumn and spring are variable. In winter we observed a sustained reduction of the area with snow persistence > 95% in periods 2 (24%), 3 (56%), 4 (25%), and 5 (42%) with respect to period 1 (Figure 5, Supplementary Material, and SupplementaryFigure 5).
When subtracting the snow persistence maps of the different periods from the map of period 1 (Figure 6), we can identify specific locations with a gain (light blue to blue colors) or a reduction (yellow to red colors) in snow persistence. Changes in snow persistence through the periods are seen for summer and autumn at the middle elevation belt (Figure 6). In spring and winter, the lower elevational belt is decreasing in snow persistence especially after 2010. The affected area in the middle elevation belt is smaller and less drastic than in autumn. Extreme changes in the persistence of snow are seen for low elevations in winter.

DISCUSSION
In this study, we have investigated the spatial and seasonal variation of snow cover and snow persistence in the Mocho-Choshuenco volcano in Southern Chile using a long record (35 years) of Landsat Satellite images. Snow cover in mountain areas has been traditionally assessed using remote sensing data and the NDSI from global (Notarnicola, 2020) to regional scale, e.g., the South American Andes (Masiokas et al., 2020), using MODIS products (250 m resolution), and to local scale (Burns and Nolin, 2014) using Landsat data (30 m resolution). Data loss is huge, since all images with a partial coverage of the study area have to be neglected. MODIS 8-day composites cope partially with this problem, but still for the 8 day temporal window the scene has to be complete in order to be used. In this manuscript, in addition to the traditional snow cover calculation, we introduced a different metric: seasonal snow persistence calculated at a pixel basis, allowing us to take full advantage of Landsat images both complete cloud-free (159 scenes) and incomplete scenes (519). Using this approach we can construct robust snow persistence change maps for different periods and across seasons with 326% more data than the traditional approach. The AUC-ROC maps for the five different periods showed a good performance of the snow persistence models for almost the entire study area with less than 0.5% of the study area presenting a "poor" performance (AUC-ROC < 0.7). We have found a clear seasonal variation in snow cover and persistence for the whole area above the treeline, as well as across different elevational belts, with the low elevation belt being the most variable in terms of snow cover and persistence. Comparing across periods, the areas of high snow persistence (75-95 and >95%) are the most variable across periods. Indeed, a clear reduction of the area of high snow persistence is seen in summer (28% reduction in the period 2015-2019 with respect to the base period [1984][1985][1986][1987][1988][1989][1990]. Evidence of reduction of snow persistence was also found when comparing periods for the middle and low elevational belt. The middle elevational belt is decreasing in snow persistence in summer and autumn while a decrease in snow persistence in winter and spring is seen in the lowest elevational belt. In addition, we found that the minimum snow cover area is decreasing through the years, particularly in the middle elevation belt. Here we compare our findings with the ones found in other mountain areas. We explore the ecological relevance of our findings and highlight the advantages and constraints of this methodology. Finally we discuss the potential ecological questions that can be addressed with this method.

Contrasting Snow Seasonality Between Elevational Belts in the Mocho-Choshuenco Volcano
The kernel density estimation applied to the snow area along the year allows us to measure the long-term snow seasonality for the whole volcano and its temporal variation, given by the frequency distribution along the most expected snow area along the year. This non-parametric approach has been proposed and used to study vegetation phenology and disturbances (Chávez et al., 2019a,b,c;Decuyper et al., 2020) and it is a useful approximation to assess the seasonality of snow cover. The seasonality of the area covered by snow varies strikingly between elevational belts. The snow lasts approximately 70 days more in the middle elevation belt compared to the low elevation belt, this is a big difference in snow duration over short elevational scale, that accentuates the environmental differences between elevations. Mountains are usually viewed and studied in terms of the air temperature decreasing across the elevational gradient. However, snow decouples plants from air temperature and thus paradoxically, plants at lower elevations are exposed to colder temperatures than higher elevation ones. Considering the effects of climate change on the seasonality of snow duration is crucial to understanding how climate change will impact alpine vegetation.

Snow Cover and Persistence Reduction Through the Years in Mocho-Choshuenco Volcano
We have found evidence of snow cover and persistence reduction in the Mocho-Choshuenco volcano. Firstly, using all the snow free scenes, a decrease in the minimum area covered by snow can be seen since 2014 onward. This was also evident when studying the seasonality of the snow cover by elevational belts, with the middle elevational belt showing the decrease in minimum snow cover. The maps of difference in snow persistence per season also show a striking variation in the area of snow persistence. The middle elevational belt is decreasing in snow persistence in summer and autumn while a decrease in snow persistence in winter and spring is seen in the lowest elevational belt. The highest elevational belt is not experiencing any changes in snow persistence with the years. A range of studies have identified that the most rapid and consistent losses of snow, in terms of depth and duration, are mid-elevation areas (e.g., subalpine zone) and those with Mediterranean/maritime climates, where mean air temperatures are close to freezing (Brown and Mote, 2009;Steger et al., 2013). In the Tibetan Plateau, areas below 3,500 m a.s.l. experienced shorter snow cover duration in a 15 years study using MODIS (Wang et al., 2017). A review of the status of the European cryosphere found a general negative trend for snow depth and snow water equivalents below 2,000 m a.s.l. and no trend at higher elevations in the Alps (Beniston et al., 2018). For Western United States, Mote et al. (2018) found that 33% of the monitored sites showed a significant decline in snowpack, these changes occur mainly in spring and at milder elevations. Recently, a study looking at snow cover in mountains around the globe, found that 78% of the global mountains are undergoing a decline in snow cover duration of up to 43 days and a decrease of snow cover area up to 13%. In this study, the area of the Andes mountains considered (between 42 • S and 29 • S) exhibited the most drastic change in snow cover, with a decrease of snow cover duration of −26.6 days between 2,500 and 4,000 m a.s.l. (Notarnicola, 2020), similar results were found by Saavedra et al. (2018) and Malmros et al. (2018). Our results show reductions on snow persistence in every season. Indeed, besides spring, the area with snow persistence > 75% is the lowest in the last of the five periods, and excluding winter, the area with snow persistence > 95% is also the lowest in the last period. The effects of a reduction in snow persistence on alpine vegetation at different elevations and for different seasons is of critical relevance if we are to make predictions of the effects of climate change on alpine ecosystems.
Climate predictions for the Andes are extremely complex, due to the sharp vertical rise of the Andes and the length of the Cordillera covering several climatic regions and the influences of additional climate systems including El Niño patterns (Kohler and Maselli, 2009). Therefore, studies on changes in snow persistence here need to encompass several sites across this mountain chain and, for example, findings on the Southern Andes cannot be extrapolated to what might occur in the tropical, Mediterranean or Patagonian Andes. Studies covering most of the climatic regions of the mountain range are needed to understand the local effects of climate change or climate systems on snow persistence in the Andes.

Potential Ecological Effects of Reduction of Snow Persistence
There is clear evidence that snow duration is changing in mountains worldwide, however, the question of how plants will respond to changes in snow conditions is still unclear and most of the studies looking at the ecological responses to variation in seasonal snow cover are done in the northern hemisphere (Slatyer et al., 2021). Advanced snowmelt could potentially increase plant fitness by prolonging the growth period and hence increasing resource allocation to seed production (Galen and Stanton, 1993). However, there is increasing evidence that early snowmelt also increases the risk of frost damage in adult life stages. We already know that naturally occurring warm spells and early snowmelt increase the mortality of flower buds due to frost damage (Molau, 1993;Inouye, 2008;Lambert et al., 2010). Further, experimental snow reduction and advancing snowmelt have been linked to reductions in biomass and fruit production, again resulting from frost damage (Gerdol et al., 2013). The effects of reduced snow persistence will depend on the season and elevation where this reduction occurs, with perhaps more detrimental effects for winter and for mid and high elevations areas where plants are acclimated to late snowmelt. Therefore, methods that allow us to look at the spatial and temporal variation in snow persistence are extremely relevant.
Snow duration determines species composition in the mountains, therefore, changes in snow conditions could promote changes in species composition and modify the dominant plant life-form (e.g., Hagedorn et al., 2014;Elliott and Petruccelli, 2018), and even impact on other trophic levels (Roland et al., 2021). In addition, a potential slow upward colonization of shrubs over areas currently occupied by perennial herbs might occur as a consequence of early snow melt in these areas. A second potential effect could be observed in the distribution of the species along the elevational gradient. In several alpine environments it has been observed a significant change in the optimal elevation of the plant species distribution. Similar changes have been observed in Western Europe (Lenoir et al., 2008). In these cases, elevational changes in the maximum value for some fitness metrics (e.g., seed production, seed germination, and growth, etc.) could be expected. The effects of early snow melt on plant fitness in alpine areas requires urgent assessment to provide key information for the conservation of these environments. Using this methodology we have found reduction of snow duration at a temporal and spatial scale in the Mocho-Choshuenco volcano. This method could be used to study the influence of snow duration on plant reproduction, phenology, seed germination, species persistence, and stress tolerance in environments where snow is a major feature. In addition, this methodology can be used to find priority areas for conservation in the mountains, such as areas where a strong pattern of snow reduction is seen or areas where plant species exist without significant potential for adaptation.

Advantages and Disadvantages of Seasonal Snow Persistence Calculations
We have already highlighted that our approach to the study of snow persistence optimizes the use of remote sensing data, including scenes with spatio-temporal gaps due to cloud cover or sensor malfunction. It can be applied to different remote sensing data such as MODIS snow products (250 m), Landsat (30 m), or Sentinel 2AB data (20 m). In this case, we were able to use three times more Landsat data than a traditional cloud-free based assessment. This allowed us to provide a robust description of snow spatio-temporal patterns and changes even for the cloudy Seasons. The main limitation of this approach is that it does not allow monitoring of snow cover on a short-term temporal basis, since several observations are needed to calculate a robust snow persistence. For example, it is difficult to calculate snow persistence on a monthly basis instead of seasonal as we did here, especially during cloudy seasons. Additionally, we advise to use this methodology for areas above the treeline (i.e., alpine areas), given that taller vegetation below the treeline will likely interfere with the snow visualization. In this study, we took careful consideration of this limitation when defining our study area: above the treeline (1,200 m.a.s.l.). Snow duration studies in the field are extremely laborious and of short span. We hope this methodology helps to understand how climate change is affecting snow cover and persistence patterns worldwide, especially for cloudy areas or seasons with patchy satellite data, and the ecological consequences of it.

DATA AVAILABILITY STATEMENT
The data and scripts supporting the conclusions of this article are made available by the authors in this link https://github.com/ JoseLastra/Frontiers_paper, without undue reservation.

AUTHOR CONTRIBUTIONS
RC and VB conceived the research article. RC, SE, and JL developed the R code for the remote sensing analyses. RC, VB, SE, and DH-P prepared and edited the manuscript. JL and RC conceived and prepared the figures. All authors have revised and approved the manuscript.