Snow Cover Change as a Climate Indicator in Brunswick Peninsula, Patagonia

Snow cover changes are assessed for the Brunswick Peninsula in southern Patagonia (52.9oS to 53.5oS), located on the transition between the wet Pacific Ocean area and the drier leeward side of the Andes. We use the Normalized Difference Snow Index (NDSI) and a new index which we call SnowPower, combining the NDSI and the Melt Area Detection Index (MADI), to reconstruct the snow cover extent and its temporal distribution for the period 2000-2016, based on Moderate Resolution Imaging Spectroradiometer Sensor (MODIS) MOD09GA products. Comparison of these satellite products with available snow duration and snow height data for 2010-2016 at the local Club Andino ski centre of Punta Arenas shows that the NDSI exhibits the best results. A reasonable and significant linear correlation is found between the MODIS NDSI snow cover extent and the mean monthly temperature at Punta Arenas Airport combined with the monthly snow accumulation at Jorge Schythe station at Punta Arenas city for the extended winter period (April to September) from 2000 to 2016. Snow cover changes within this time series are extended to 1972 and 1958 based on historical climate data of Jorge Schythe and Punta Arenas airport, repectively. The results show a significant decreasing trend of snow extent of 19% for Brunswick Peninsula for the 45-year period (1972-2016), which can be attributed to a statistically significant long-term warming of 0.71oC at Punta Arenas during the extended winter (April to September) in the same period. Multiple correlation with different climate variables indicates that solid precipitation has a relevant role on short-term snow cover variability, but is not related to the observed long-term snow cover decrease.


INTRODUCTION AND REGIONAL SETTING
Climate change particularly affects snow-covered areas (Hughes, 2003;IPCC, 2013), since even small changes in the surface-atmospheric energy balance may cause significant decreases in snow cover, affecting the local hydrology, associated ecosystems, hydroelectric power generation, and tourism activities (Gobiet et al., 2014;McGowan et al., 2018). Therefore, if climate trends can be detected, valuable information can be provided to decision makers and stakeholeders. Snow cover is a critical component of the cryosphere and climate system at both local and global scales (Chen et al., 2018). Snow also affects other components of the Earth system at different scales by its radiative and thermal properties which modulate transfer of energy and mass at the surface-atmosphere interface (Frei et al., 2012). Snow cover can affect glaciers through changes in albedo and mass balance and also affects permafrost through its thermal insulation properties and meltwater input (Beniston et al., 2018). Additionally, snow is critical for water resources during the melt season, especially in areas lacking glaciers. Snow occurrence is also relevant for regional planning, including vehicle traffic, pasture and livestock management, its retarding effect on flash floods and slope erosion during extreme precipitation events and winter tourism and sports. In fact, snowbased tourism and sport was identified as particularly vulnerable to climate change. This climate dependency together with its economic impact on ski tourism (Steiger and Abegg, 2018), has prompted the authorities and private sector of the Magallanes Region in austral Chile to re-think the future of the current ski resort "Club Andino" located at Cerro Mirador, close to Punta Arenas city, which has barely operated for 10 days in the last two snow seasons (2016)(2017). A possible alternative is a new ski resort planned to be located at Tres Morros, 25 km southwest and 200 m higher than Cerro Mirador.
The climate along the west coast of southern South America, including southern Patagonia, is determined by the permanent high pressure anticyclonic circulation (the subtropical South Pacific High, SPH) located in the south-eastern Pacific Ocean centered around 30 • S, the circumpolar low pressure belt that surrounds the Antarctic continent around 60 • S, and the resulting westerly wind regime that prevails in mid-latitudes (40-55 • S). The seasonal north-south displacement of the SPH and the circumpolar trough modulate the westerly flow, and therefore, defines the eastward trajectory of the low-pressure frontal systems that cross southern South America. In particular, the west side of southernmost Patagonia is frequently affected by passing fronts and low-pressure systems that results in a high number of cloudy (77-86%) and precipitation days (81-88%) throughout the year (Carrasco et al., 2002). In addition, the complex topography defines a distinct climate regime at both sides of the Andes, with a relatively warmer and moister environment on the western side in comparison with the eastern side (Schneider et al., 2003;Weidemann et al., 2018). The precipitation is much higher and increases with altitude on the western slope due to orographically induced upward motion, and rapidly decreases toward the east due to the foehn mechanism (Weidemann et al., 2013).
The scarcity of rain-gauge stations precludes exhaustive regional and temporal analysis. The last IPCC report (Hartmann et al., 2013) reveals this difficulty. Their assessment indicates a negative trend of the annual precipitation ranging between 5 and 10 (2.5 -5) mm per decade around Punta Arenas for the 1901-2010 (1951-2010) period using the GHCN (Global Historical Climatology Network of NOAA) dataset on the western coast of Patagonia (see their Figure 2.29). Earlier, Carrasco et al. (2008), analyzing the 1958-2006 period, showed an overall increase of the precipitation registered at the Punta Arenas airport from around 1967 to 1991, but no trend thereafter. Later, González-Reyes et al. (2017) with data updated to 2014, indicated that large interannual and inter-decadal precipitation variability has been observed for the last 100 years at Punta Arenas. Although a long-term trend is not evident from this record, a significant precipitation reduction has been observed after 1990, especially in spring and summer, while the winter precipitation has shown a significant increase.
Furthermore, an overall temperature increase has been detected in southern South America. Early work by Ibarzábal et al. (1996) and Rosenblüth et al. (1997) suggested that the warming has been higher in the eastern side of the Andes than on its western side. Falvey and Garreaud (2009) showed a strong temperature contrast in central and northern Chile (17 • -37 • S) for the 1979-2006 period, with a coastal cooling of −0.2 • C/decade and warming of +0.25 • C/decade only 100-200 kilometers further inland. This cooling could be linked to an intensification of the South Pacific Anticyclone during recent decades, that enhances ocean upwelling off the coast of Chile. However, a recent study suggests that this recent period of coastal cooling has ended, modulated by a phase change of the Interdecadal Pacific Oscillation (Burger et al., 2018).
Therefore, there is a need for systematic studies to characterize the spatial and temporal distribution of snow cover in the Southern Hemisphere, and particularly in southernmost Patagonia, since all of the scarce long-term in-situ records from South America are from the central Andes (Vaughan et al., 2013). Here, we reduce this gap by using a remote sensing approach based on MODIS satellite data to reconstruct a snow cover time series and to relate snow conditions to present and past climate variations. This study is framed within a broader science programme that aims at contributing snow and climate information for regional planning in the Magallanes Region in southernmost Chile. The main problem we wish to target is whether snow cover has shown significant changes during the last few decades at Brunswick Peninsula, and if so, which are the main forcing factors. As such, this article will contribute to predict the future snow evolution in the area of Punta Arenas, the largest Chilean city in southern Patagonia.

