Decrease in Snow Cover over the Aysén River Catchment in Patagonia, Chile

The evidence for global warming can be seen in various forms, such as glacier shrinkage, sea ice retreat, sea level rise and air temperature increases. The magnitude of these changes tends to be critical over pristine and extreme biomes. Chilean Patagonia is one of the most pristine and uninhabited regions in the world, home to some of the most important freshwater reservoirs as well as to evergreen forest, lakes and fiords. Furthermore, this region presents a sparse and weak network of ground stations which must be complemented with satellite information to determine trends on biophysical parameters. The main objective of this work is to present the first assessment on snow cover over the Aysén basin in Patagonia-Chile by using Moderate Resolution Imaging Spectroradiometer (MODIS) data from the period 2000–2016. The MOD10A2 product was processed at 500 × 500 m spatial resolution. The time-series analysis consisted in the application of non-parametric tests such as the Mann–Kendall test and Sen’s slope for annual and seasonal mean of snow covered area (SCA). Data from ground meteorological network and river discharges were also included in this work to show the trends in air temperature, precipitation and stream flow during the last decades. Results indicate that snow cover shows a decreasing non-significant trend in annual mean SCA with a −20.01 km2·year−1 slope, and neither seasonal mean shows statistical significance. The comparison with in situ data shows a seasonal decrease in stream flows and precipitation during summer. The hydrological year 2016 was the year with the most negative standardized joint anomalies in the period. However, the lack of in situ snow-monitoring stations in addition to the persistence of cloud cover over the basin can impact trends, creating some uncertainties in the data. Finally, this work provides an initial analysis of the possible impacts of global warming as seen by snow cover in Chilean Patagonia.


Introduction
The warming of the climate system and radiative forcing observed since the 1950s have increased the temperature of the atmosphere and oceans, reduced the amount of snow and ice cover, and raised the sea level (3.2 mm·year −1 ) [1]. Snow has decreased in most regions, especially during the spring and summer seasons, due to a positive feedback with the air temperature trend [2]. Over the last three decades in the northern hemisphere, the duration of snow cover has remained stable in North America, while in Eurasia it has decreased (−2.6 ± 5.6 days·decade −1 ) [3]. In the southern hemisphere, however, the general lack of available records makes it impossible to establish a trend [2]. There is evidence, however, of retreating glaciers in the Andes mountain range [4][5][6][7], although no detailed analysis has yet been carried out on the variability in snow cover, which plays a crucial role in water supply for agriculture and hydroelectric energy production, mainly in the southern area of the Andes [8,9].
Monitoring snow and its properties can be very problematic since conventional methods of measuring at high altitude are limited by rugged topography and adverse weather conditions. For this reason, satellite images are useful because they provide information on snow cover in mountainous areas at regular intervals of time [10]. Due to the physical properties of snow, it has a high reflectance in the visible spectrum (0.5-0.7 µm) and high absorption in the shortwave infrared (1.0-3.5 µm), which provide the basis of the Normalized Difference Snow Index (NDSI) for differentiating snow from other types of cover. At the same time, this reflectivity depends on factors such as the snow grain size and shape, the content of liquid water, the depth of the snowpack and the presence of impurities [10,11].
Some of the main snow products that have been derived from remote sensors include snow extent maps and snow albedo, snow grain depth and size, and the snow water equivalent (SWE) [12,13]. At global level, the northern hemisphere has the largest number of available studies and satellite products since it has continuously monitored snow cover since 1960, while in the southern hemisphere, this has only been done on specific regions and in the last decade (excluding Antarctica) [14].
The largest part of the Andes mountain range is in Chile, where it stretches 8000 km and has an average height of 4000 m, nearly all of which is covered in snow during the winter season [15]. Some of the most notable studies on snow levels include those by [16] in the northern zone and their study of how wind impacts snow cover duration patterns in Pascua Lama and the Andes through the use of ablation simulation models. At mid-latitudes, Masiokas et al. [17] analyzed the central region of the Andes and identified a positive correlation with warm phase of the El Niño/Southern Oscillation (ENSO), showing the importance of western air masses in the regulation of mountain snowfall. Between the northern and southern zone, Stehr et al. [18] have studied spatiotemporal dynamics of snow cover over five Andean watersheds, detecting an important decline of snow cover. In the southern zone, [19] have studied the snowpack on the Northern Ice Field using satellite images and meteorological data. However, there is a lack of research on the amount of snowfall in the region of Aysén over the past few decades and how this variation might impact the water scheme of the region's most important rivers, which have been under consideration for several hydroelectric projects [20]. Likewise, this variability might affect how water is managed regionally for different uses, like drinking water and agriculture, thus leading to unsustainable water-use conditions under future scenarios of lower resource availability. In fact, a decrease in precipitation has been estimated for the coming decades [21,22], which might potentially lead to greater vulnerability and might require water users' organizations to be created for resource use under a scenario of scarcity. It is for this reason that the main objective of this work is to analyze the spatial and temporal variability of snow cover in the Aysén river basin using satellite images and in situ measurements. The structure of this paper is as follows: Sections 2 and 3 detail the study area and data used. Section 4 presents the methodology, Section 5 the results and, lastly, Sections 6 and 7 provide a brief discussion and conclusions.