STUDY AREA
The study area is located on Brunswick Peninsula in southern Patagonia between 52.9 • and 53.9 • S and 73.0 • and 70.4 • W. This peninsula is part of the southern extension of the Patagonian Andes where Punta Arenas is the largest city and capital of the Magallanes Region, the southernmost region of Chile. It is surrounded by the Strait of Magellan to the east, and by Seno Otway to the west. Most of the southwestern part is dominated by small mountains of the southern Andes, with many summits between 500 and 1,000 m, but all lower than 1,200 m a.s.l. Contrastingly, the northeastern part exhibits mostly a flat topography below 200 m a.s.l., which is part of the Patagonian steppe ("pampa"). The first mountains to the south are located west of Punta Arenas, with a maximum altitude of 590 m a.s.l. at Cerro Mirador. On the eastern slope of this summit the skiing resort "Club Andino" is located at 53.15 • S, 71.05 • W. A future skiing resort is planned at "Tres Morros" at 53.31 • S, 71.29 • W, 866 m a.s.l., located 25 km SW of Cerro Mirador. The area has a strong Pacific Ocean influence, through the prevailing westerly circulation regime. This in combination with the unique topographic distribution of the Patagonian Andes creates one of the most pronounced precipitation gradients on earth (Carrasco et al., 2002;Schneider et al., 2003;Smith and Evans, 2007;Garreaud et al., 2013;Weidemann et al., 2018). Due to the north-south orientation of the Andes almost perpendicular to the westerly airflow, the regional west-to-east variation along the 53 • S latitude and across the southern Andes shows an annual precipitation of about 6,000 mm on the western slopes of the mountain range decreasing to <1,000 mm on the eastern slope (Schneider et al., 2003) and to 416 mm at Punta Arenas (González-Reyes et al., 2017). Figure 1 shows the geographical location and the three different specific study areas where the snow-cover variability has been reconstructed. The first chosen area is the "Macro zone" of 1,637 km 2 covering a large part of the Brunswick Peninsula, used here as representative of the snow cover distribution on the entire peninsula. Two smaller embeded zones are respectively centered on the existing ski resort "Cerro Mirador" (302 km 2 ) and the projected ski resort "Tres Morros" (291 km 2 ). These two locations, during the time frame of the analysis, proved to have less cloud coverage than the southern Macro zone area, allowing adequate snow sensing by visible and near-infrared MODIS imagery.