Study Area
The Aysén river basin ( Figure 1) is located in the XI Region of Aysén del General Carlos Ibáñez del Campo between 45 • and 46 • south latitude. It has a surface area of 11,540 km 2 , of which 38.3% is covered by native forests, 22.6% by scrub and meadows, and 13.6% by glaciers and permafrost [23].
Its climate is cool oceanic with low temperatures, heavy precipitation and strong winds. However, the different terrain within the basin creates micro-climates: the western sector is influenced by islands and archipelagos and the western slope of the Patagonian Andes, with average annual precipitation between 3000 and 4000 mm and an average annual temperature varying between 8 and 9 • C. The central part of the basin presents a cold steppe climate on the eastern slope of the Patagonian Andes. On the eastern side, rainfall drops as far as 621 mm per year in Balmaceda and down to 1385 mm per year in the city of Coyhaique [24]. Its climate is cool oceanic with low temperatures, heavy precipitation and strong winds. However, the different terrain within the basin creates micro-climates: the western sector is influenced by islands and archipelagos and the western slope of the Patagonian Andes, with average annual precipitation between 3000 and 4000 mm and an average annual temperature varying between 8 and 9 °C. The central part of the basin presents a cold steppe climate on the eastern slope of the Patagonian Andes. On the eastern side, rainfall drops as far as 621 mm per year in Balmaceda and down to 1385 mm per year in the city of Coyhaique [24].
The basin has an elevation ranges from the sea level to a high of 2227 m and an average slope of 32%. The western part is where the Andes mountain range intersects and the highest elevations and steepest slopes are found, and where the typical vegetation is temperate, lenga beech forest (Nothofagus pumilio) [25]. The eastern part, on the other hand, has wide valleys and gentle slopes with the Patagonian steppe [26].
The river basin's predominant hydrological regime is pluvial-nival ( Figure 2), with three subregimes: mixed nivo-pluvial in trans-Andean riverbeds, mixed nivo-pluvial and nivo-pluvial regulated by lakes and glaciers. The first sub-regime is characterized by flow coming mainly from melting; the second by the predominance of rainfall as a result of orographic factors; and the third, by the regulation of lakes and glaciers. In this last sub-regime is the Blanco River and 10 other lakes that produce an increase in flows during the dry season, in addition to a several peaks with permafrost [24]. The major tributaries of the main river are the Mañihuales, Simpson, Blanco, Coyhaique and Los Palos rivers, which provide the basin with an average flow of between 250 and 800 m 3 ⋅s −1 [27,28]. The basin has an elevation ranges from the sea level to a high of 2227 m and an average slope of 32%. The western part is where the Andes mountain range intersects and the highest elevations and steepest slopes are found, and where the typical vegetation is temperate, lenga beech forest (Nothofagus pumilio) [25]. The eastern part, on the other hand, has wide valleys and gentle slopes with the Patagonian steppe [26].
The river basin's predominant hydrological regime is pluvial-nival ( Figure 2), with three sub-regimes: mixed nivo-pluvial in trans-Andean riverbeds, mixed nivo-pluvial and nivo-pluvial regulated by lakes and glaciers. The first sub-regime is characterized by flow coming mainly from melting; the second by the predominance of rainfall as a result of orographic factors; and the third, by the regulation of lakes and glaciers. In this last sub-regime is the Blanco River and 10 other lakes that produce an increase in flows during the dry season, in addition to a several peaks with permafrost [24]. The major tributaries of the main river are the Mañihuales, Simpson, Blanco, Coyhaique and Los Palos rivers, which provide the basin with an average flow of between 250 and 800 m 3 ·s −1 [27,28].

Remote Sensing
The Moderate Resolution Imaging Spectroradiometer (MODIS) is a terrestrial remote-sensing

Remote Sensing
The Moderate Resolution Imaging Spectroradiometer (MODIS) is a terrestrial remote-sensing system which employs a cross-track mirror and a set of individual sensors to provide images of the Earth surface and clouds in 36 spectral bands in a wavelength range between 0.405 to 14.385 µm. The main purpose of MODIS is to aid in the study of vegetation and its properties, the concentration of aerosols, albedo and surface temperature, as well as snow and ice cover. Snow cover is identified using the Normalized Difference Snow Index (NDSI) [11,29], while snow cover over densely forested areas is identified using thresholds of NDSI along with the Normalized Difference Vegetation Index (NDVI) [30]. MODIS snow cover products are provided daily (MOD10A1) and as 8-day composites (MOD10A2), both at a resolution of 500 m with a sinusoidal map projection. MOD10A2 is a composite of MOD10A1 generated to show maximum snow cover extent [31]. For this study, MOD10A2 (V005) was used to reduce the effects of cloud cover on the study area. The sample period covers the years between 2000 and 2016, with a total of 770 images.

In Situ Information
In order to compare variations in the snowpack with hydro-meteorological parameters, daily records of average flows (m 3 ·s −1 ), average air temperature ( • C) and accumulated precipitation (mm) from stations belonging to the General Directorate of Water (Dirección General de Aguas, DGA) and the Meteorological Directorate of Chile (Dirección Meteorológica de Chile, DMC) were used. Table 1 shows the hydro-meteorological stations and their data availability during the time period.