DATA AND METHOD
In this section we describe the different data and methods used to reconstruct the inter-annual snow cover variability and the spatial and temporal evolution. Data include MODIS satellite information, ground stations, radiosonde, and MERRA-2 reanalysis. Calculations including snow reconstruction and statistical applications were made using the open source programing language Python v._3.5 (http://www.python.org). Maps (Figures 1, 2) were depicted using the open source QGIS v_2.18 1 A couple of figures (9 and 13B) were made using Matplotlib Python graphic package (Hunter, 2007).

MODIS Data and Snow Cover Detection
To reconstruct the extent of the seasonal snow cover at Brunswick Peninsula, we take advantage of the snow spectral reflectance properties given by the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data at 500 meters spatial resolution. The high temporal resolution (<2 days) of MODIS data facilitates selection of cloud-free imagery at the latitude of Brunswick Peninsula. Additionally, its multispectral capabilities make it one of the most advanced optical sensors available for large-scale terrestrial applications (Salomonson et al., 1989).
Several MODIS products are available, including the MOD10 snow cover data in its version MOD10A1 (daily) and MOD10A2 (8-day composite). Both products use the Normalized Difference Snow Index (NDSI) algorithm based on MODIS reflectance data of bands 4 and 6 to discriminate snow from other targets (Winther and Hall, 1999;Sirguey et al., 2009). This method is usually preferred instead of other classification methods such as the supervised classification (Wang and Li, 2003).
The MODIS cloud product (MOD35) is also incorporated into the algorithm to construct the snow products. Nevertheless, inaccurate detection of clouds show some problems in high elevation regions such as the Southern Alps of New Zealand (Hall et al., 2002;Sirguey et al., 2009) and the Northern Patagonia Icefield (Lopez et al., 2008;Sirguey et al., 2009). Some corrections have been implemented in updated versions of the snow cover product but snow and cloud discrimination problems persist (Hall and Riggs, 2007). Recurrent misclassifications limit the suitability of the MOD10 snow product for detecting the snow cover in frequently cloud-covered regions, such as southernmost Patagonia.
Therefore in this study we use the MODIS land product MOD09GA from the Terra satellite in its version 6 (Vermote and Wolfe, 2015) at daily temporal resolution. We did not include the MODIS land product derived from the Aqua satellite due to malfunctioning of one of the spectral bands (Saavedra et al., 2018). The product MOD09GA includes corrected reflectance with a spatial resolution of 500 m (available from the Land Processes Distributed Active Archive Center, http://lpdaac.usgs.gov/) for the period 2000-2016. To optimize the download process, we used the script "order_MWS.pl" (available from https://ladsweb.modaps.eosdis.nasa.gov/toolsand-services/), downloading a total of 6079 images. We downloaded the satellite data covering approximately an area of 31,500 km 2 , in order to obtain a general view and select image subsets.
To better discriminate potential cloud cover we conducted a visual inspection of the all images, which resulted in a selection of 484 could-free images. Sixty percent of the selected images (290) correspond to the period between April to September, temporal frame that we name "extended winter". At least, there is one cloud-free image for each extended winter month between 2000 and 2016, except for May and August 2000. Figure 2 shows an example as a comparison between true-color (RGB: spectral MODIS bands B1, B4, B3) and false-color (RGB: spectral MODIS bands B3, B6, B7) compositions to help selecting cloud-free imagery.
A water mask was generated using the coastal margin of the Regional Government of Magallanes (http://www. goremagallanes.cl) and the information of lakes from the Dirección General de Aguas, (DGA, General Water Management of Chile, http://www.dga.cl/productosyservicios/mapas/ Paginas/default.aspx#nueve) to exclude any possible interference between the snow detection and water bodies during in the analysis.

Snow Detection Indices
There are several alternative methods for snow detection using different sensors and wavelengths. In this work, we computed two different indices from MOD09GA bands 1, 4, 6 and 7 (specified below): the NDSI and the Melt Area Detection Index (MADI). Both indices plus the calculation of the new proposed Index (Snow power index), were applied to each selected image and subsequently grouped to monthly means to obtain a monthly representation of the snow cover.

Normalized Difference Snow Index
NDSI is the most common method used for detection of the snow cover (Frei and Lee, 2010). It uses the contrast between the high reflectance (ρ) of snow in the green part of the visible spectrum (MODIS band 4, 545-565 nm) and the low reflectance in the short-wave infrared part of the spectrum (MODIS band 6, 1,628-16,52 nm). We used the threshold value NDSI ≥ 0.4 (Hall et al., 2002) for identification of snow pixels based on: (1)

Melt Area Detection Index
The MADI approach uses the reflectance ratio between the red portion of the visible spectrum (MODIS spectral band 1, 0.62-0.67 µm) and the near-infrared (MODIS spectral band 7, 2.105-2.155 µm). Its principle is based on identifying liquid water layers that coat snow and ice when melt occurs to discriminate between frozen and liquid water (Chylek et al., 2007) according to: MADI has been usually used to detect wet and dry snow over regions of known snow cover only. A study carried out by Bormann et al. (2012) in Australia considered a suitable threshold for snow cover detection similar to the NDSI method. This method establishes ranges for the MADI values using in-situ observations for snow detection. If the MADI value exceeds 6, then the pixel is characterized as snow. Additionally, Bormann et al. (2012) used the MADI values to infer snow depth validated with in-situ observations. We combine the MADI and NDSI according to the following two steps procedure (see Figure 2): 1. We detect snow cover area (km 2 ) inside of each study zone.
For this purpose, the NDSI is used as a binary snow cover detector representing snow-cover extent. 2. For the area detected above, we apply the MADI to each pixel classified as snow. Using a zonal mean we obtain a zonal representation expressed as a weight indicator of snow depth, according to Bormann et al. (2012).
With these two different indices we reconstruct and compare the inter-annual variability of snow cover for the period 2000-2016. Therefore, the NDSI is used to detect snow pixels and MADI is incorporated in a newly proposed index, called Snow power , that incorporates the performance of both indices. We used the following equation for computing Snow power index: Where NDSI mask (i) is the monthly average of the snow cover extension (pixels) for month (i), MADI zonal average (i) is the zonal average MADI value for all snow-covered areas for the month (i), and max NDSI mask is the maximum snow cover monthly NDSI value for all months considered.
Thus, this index is the result of the percentage of snow cover presence (NDSI) in a given area scaled by the zonal average MADI values for all the snow cover extent, divided by the present snow cover and its maximum monthly value for all months considered. We postulate that this relationship represents better the snow cover since it incorporates both snow cover fraction, in each zone of interest, and a weight linking it to snow depth obtained by MADI so that we derive the variability of the entire snowpack. In this study, we use both the snow cover extent using only NDSI, and the Snow power index as a combination of NDSI and MADI indices.

Weather Station Data
A major difficulty in Patagonia is the scarcity of locally observed meteorological data, since weather stations are restricted to the main settlements and many stations have significant data gaps. A main weather station is Carlos Ibáñez Airport (53.00 • S, 70.53 • W, 38 m a.s.l), 18 km to the north of Punta Arenas city center. It has permanent atmospheric observations since 1958 and belongs to the Dirección Meteorológica de Chile (DMC, the National Chilean Weather Service). In 1970 the Instituto de la Patagonia deployed the weather station "Jorge Schythe" (Santana et al., 2009) at 53.07 • S and 70.52 • W at 5 m a.s.l, within the city of Punta Arenas, 14 km south of the Punta Arenas Airport. In 1973 DGA (Dirección General de Aguas) installed at Jorge Schythe station a pluviometer, a thermometer and a relative humidity sensor. Data are available for Jorge Schythe station at a monthly level since 1972, and at a daily level since 1980. Another long-term meteorological record (more than 100 years) is also available from the Navy station Evangelistas Lighthouse (52.24 • S, 75.06 • W, 55 m a.s.l.) located at the western entrance of the Strait of Magellan.
The available station data allow studying the inter-annual and decadal pattern of air temperature and precipitation in relation to forcing mechanisms such as the El Niño Southern Oscillation (ENSO), the Antarctic Oscillation (AAO) or the Southern Annular Mode (SAM), and the Pacific Decadal Oscillation (PDO). Other stations have been deployed in recent decades in southern South America by DGA (e.g., Teniente Merino, San Gregorio, Pampa Huanaco and Bahía San Felipe), by DMC (Porvenir station) and by specific research projects such as Bahamondes station (Schneider et al., 2003;Weidemann et al., 2018). Table 1 lists the above stations used in this study.

Radiosonde Data Punta Arenas
At Punta Arenas airport, radiosonde balloons are launched every day at 12 UTC (Universal Time Coordinated) since April 1975. This upper-air station measures a daily vertical profile of air temperature, relative humidity, wind direction and wind speed at different standard levels of the atmosphere. The zero-isotherm altitude (ZIA) can be determined from the air temperature profile. Monthly radiosonde data can be obtained from the Anuario Climatológico (Climatological Yearbook) available online from www.meteochile.gob.cl. Here, monthly averages of surface air temperature and 850 hPa (hectopascal) levels were used to determine the monthly altitudinal lapse rate which along with the geopotential height at the 850 hPa allowed to calculate the monthly and annual ZIA.

Reanalysis Data
The monthly snow cover product from MERRA-2 reanalysis data ("fractional snow-covered area FRSNO") for April-September (Gelaro et al., 2017) period were used for comparison with the MODIS snow cover reconstruction from 1980 to 2016. MERRA-2 data were downloaded from https://giovanni.gsfc.nasa.gov. MERRA-2 has a resolution of 0.5 • × 0.625 • latitude-longitud grid with 72 model layers interpolated to 42 standard pressure levels (Bosilovich et al., 2015). It uses the modified three-layer snow model of Lynch- Stieglitz (1994). This snow hydrology model resolves for snow melting and refreezing, as well as dynamic changes in snow density, snow insulating properties, and other physics relevant to the growth and ablation of the snowpack (Stieglitz et al., 2001). The snow model component of the Land Surface Model (LSM) uses the minimum snow water equivalent (WEMIN) before it considers the surface to be snow covered (Reichle, 2012). In this sense, if the amount of snow (SNOMAS) on a given model-surface area is insufficient to cover the entire area with at least an amount of WEMIN, the snow cover fraction (FRSNO) will be the fraction of the area that is covered by an amount equal to WEMIN, that is expressed as, FRSNO = min(1, SNOMAS/WEMIN) (Reichle, 2012;Reichle et al., 2017). No studies have so far been conducted on how well the MERRA-2 can simulate snow cover in Patagonia.
The MERRA-2 snow cover values were extracted for the grid point 53.3 • S, 71.1 • W within the Macro zone. The values are expressed as a monthly snow cover fraction, with a value of 1 representing 100% of snow cover. For example, the maximum monthly FRSNO value for the period 2000-2016 was 0.38 for June 2002, when the NDSI snow cover extent indicated an area of 1,443 km 2 for the Macro zone, equivalent to 88% of the total domain area of 1,637 km 2 .

Climate Indices
ENSO is the most significant mode of inter-annual climate variability in the Southern Hemisphere with world-wide impacts. El Niño (La Niña) corresponds to the warm (cold) phase which is characterized by weaker (stronger) trade winds and above (below) normal sea surface temperatures (SST) in the equatorial Pacific Ocean. Also, a weaker (stronger) SPH is encountered during El Niño (La Niña) with stronger (weaker) subtropical and weaker (stronger) polar jet-streams (Chen et al., 1996). In southernmost South America, a weak correlation between ENSO and precipitation occurs to the west of the Patagonian Andes between 45 • and 55 • S (Schneider and Gies, 2004), i.e., with below (above) normal precipitation during El Niño (La Niña) (Montecinos and Aceituno, 2003). The PDO, also called Interdecadal Pacific Oscillation (IPO) has a cycle of 15-30 years that shows similar SST and sea-level pressure anomaly patterns as ENSO (Garreaud and Battisti, 1999;Mantua and Hare, 2002). During the warm (or positive) PDO phase, El Niño events prevail in frequency and intensity. On the contrary, during the cool (or negative) phase, the opposite pattern occurs and La Niña episodes are more frequent and last longer (Andreoli and Kayano, 2005). In southern Patagonia the prevailing large-scale pattern associated with a positive PDO is an anticyclonic circulation centered in the southeastern Pacific Ocean to the west of the Antarctic Peninsula (Garreaud and Battisti, 1999). This pattern might reduce precipitation along the southern tip of South America, so that the storm systems may be steered northwards (Garreaud and Battisti, 1999).
Another natural internal mechanism for interannual variability is related to the SAM (also known as the Antarctic Annular Oscillation AAO) a large-scale alternation of atmospheric pressure between the mid-latitudes and high latitudes at the surface pressure level (Gong and Wang, 1999). These changes are linked with changes in the strength of the Southern Hemisphere westerly winds due to the latitudinal shift in the area of maximum zonal flow (∼45 • S) or in its intensity (surface wind speed) (Moreno et al., 2010). Positive correlations between the westerlies and precipitation prevail along the western margin of the southern South America between 40 • and 55 • S, while correlations tend to be negative at similar latitudes along the Atlantic margin (Moreno et al., 2010). Thus positive (negative) AAO means southward (northward) displacement of the westerlies, and therefore, of the storm track in the Pacific Ocean. This can also favor (inhibit) precipitation events in the western coast of southern South America, although it depends how far the southward (northward) displacement extends and whether the area of interest is governed by orographic rain or rather by foehn effects on the east side of the Cordillera.
These different mechanisms are expressed by their respective indices. The interannual variability of ENSO can be studied by using SOI, a normalized sea-level pressure difference between Tahiti and Darwin stations (Chen, 1982;Trenberth, 1997). Alternatively, among other indices, the SST monthly anomalies in the central-equatorial sector of the Pacific Ocean, known as Niño 3.4 can be employed as an ENSO index (Trenberth, 1997). The PDO index is calculated by the leading principal component applied to the spatial monthly averaged SST of the North Pacific Ocean north of 20 • N, using mainly the October to March period so that the interannual fluctuation mostly takes place in winter (Mantua and Hare, 2002). The AAO is calculated by principal component analysis so that the first three modes explain 33, 11 and 9% of the variance, respectively, and it can be constructed using the 700 hPa height anomalies poleward of 20 • S (Gong and Wang, 1999). The above indices were downloaded from the Earth System Research Laboratory www.esrl.noaa.gov/ psd/data/climateindices/list/, specifically, the SOI from www.cpc. ncep.noaa.gov/data/indices/soi, the PDO from www.esrl.noaa. gov/psd/ data/correlation/pdo.data and the AAO from www.esrl. noaa.gov/psd/data/correlation/aao.data.

Statistical Analysis
In order to analyse the variability of snow cover extent, statistical correlations were performed with different climate variables. We use the Pearson correlation quantified by the coefficient of determination (R 2 ), the Durbin-Watson statistics for detecting autocorrelation between variables, and significance tests with t-statistic using a p-value at 5% significance. The structure of the residuals is also analyzed to detect any biases. Based on the best relationship we express the satellite-derived snow cover as a function of climate variables, which allows to reconstruct the snow cover to 1965. In order to detect long-term trends we use an exponential filter (Rosenblüth et al., 1997) to reduce the interannual variability and apply a linear regression. Python Statsmodels statistic package was used (Seabold and Perktold, 2010).

Precipitation and Temperature From Surface Stations
The average monthly precipitation 2000-2016 for 8 weather stations aligned from west to east around 53 • S is shown in Figure 3. Precipitation is much higher (around one order of magnitude) in the western stations Evangelistas and Bahamondes located in the windward side of the Andes than the stations located in the leeward side (for example Punta Arenas).
Precipitation occurs year-round with only little seasonality, although the minimum (maximum) tends to occur at the end of winter (summer)-beginning of spring (autumn). In general, precipitation increases from north to south although Bahía San Felipe station shows higher accumulation than Pampa Huanaco located further south. The annual variability during the 2000-2016 period is displayed in Figure 4. There are no correlations among different stations except between the two nearby ones located in Punta Arenas (airport and city). This reveals the large spatial variability of precipitation in southernmost South America. Also, no common and significant trends can be observed in the analyzed stations. Thus, while Evangelistas suggests an increase, Punta Arenas city and Bahía San Felipe suggest instead a decrease in precipitation in the last 17 years. This finding was already mentioned by González-Reyes et al. (2017) who analyzed the long-term behavior for the 1900-2014 period using precipitation data from Observatorio Meteorológico Monseñor Fagnano located in Punta Arenas city center (where the Salesian mission began measurements in 1887) and Punta Arenas Airport where observations started in 1930. The Airport station was first located in Bahía Catalina until 1958, 4 km north of the city center, and then moved 14 km north to Punta Arenas Airport (Presidente Carlos Ibáñez del Campo) where DMC continues the measurements at present. González-Reyes et al. (2017) found a significant reduction in the annual precipitation sums after 1990, especially in spring and summer, while the winter precipitation showed a significant increase during that period.
The other long-term precipitation record is from Evangelistas with data available from 1899 onward, operated by the Chilean Navy. Monthly data and some complete years are missing which preclude exhaustive analysis. Nevertheless, 80% of the data are available and only 8 years have more than 6 months missing, with 3 complete years missing. Figure 5 shows the annual and extended winter (April-September) annual variability and longterm trend of precipitation. It was constructed using the Evangelistas data, the missing monthly averages were replaced by the respective long-term monthly averages to facilitate data management. The vertical lines on the annual data are the sum of the 1st standard deviation of the missing months which assigns an upper and lower boundary for the most probable annual precipitation. The solid curve results by applying an exponential filter to the data to eliminate the inter-annual variability and highlighting the interdecadal variability (Essenwanger, 1986;Rosenblüth et al., 1997). This filter is similar to a moving mean, but it has the advantage of keeping the complete series for analyzing relatively short periods.
The annual data show an overall decrease in precipitation from 1899 until 1980, after which a significant increase took place between 1984 and 1998. The precipitation decline observed by González-Reyes et al. (2017) during 1990 at Punta Arenas is not completely captured by the Evangelistas data, although an overall decrease could have taken place during 1997-2007. The significant precipitation increase from around 1984-1998 has been a matter of discussion of whether it is a real signal or rather a result of an instrument change or relocation of instrument or change in staff and operation mode (Carrasco et al., 2002;Aravena and Luckman, 2009). According to the Servicio Meteorológico de la Armada (Navy Weather Service) the weather station has not been changed in location, nor the instrument type or operation mode, only the staff is replaced every one or 2 years (personal communication Ismael Escobar, Chilean Navy). The insert in Figure 5 shows an enlargement of the 1975-2000 period, suggesting a rapid but gradual increase of precipitation. However, the magnitude of the increment (of about 70% between 1984 and 1998) should be taken with caution, even though Evangelistas is located just to the south of the maximum precipitation region (Miller, 1976) and a slight southward displacement of the westerlies could indeed trigger such profound change in overall precipitation. Figure 6A shows the monthly average of the mean air temperature at Evangelistas, Bahamondes and Punta Arenas Airport over the 2000-2016 period. The ocean and continental influences are clearly seen by comparing the annual cycle at each station located in the latitudinal transect around 53 • S. The average winter (summer) air temperature at Punta Arenas Airport is about 3 (2) • C colder (warmer) than at Evangelistas. The inter-annual variability and the long-term trend (solid lines) of the extended winter (April-September) mean temperatures for Evangelistas and Punta Arenas are shown in Figure 6B for the available 1965-2016 and 1958-2016 periods, respectively. The temperature trends within both periods show a significant warming in all seasons (especially in the winter) including the montly mean and minimum temperatures. A total warming of 0.54 • C (0.11 • C/decade) and 0.72 • C (0.12 • C/decade) is found at Evangelistas and Punta Arenas Airport over the 1965-2016 and 1958-2016 periods, respectively, computed using the interannual extended winter data (smoothed data).

Upper Air Sounding
The analysis of upper air data from Punta Arenas Airport is restricted for the lower half troposphere given by the standard levels at 850, 700, and 500 hPa which respectively correspond to the annual (January-December) mean altitudes of 1,338, 2,856, and 5,377 m. This corresponds to an annual mean air temperature of −1.5, −9.7, and −26.1 • C, respectively. The chosen levels show positive linear trends of +0.01, +0.25, and +0.45 • C (non-statistically significant) over the 1976-2016 period at 850, 700, and 500 hPa, respectively, revealing a moderate warming trend in the mid-troposphere. Also, the zonal component of the wind (u-wind) at surface, 850, 700, and 500 hPa shows positive linear trends of 0.81, 1.73, 2.39, and 2.45 ms −1 for the 1976-2016 period, respectively (Figure 7). The increase in wind speeds is not statistically significant. However, this finding might be indicative of southernmost South America   becoming more and more affected by a windier environment due possibly to the enhancement and southward displacement of the westerly winds associated with positive SAM during recent decades (Marshall, 2003). Following the analysis of Schneider and Gies (2004) and Garreaud et al. (2013), the correlation between precipitation and u-wind at 850 hPa reveals that with stronger zonal wind there is an increase in precipitation westward from the Patagonian Andes and a decrease over the eastern side. The overall different behavior of the precipitation observed at Evangelistas and Punta Arenas may be the result of the overall increase in the u-wind component.
The mean annual freezing level calculated from the radiosonde data shows a 49 m non-statistically significant increase during the 1976-2016 period. The highest freezing level recorded in the period occurred in 2016 with an altitude of 1,304 m, 263 m higher than the average in the whole period (see Figure 7). The year 2016 also showed the lowest annual precipitation since 1970 according to the Punta Arenas Airport station record.

Snow Cover Reconstruction and Variability
In Figure 8, we show the monthly average of the snow cover area using the NDSI and the MADI indices for each domain namely the Macro zone, Cerro Mirador, and Tres Morros (Figure 1). We assume that the beginning and the end of the snow season occur when the monthly snow cover extent rises or drops below 20% of the multi-year maximum montly snow extent, that is 327 km 2 for the Macro zone, 58 km 2 for Cerro Mirador and 60 km 2 for Tres Morros.
With this approach, we can identify the annual duration of the snow season which usually begins in May and ends in September for Cerro Mirador and in October for Tres Morros. The relation between the snow cover area and MADI values is shown in Figure 8. The monthly averaged snow cover duration is 4 weeks longer at Tres Morros compared to Cerro Mirador due to higher altitudes. Similar Snow power index values can be observed between Macro and Cerro Mirador zones, while Tres Morros shows higher values due to its higher altitude. Figure 9 shows scatter plots between the Snow power index and the snow cover extent derived from the NDSI and MADI indices. Coefficient of determination (R 2 -values) higher than 0.86 are obtained for the Snow power index between each pair of zones especially between the Macro zone and Cerro Mirador (R 2 = 0.96, Figure 9A). This is reasonable considering that the Macro zone includes the other two zones and considers a more varied topography. A slightly smaller Snow power index R 2 = 0.86 is obtained between Cerro Mirador and Tres Morros (Figure 9C), which can be attributed to its higher altitude which results in a larger and more permanent snow cover extent, unlike to the case of Cerro Mirador.
At Cerro Mirador the ski club "Club Andino" has recorded information of the ski season duration based on the opening and closing dates of the ski lifts and the snow height average during the ski season from a stake located close to the main hut. There is no artificial snow production at the ski club so the start and end of the ski season depends largely on snow precipitation and air temperature, with a minimum natural snow height requirement of about 20 cm without compaction to start the skiing operation. Since 2010 the longest season occurred in 2015 with 95 snow-days (representing the duration of the skiing season) and an average snow height of 140 cm. Figure 10A displays the relation between the duration of the ski season for the Macro zone and Cerro Mirador and its reconstructed snow cover extent from the NDSI index. We use the June-August period for the ski season, when most of the skiing activity is concentrated. The linear correlation between the ski season duration (days) and the Cerro Mirador snow cover ( Figure 10B) is R 2 = 0.94, and R 2 = 0.87 for Macro Zone (Figure 10C). The correlation of the ski season duration with the Snow power index is also good but somewhat smaller, with R 2 = 0.77 for Cerro Mirador and R 2 = 0.83 for the Macro Zone, respectively. The correlation between of the average snow height from Club Andino with the NDSI snow cover extent (Snow power index) is not as good as the correlation with the ski season duration, with R 2 = 0.72 (0.80) and R 2 = 0.75 (0.67) for the Macro Zone and Cerro Mirador, respectively. Based on these results, hereafter we consider the NDSI snow cover extent as a better indicator of snow cover variability than the Snow power index. Figure 11 shows the climate indices and snow cover of the Macro zone at Brunswick Peninsula after smoothing the interannual variability using the exponential filter (Rosenblüth et al., 1997) over the 2000-2016 period. It reveals the overall positive trend of the AAO, i.e., the southward displacement of the westerly winds, which could be linked with the decrease of   cover behavior, although strong El Niño (La Niña) events tend to concur with less (more) snow cover and precipitation in the study area. This could be due to the interannual AAO that can also suppress or favor precipitation in the southern tip of South America. Figure 12 shows the number of precipitation days with air temperature <2 • C during April and September (extended winter) for each year for the 1958-2016 period at Punta Arenas Airport, along with the average of air temperature for these days. Also displayed is the air temperature average during days following a snow event. Results indicate a warming of 0.65 • C (p < 0.05) for precipitation days over the 1958-2016 period. On average, 11.3 ± 4.4 (±1 standard deviation) of precipitation days with air temperature <2 • C occur during the extended winter in the period 1958-2016. These days are regarded as "snow days", revealing a decline of 11% (1.2 "snow days, " p < 0.05) over the 1958-2016 period. The long-term behavior obtained after applying the exponential filter to the interannual variability shows an increase from 1958 to 1971 and from 1982 to 1999, as well as a decrease from 1972 to 1981 and also from 2000 to 2016. During 2000-2016 the linear trend indicates a decrease of 18% (2 "snow days, " p < 0.05). The mean air temperature for days following the snow events reveals a slight positive trend of 0.27 • C (p < 0.05) over the 1958-2016 period, of 0.52 • C (p < 0.05) over the 1973-2016 period and of 0.55 • C (p < 0.05) over the 2000-2016 period. In summary, this analysis indicates an overall decrease of snow precipitation days along with an overall increase of the air temperature, mostly concentrated in recent decades (Figure 12). Figure 12 includes as well the snow fraction obtained from the MERRA-2 reanalysis data. Correlations of R 2 = 0.16 (p < 0.05) are obtained between the snow fraction and air temperature after snow events; R 2 = 0.41 (p < 0.05) between the snow fraction and air temperature during precipitation days; and R 2 = 0.42 (p < 0.05) between snow fraction and air temperature during "snow days."

DISCUSSION
The snow cover reconstruction can be analyzed in relation to climate variables at a montly and annual level, such as freezing level (upper air sounding at Punta Arenas Airport), snow accumulation (measured daily at Jorge Schythe and modeled at Punta Arenas Airport), and minimum and average temperatures TABLE 2 | Result of the different linear regressions between Snow power index and NDSI snow cover (km 2 ) for the Macro zone, with the Punta Arenas Airport freezing level from radiosonde data, precipitation (mm), minimum, and mean temperature (monthly averages) from Punta Arenas Aiport and Jorge Schythe stations, for two different periods in 2000-2016: April-September (monthly and annual values) and June-August (annual).

April-September
April-September June-August Precipitation Airport ( Frontiers in Earth Science | www.frontiersin.org (at Jorge Schythe and Punta Arenas Airport). These variables were chosen because of the better correlation with the snow cover (summary shown in Table 2). The measured daily snow accumulation data are limited to the Jorge Schythe station only, published here for the first time. We use a threshold relation with air temperature to derive the solid/liquid precipitation fraction to infer snow accumulation [c(i)] at Punta Arenas Airport, as indicated by Schaefer et al. (2013), as follows: Where q(i) is the solid part of the total precipitation, prescribed by temperature as: Where T(i) is the mean air temperature ( • C), calculated as the average between the daily maximum and minimum air temperatures, P(i) is the daily total (liquid and solid) precipitation accumulation (mm), Ut and Lt are upper and lower temperature thresholds ( • C).
Ut and Lt were estimated as 1.6 and 0.9 • C, respectively. These values were determined through a sensitivity analysis using the daily snow accumulation and mean air temperature at Jorge Schythe for 1980-2016. In our case Ut and Lt are slightly larger than the values of 1.5 and 0.5 • C assumed by Schaefer et al. (2013) for the Northern Patagonia Icefield. The modeled precipitated snow has a correlation of R 2 = 0.48 with the measured snow values at a daily level, R 2 = 0.76 for monthly values, and R 2 = 0.88 for annual values. Equation (5), based on daily values, was applied at Punta Arenas Airport to infer its monthly snow accumulation.
In this analysis we only used the Macro zone indices due to the high level of correlation with the other two study zones (Figure 9). This allows to extend back the period of analysis earlier than the year 2000 in order to identify long-term trends of snow cover behavior. Based on these results, to find a better relationship for explaining snow cover variability, we tested both simple and multiple linear regressions.
The multiple regressions with the best significant correlation for the MODIS snow cover reconstructed for the Macro zone in the period 2000-2016 was found using the monthly snow accumulation and the mean monthly temperature, with R 2 = 0.62 (Equation 6, which allows to model the snow cover for the period 1972-2016), and R 2 = 0.62 (Equation 7, which allows to model the snow cover for 1958-2016 period), both statistically significant at 95% (Figure 13). Temperatures were chosen at the Airport station to avoid possible urban effects at Jorge Schythe station located within the city of Punta Arenas. Correlations for the "extended winter" period from April to September were better than the annual (January-December) and the winter (June-July-August) correlations.
Inspection of Table 2 shows that the snow cover extent reconstructed from the NDSI index derived from MODIS satellite data shows better correlations with temperatures (minimum and average) and freezing level than the snow cover derived from the Snow power index at monthly level. However, the Snow power index shows a slightly better correlation with total precipitation (at Jorge Schythe and Punta Arenas Airport stations) and snow accumulation (at Punta Arenas Airport), at seasonal level, compared with the NDSI index. This suggests that temperature may have a greater impact in the snow cover extent, while the winter precipitation may cause a more direct influence on the snow depth as reflected in the Snow power index. In-situ measurements are needed to verify these relations.
• Modeled Snow cover (SNC 1 ) in km 2 for 1972-2016: • Modeled Snow cover (SNC 2 ) in km 2 for 1958-2016: Where T • Apt ( • C) is the monthly mean air temperature measured at Punta Arenas Airport during April to September. SN Sch (cm) is the monthly snow accumulation at Jorge Schythe during April to September. SN Apt (cm) is the monthly snow accumulation modeled at Punta Arenas Airport during April to September.
When we correlate the NDSI snow cover extent for the Macro zone for 2000-2016 with the monthly snow accumulation at Jorge Schythe station a correlation of R 2 = 0.26 is obtained, which increases to R 2 = 0.59 when correlated with the monthly mean temperature at Punta Arenas Airport. This means that the snow precipitation alone accounts for 26% of the variance of the NDSI snow cover extent, while the air temperature by itself accounts for 59% of the variance of the NDSI snow cover extent. Therefore, temperature is the main driver of snow cover changes, while the snow precipitation has a less relevant influence. The air temperature acts both as a forcing factor for solid/liquid precipitation, and also for melting the snow cover after precipitation events. As mentioned earlier, a clear warming trend of 0.52 • C can be detected in the air temperature immediately after snow events at Punta Arenas Airport station in the period 1973-2016 ( Figure 12A). The modeled snow cover (Modeled Snow cover-1, Equation 6), shows a decreasing trend for the Macro zone of 19% between 1972 and 2016, from 4,162 to 3,360 km 2 , with a 95% significance (Figure 13). This long-term snow decrease is attributed mainly to the significant long-term warming of 0.71 • C at Punta Arenas for the same period (Figures 6, 12).
When reconstructing the snow cover with the Modeled Snow cover-2 (Equation 7) for the period 1958-2016, the curve for the period 1972-2016 (Figure 13) turns out to be very similar to the Modeled Snow cover-1 (Equation 6). Expanding the period to 1958 allows seeing a slope change at the beginning of the 1970's, showing a snow cover increase from 1958 to 1973. Interestingly, this snow extent increase coincides with a Southern Hemisphere cooling of 0.2 • C (annual average temperature) as reported by Bindoff et al. (2013) and Harris et al. (2014) for the same period.  (6) (snow cover as a function of monthly mean temperature at Punta Arenas Airport and snow accumulation at Schythe station) is shown in light blue for the 1972-2016 period, while the red dashed line represents its linear trend. The annual (April-September sums) of snow cover modeled with Equation (7) (snow cover as a function of monthly mean temperature and snow accumulation modeled at Punta Arenas Airport) is shown in orange for the 1958-2016 period. The trend of the modeled snow cover is −19% for the period 1972-2016 (p < 0.05). Lower section. The annual (April-September) snow accumulation at Jorge Schythe station is shown as a purple curve. The inset (B) shows a scatter plot between the monthly values of modeled snow cover area on the Macro zone and the reconstructed snow cover from MODIS data for the 2000-2016 period, with R 2 = 0.62 and p < 0.05. This is also evident at Punta Arenas Airport, which shows a cooling of 0.4 • C for 1958-1973 (extended winter). Santana et al. (2009) reported that the annual average of the air temperature at Punta Arenas city decreased around 1 • C from 1932 to 1971. This negative trend is part of the global cooling that took place from mid-40s to mid-70s (Jones and Moberg, 2002).
Furthermore, the Modeled Snow cover-1 for 1972-2016 is compared with the monthly snow fraction obtained from MERRA-2 reanalysis data for the period 1980-2016, resulting in a significant correlation of R 2 = 0.48 and p < 0.05. In spite of the significant correlation between the modeled snow cover and the MERRA-2 reanalisys data, the MERRA-2 data show no trend within the 36-year period ( Figure 12B). In order to check the temperature trends within the MERRA-2 reanalysis dataset, the mean temperatures of the extended winter period (April-September) were analyzed for the same grid point 53.3 • S, 71.1 • W within the Macro zone at the near-surface and 850-hPa levels. Figure 6B displays the observed and the MERRA-2 simulated air temperatures for both near-surface and 850-hPa levels. It shows the overall good agreement in the interannual variability, although the MERRA-2 data show a cold bias. The long-term behaviors given by the exponential filter show positive trends for both the near-surface and 850-hPa and for the observed and MERRA-2 data during the 1980-2016 period. All curves reveal a slight cooling during the 1980-1993 period and a slight warming during the 1994-2016 period. The results show a significant correlation of R 2 = 0.83 (p < 0.05) between the observed near-surface air temperature at Punta Arenas Airport and the temperature obtained at the nearest gridpoint of MERRA-2 reanalysis data for the 1980-2016 period. Also a significant correlation of R 2 = 0.71 (p < 0.05) is found between the observed 850-hPa monthly mean air temperature obtained from the radiosonde data (12 UTC) from Punta Arenas Airport and the 850-hPa MERRA-2 data for the same period ( Figure 6C). However, there are some differences between the observed and reanalysis data. The observed near-surface warming is 0.60 • C (0.52 • C), larger than the 0.13 • C (0.17 • C) MERRA-2 warming for the1980-2016 (1994-2016) period, while the MERRA-2 warming of 0.40 • C at 850-hPa level is larger than the 0.05 • C given by the radiosonde data. Although the MERRA-2 data show good general correlations with the surface and radisonde data, its use for trend analysis is limited, as shown above. The main reason being that only few surface and upper air data are available in Patagonia for computing the reanalysis data.
Interannual precipitation changes affect the short-term snow cover variability as occurred in 2016, the year with the smallest winter precipitation for the last 50 years at Punta Arenas Airport and also at Punta Arenas city (90 mm and 165 mm, respectively, during April-September). The negative SOI phase that prevailed between 2008 and 2015 allowed the predominance of La Niña events. Therefore, the large-scale atmospheric circulation favored an increase in precipitation in the region (Figure 5). This period ended with the development of a strong El Niño 2015-2016. The year 2016 presented some extreme climate events in Patagonia and the southern ocean which has been explained through climatecouplings. Garreaud (2018) reported for 2016 a significant precipitation deficit larger than 50% in western Patagonia while at the same time the lowest Southern Hemisphere (SH) sea ice extent of the satellite record occurred around Antarctica (Stuecker et al., 2017). Both studies agree that the coupling between the ENSO and AAO was the main causative factor of the 2016 event. This shows the sensitive response of the local and regional climate to wider-scale climate conditions.
The AAO index has been increasing in recent decades (Marshall, 2003), which points to windier conditions in southern South America, this as a result of the southward displacement of the westerly winds and the associated trajectory of the frontal systems that are the main contributor of precipitation in the region. The behavior of the AAO coincides with a nonsignificant wind increase at Punta Arenas Airport, and with the precipitation increase at Evangelistas. The precipitation increase at Evangelistas is probably an evidence of the increase of the intensity of the westerlies, while the precipitation decrease at Punta Arenas Airport, particularly in spring and summer, could be an evidence of the enhanced leeward effect of the wind increase which results in drier conditions (Schneider and Gies, 2004;Moreno et al., 2010;Garreaud et al., 2013).
A recent study of snow cover variability in Aysén basin, northern Patagonia, using MODIS data for period 2000-2016 period (Pérez et al., 2018), shows a decreasing but non-significant trend of 20 km 2 year −1 in snow covered area. Similar trends have been detected in the Australian Alps where the maximum snow depth has declined by up to 15% since the 1960s, particularly in spring (McGowan et al., 2018). Model predictions suggest this trend will continue along with a reduction in length of the snow season in response to global warming (Hennessy et al., 2008;McGowan et al., 2018). In the European Alps, longterm observations show negative trends in snow depth and snow duration over the past decades (Beniston et al., 2018). The changes are typically elevation dependent, with more (less) pronounced changes at low (high) elevations due to changes in snowline altitude driven by warmer temperature. In Canada, the Canadian Sea Ice and Snow Evolution (CanSISE) Network has made an assessment of snow cover trends (Mudryk et al., 2018). The results show a decrease in terrestrial snow cover fraction in fall (delayed snow cover onset) and spring (earlier snow melt), generating a shortening of the snow season. These changes are driven by warmer temperatures and an associated modification of the elevation-dependent rainfall vs. snowfall ratios (Mudryk et al., 2018). These studies, as well as those published for other mountain regions of the world (Vaughan et al., 2013), agree that temperature increase is the main variable that drives the snow cover reduction.
By the end of this century, our study area in southern South America is predicted to warm by about 1 • C (Representative Concentration Pathway RCP2.6) to 2.5 • C (RCP8.5) (Magrin et al., 2014). An increase in westerly winds due intensification of the AAO and therefore of precipitation is expected to occur in western Patagonia in the future. However, the warming trend should result in a decrease of snow precipitation and therefore snow cover reduction at low elevation sectors of the Patagonian Andes (less than ca. 1,500 m a.s.l.) such as Brunswick Peninsula.

CONCLUSION
In this paper, we assess the snow cover changes for Brunswick Peninsula in southern Patagonia (53 • S), where the Andes attain maximum elevations of about 1,100 m. The NDSI index applied to the MODIS satellite data at 500 m resolution provides an adequate estimation of the snow cover extent for the period 2000-2016. The MODIS snow cover data were compared with MERRA-2 reanalysis data and the snow data from the Andino ski center close to Punta Arenas city. The MODIS-derived snow cover extent compares well with the snow duration and snow height recorded at Club Andino in 2010-2016, which helps to validate the satellite data.
Two significant correlations between the MODIS NDSI snow cover extent and climate variables are found for the 17-year period 2000-2016. The first is a linear relation between the monthly mean temperature at Punta Arenas Airport station and monthly snow accumulation at Jorge Schythe station (Modeled snow cover-1, Equation 6), and the second is a similar correlation between the monthly mean temperature and monthly modeled snow accumulation at Punta Arenas Airport (Modeled snow cover-2, Equation 7). The snow cover series is extended with these linear models to 1972 and 1958, respectively. The first regression (Equation 6) shows a significant 19% snow extent reduction in the 45-year period 1972-2016 of 4.2% per decade, while the second regression (Equation 7) shows two different slopes, a significant 32% snow extent increase for the 1958-1973 (20% per decade) and a significant 16% snow extent reduction for the 1973-2016 period (3.6% per decade). The slope change in 1973 coincides with a cooling observed before mid-70's in Patagonia and in Southern Hemisphere temperatures. According to Equations. (6) and (7), the 1972-2016 snow cover reduction can be attributed mainly to a long-term atmospheric warming during the winter season.

AUTHOR CONTRIBUTIONS
RJ, FA, GC, TS, and CS jointly designed the study. FA and JC processed and analyzed the climatological data. FA, RJ, and GC processed and analysed MODIS satellite data. NB provided and analyzed the Jorge Schythe station meteorological data. EG, RA, and KG compiled snow depth and meteorological data. TS and CS among others designed, installed, and maintain the measurement network at Bahamondes station. All authors discussed the results and jointly worked on the manuscript.

ACKNOWLEDGMENTS
This study has been funded by the FNDR GORE Magallanes 2048 Programme Transferencia Científico Tecnológica Modelamiento Climático Planificación, XII Región, Código BIP N o 30462410, and by the Chilean Antarctic Institute MoCliM Programme. Partial funding has been contributed by the CONICYT-BMBF project GABY-VASA grant N • BMBF20140052 (Chile) and N o 01DN15007 (Germany); by the German Research Foundation project Gran Campo Nevado grant N• SCHN 680/1-1; by the German Research Foundation project MANAU grant No SA 2339/3-1; and Conicyt-Anillo ACT-1410. Data from Evangelistas station were kindly provided by Servicio Meteorológico de la Armada de Chile (SERVIMET). Data from Punta Arenas airport and Porvenir were downloaded from the institutional webpage of the Dirección Meteorológica de Chile, while data from San Gregorio, Pampa Huanaco, Bahía San Felipe, and Teniente Merino are from Dirección General de Aguas. The authors gratefully acknowledge the valuable comments and suggestions from both reviewers and the editor that helped improving the manuscript.