Remote-Sensing Processing
To cover the study area, two MOD10A2 (h12v13, h13v13) tiles for the period between 2000 and 2016 were acquired from the National Snow and Ice Data Center (NSIDC) (https://earthdata.nasa. gov/). The images were reprojected into the Universal Transverse Mercator System (UTM) Datum WGS-84 Zone 19S. The original values of the MODIS product were classified into three classes: snow, no-snow (land, continental bodies of water and ocean) and clouds, where the latter corresponded to pixels with no information corresponding to clouds and classes with missing values and errors generated by the sensor (Table 2).

Snow Covered Area (SCA)
The snow cover corresponding to the reclassified product's snow class was converted to area by multiplying by a conversion factor of 0.25 km 2 . This area was presented as snow covered area (SCA).

Seasonal Trend Analysis
In order to analyze snow cover variation, the average monthly anomalies in SCA were calculated using the method by [32] (Equation (1)): where xa i is the SCA anomaly, x i is to the annual/seasonal SCA value, and x is the average for the seasonal/annual period for the years in the time series. The trend analysis was performed on the anomaly estimation indicated above, in which the slope among all the ordered pairs on a scatter plot was calculated using Sen's slope estimator (Equation (2)). In the time series, time will be on the x-axis and SCA on the y-axis. Then, the median of the created data set (m q ) was determined and the result represented the non-parametric trend of the sampled time series [33,34], where m q is the series of slopes among all the combinations of ordered pairs which are natural numbers from 1 to k (Equation (3)): where n is the total number of ordered pairs. To estimate whether there is a significant trend in the time series, the non-parametric Mann-Kendall test was performed [35,36] since it has been identified as one of the most robust techniques available for estimating trends in environmental variables [37]. Equation (4) gives the statistics on the sign of the trend for the sampled series, where a positive value of t indicates an increasing trend and vice versa: where signo(x) extracts the sign of the expression entered (Equation (5)).
When the t statistic is calculated with more than 10 data points, it presents a normal distribution, which is what makes it possible to run it through a hypothesis test, for which the z statistic must be calculated as a function of t and its variance (Equation (6)): To calculate the z-statistic, Equation (7) was used, The significance level for this work is 0.05, and because it is bilateral, the rejection region of the t statistic will be [Z 0.025 ; Z 1-0.025 ].

SCA and In Situ Data
The stations whose data availability was equal to or greater than 85% in the period 2000-2016 were used. Hydro-meteorological data were added annually and measurements from the stations were averaged to obtain a representative value at hydrological basin level. As far as precipitation and average air temperature, the annual cumulative and the annual average were used, respectively, as in Equation (8): where T is the annual mean temperature at basin level, Ti is a specific weather station, i is the air temperature measurement and n is the number of stations in the basin. P is the average accumulated precipitation at basin level and Pi is the annual accumulated precipitation measured by one weather station. As far as stream-flow, the annual average of the gauging station at the outlet of the basin was used.           The annual mean SCA shows a non-significant decreasing trend at a rate of −20.01 km 2 ·year −1 at 95% confidence. The autumn and winter seasons show decreasing rates of −14.92 and −19.44 km 2 ·year −1 with no statistical significance (p > 0.05). In addition, in winter, beginning in 2014, there are negative anomalies of the annual mean SCA of about 1000 km 2 below the mean for the 2000-2016 period. Likewise, the spring and summer seasons show non-significant decreasing rates with a slope of −4.60 and −9.99 km 2 ·year −1 , respectively.  Seasonal trends of mean air temperature anomaly are positive over the study area ( Figure 8). For autumn, the time period presents a statistical significance trend (+0.10 °C⋅year −1 ). It can be seen that positive air temperature anomalies are in relation to the negative SCA anomalies. This effect could be related to the global warming causes, specifically an increase of the 0 °C isotherm, resulting in decreased snow cover in mountainous regions or a reduction in winter snowfall, while also pushing Seasonal trends of mean air temperature anomaly are positive over the study area ( Figure 8). For autumn, the time period presents a statistical significance trend (+0.10 • C·year −1 ). It can be seen that positive air temperature anomalies are in relation to the negative SCA anomalies. This effect could be related to the global warming causes, specifically an increase of the 0 • C isotherm, resulting in decreased snow cover in mountainous regions or a reduction in winter snowfall, while also pushing early the thaw in spring [38,39].

Discussion
The MODIS daily snow product indicates an accuracy of over 90% identifying snow cover under clear sky conditions and its observations are matched to measurements made in situ or taken by other satellites [31,[40][41][42][43]. Likewise, in mountain areas with the presence of evergreen forest, the precision of the product decreases up to 75-80% for winter months [44], because thick canopy does not allow the snow to be identified [45]. In this research, a decreasing trend in the annual and seasonal mean of SCA was determined in the Aysén river basin in the period 2000-2016, although with no statistical significance. However, it is necessary to validate these results in future work by incorporating data about snow height and cover measured in situ, which are currently non-existent for the study area. This would improve the estimation of the amount of snow, and not just its surface area, thus improving available estimates.
The last five years in the time period show negative anomalies in the precipitation, streamflow and SCA, with 2016 as the year with the most negative set of anomalies. The variations in snow cover are closely related to the stream flow hydrological response in seasonal cycles. In this sense, generally a positive anomaly of SCA corresponds to a positive anomaly in river discharge in the periods of snow fusion [46]. The SCA shows wide intra-and inter-annual variability, which is consistent with research done in other studies on central Chile and Asia [17,47,48]. In fact, compared to results obtained by [18] in Chile's Andean basins between latitudes 32.0-39.5 • S, there is a decreasing trend in SCA and no significance in the southernmost basins, which is consistent with SCA variability in the Aysén river basin. However, [17] identified a non-significant increasing trend of snowpack over central Chile in the period 1951-2005, which differs from the results of this study due to the fact that it comprises a more extensive data series (54 years) and does not include austral latitudes.
According to [21], a drop in SCA is expected over the Andes mountain range, as well as an increase in the intensity and frequency of extreme climate events related to prolonged droughts and floods. In line with this, central Chile recently experienced the so-called mega-drought in the 2010-2015 period, with a deficit of between 25% and 45% in annual rainfall [49]. During the period from 2000 to 2016, there were five occurrences of El Niño or warm phase of ENSO (2002, 2004-2005, 2006-2007, 2009 and 2015-2016) and three of La Niña (1999)(2000)(2007)(2008) and 2010-2011) (http://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php). In keeping with this, the maximums in annual mean SCA coincide with the periods of La Niña, with the exception of the year 2005. El Niño events are often associated with increases in precipitation, but can also increase the air temperature, which together can have a pronounced altitude effect on the snow cover during spring and fall [50]. In this sense, an increase in air temperature causes 0 • C isotherm to ascend to higher elevations resulting in a larger proportion of precipitation falling as rain as opposed to snow [51,52]. However, due to the short time series of sampled satellite images, it would be necessary to carry out complementary studies to confirm the aforementioned relationship.
The precipitation data used show a limited spatial distribution with respect to the diversity of geographic zones in the study area. In fact, most of the stations are located in sectors near the mouth of the Aysén River, making evident the lack of data in areas of snow accumulation, height, lakeside areas, among others. Likewise, precipitation stations consist of basic rain gauges lacking any adaptations for measuring or determining snowfall. This means that there could be a bias in the data when combining solid and liquid precipitation at the same station. In addition, this study analyzed SCA and not the height or volume of snow, so a greater SCA does not necessarily mean a greater amount of water stored as snow [18].
In terms of management, variability and decreasing precipitation and SCA might potentially require a water-management plan for the basin because a negative trend in snow cover would mean a relatively greater shortage of available water in the basin, thus affecting ways of dealing with a demand that has not been prioritized as of yet in regional management plans and programs. In fact, the region is adapted to a pluvio-nival hydrological regime, and therefore water availability is abundant in the dry season. This means that the greater water demand is offset by an abundance of the resource. However, a decrease in the snowpack might impact the pluvio-nival regime by altering the water supply's seasonality. In fact, forecast studies show that in the basins of central Chile, significant changes have already been detected in the regimes of the pluvio-nivales rivers, thus leading to a change in its status to pluvial, modifying the time balance between water supply and demand, affecting agricultural and human consumption [53].

Conclusions
In this study, an initial approach is presented for evaluating snow cover variability in the Aysén river basin in Chilean Patagonia. Using MODIS satellite images, the region could be monitored over the last 17 years and SCA trends were estimated. SCA values showed a significant decrease between 2010 and 2016 as compared to the historical record, totaling 1050.48 km 2 for the entire Aysén river basin. The decreasing trend in snow cover is dropping, −20.01 km 2 ·year −1 (p > 0.05), and variables such as temperature were estimated to be positive, +0.07 • C·year −1 . Seasonal trends for precipitation and stream flows were negative estimated by −8.66 mm·year −1 and −23.78 m 3 ·s −1 ·year −1 for spring, respectively. Despite the temporal composite method applied to the MODIS imagery, the persistence of cloud cover during a winter might affect the estimates of SCA and, therefore, these retrievals could impact on seasonal SCA trends. Moreover, the data's scarce spatial representativeness and lack of stations over mountains impact on the validation of the SCA trend and its relation to precipitation, stream flows and air temperature. Nevertheless, this is the first work to contribute to the analysis of snow cover in one of Chilean Patagonia's most important river basins, where the impact of global warming could have consequences in the water balance and, therefore, impact on the sustainable management of water resources.
Author Contributions: T.P. and C.M. conceived and designed the research. T.P. performed the data processing and analysis. C.M. and R.F. assisted with interpretation and data analysis. All authors contributed to the manuscript write-up and editing.