Remote Sensing Study of Glacial Change in the Northern Patagonian Icefield

The Patagonian Icefield has the largest temperate ice mass in the southern hemisphere. Using remote sensing techniques, this study analyzed multi-decadal glacial retreat and expansion of glacier lakes in Northern Patagonia. Glacial boundaries and glacier lake boundaries for 1979, 1985, 2000, and 2013 were delineated from Chilean topographic maps and Landsat satellite images. Aster stereo images were used to measure mass balance from 2007 to 2012. The highest retreat was observed in San Quintin glacier. The area of glacier lakes increased from 13.49 km2 in 1979 to 65.06 km2 in 2013. Four new glacier lakes formed between 1979 and 2013. Between 2007 and 2012, significant glacial thinning was observed in major glaciers, including HPN1, Pared Norte, Strindberg, Acodado, Nef, San Quintin, Colonia, HPN4, and Benito glaciers. Generally, ablation zones lost more mass than accumulation zones.


Introduction
Glaciers are key indicators for assessing climate change [1]- [3].Beginning in the nineteenth century, glaciers in many parts of the world retreated significantly, which was a clear indicator of climate warming [3]- [7].Therefore, an accurate estimation of mass balance, glacial retreat, and glacier lake expansion in mountain glaciers will contribute to the understanding of the global impacts of climate change.Such data will also help develop an early warning system for Glacier Lake Outburst Floods (GLOF), which can devastate vulnerable downstream localities.
Further, documentation of glacier measurements in major glacial regions is not comprehensive.For example, the European Alpine glaciers have been well documented using both in situ and remote sensing methods.The Himalayan glaciers have also attracted more attention in recent decades.Conversely, South American glaciers are not studied as intensively.Like other glacier regions, the reduction in glacier mass is also very evident in Patagonia [8].However, documentation of the status of glaciers in the Chilean Andes has been very limited.
The glacier mass balance is the total loss or gain in a glacier's mass over the course of a hydrological year [9].There are both direct and indirect methods to measure mass balance.The direct method involves an in-situ collection of glacier data to measure retreat and mass balance taken at lower ablation regions [10] [11].In recent years, remote sensing satellite images have helped researchers study even the world's most inaccessible glaciers [12]- [15].The indirect method of mass balance measurement requires the use of remote sensing techniques, which involves the use of satellite images, aerial photos, and photogrammetric methods [15] [16].
In this study, we first analyzed the spatial and temporal positional shifts in the terminus of the glaciers in the North Patagonian Icefield (NPI).We also analyzed the expansion of glacier lakes from 1979 to 2013.Then, using Aster satellite stereo images from 2007 to 2012, we measured the mass balance of some of the major glaciers in the NPI.

Study Area
The Chilean cryosphere can be divided into four parts based on climate and glacier characteristics: 1) the Arid North (18˚S to 32˚S); 2) the Central Andes (32˚S to 38˚S), which has high precipitation in the winter and dry summers; 3) Southern Chile (38 to 41˚S), which has an active volcano region; and 4) the Northern and Southern Patagonian ice fields (41˚S), which has abundant precipitation [14].The NPI and the Southern Patagonia Icefield (SPI) together have the largest temperate ice mass in the southern hemisphere [17].The temperatures of these temperate glaciers are always at melting point throughout the glacial area [18].These glaciers respond to climate changes after only short periods [18].Therefore, such glaciers must be studied in order to understand the dynamics of climate change.
Patagonian glaciers have long been neglected for mass-balance studies, although researchers have undertaken other glacial studies in the area [17].Due to the remoteness, inaccessibility, and bad weather of Patagonia, field-based mass-balance studies are lacking [18].Because glaciers in this region have not been as closely studied, researchers are uncertain about the degree that these glacier-melts contribute to sea level rise [19].Glaciers are the most prominent water source in this area; therefore, the lack of records of water balances and annual mean discharges in the region is of great concern [20].In recent years, however, several studies have been initiated to study glacier changes in this region [21] [22] and have clearly shown glacial retreat and thinning in the Patagonian ice field.
The NPI is in the Central Andes region [23].The NPI receives abundant precipitation: about 2000 to 11,000 mm of water equivalent per year [19].The eastside of NPI, however, receives precipitation of only approximately 1500 mm per year [24].The precipitation on the east side is predominantly attributable to mid-latitude storms that move westward between 40˚S and 50˚S, as well as the mean latitudinal locations of the storm tracks, which closely follow the axis of the westerly winds in the middle and upper troposphere [25].
The west wind becomes very strong and nearly permanent when the mean air temperature at sea level decreases gradually from 13.7˚C at latitude 37˚23'S to 6.5˚C at latitude 53˚10'S [23].The NPI is also characterized by the Föhn winds, which are dry, warm, and down-sloping winds.These winds lead to dry conditions and high temperatures on the eastern side of the icefield.They also increase the precipitation gradient in the west-east direction [26].But precipitation in the NPI shows no significant seasonality [27].The elevation of the limit between rainfall and snowfall during precipitation events is the main climatic factor that causes glacier fluctuations in the NPI [23].
The NPI is located between 46˚30'S and 47˚30'S, and it stretches approximately 100 km along 73˚30'W, covering a total area of 4200 km 2 (Figure 1).It has a width of approximately 71 km in the west-east direction, from the tongues of San Quintin glacier to Soler glacier [21].Most of the glaciers lie at altitudes of 1000 m in the west and 1500 min the east.A northwest-to-southeast line of nunataks divides the northeast icefield from the southwest icefield.The NPI has twenty-eight outlet glaciers [17].San Quintin and San Rafael glaciers are the largest glaciers in the NPI, each with an area of about 760 km 2 , followed by Steffan and Colonial glaciers [28].San Rafael is the only tidewater-calving glacier in the NPI.Most of the other glaciers have freshwater-calving fronts [29].Both calving-glacier types respond directly to climate change [30].
Equilibrium Lines of Altitude (ELA) of about 900 to 1350 m separate the NPI into a 2578 km 2 accumulation zone and a 1550 km 2 ablation zone [28].The summit of Mount San Valentin, at the northeast end of the ice field, has the highest elevation, with 4032 ± 1 m above sea level.The tidal water of Laguna San Rafael, located on the western side of the ice field, has a minimum altitude at sea level [21].Due to climatic conditions, topography, and the nunataks that divide the icefield, the western outlet glaciers are larger and more active than the eastern glaciers [21] [31].Over the last 60 years, [32] indicated a general retreat of 21 glaciers in the NPI.After the 1990s, there was an increase in the rate of retreats, with the western side experiencing greater retreat than the eastern side [33].Glacier Lake Outburst Floods (GLOF) are a disastrous discharge of large volumes of water, caused by the breaking of moraine-dammed lakes [34].Glacier lakes form from melting glacier ice and snow.Higher rates of melting have increased the threat of glacier lake outburst in the world's high mountain ranges [35].Various natural phenomena, such as avalanches, earthquakes, or even the weight of the water itself, break the moraine dams.Some of the most devastating GLOF occurred in the Peruvian Andes.Researchers [36] have described the natural hazards associated with the GLOF of Palcacocha Lake in Cordillera Blanca.In 2000, a rock fall partially displaced a moraine-dammed lake in Lago Calafate glacier and caused a GLOF in the Northern Patagonia Ice field [37].

Data and Methods
We acquired topographic maps of the study area with a 1:50,000 scale created by the Istituto Geografico Militar of Chile (IGM) in 1979.The contour interval in these maps is 50 m.The vertical random error of these maps is 17 m [38].The National Imagery and Mapping Agency [39] estimated the horizontal error in the maps as 15 m.The maps were projected to WGS 1984 UTM Zone 18S.
To delineate the glaciated area from the non-glaciated area, we digitized the glacier boundaries on the ASTER 2012 and 2007 images in Arc GIS 10.1.We digitized the region that overlapped between the 2012 and 2007 images and used the resulting image for the mass balance analysis.Fourteen glaciers fell within the overlapping region: San Quintin, Fraenkel, Strindberg, Colonia, Acodado, HPN 1, HPN-4, Cachet, Cachet Norte, Nef, Arco, Benito, Pared Norte, and U5 (Figure 1).
We used three Landsat images and 1:50,000 topographic maps to measure the retreat of the NPI glaciers and to analyze the expansion of glacier lakes.The Landsat images were from 1985 (Landsat 5 TM), 2000 (Landsat 7 ETM+), and 2013 (Landsat 8).We downloaded the Landsat images from the USGS EROS website (http://earthexplorer.usgs.gov).Both image types were re-projected to WGS 84 UTM Zone 18S.
To clearly delineate the glacier lakes, we created the Normalized Differenced Water Index (NDWI).We also visually inspected the images of the lakes using various band combinations.The NDWI uses (equation 1) -the (Near Infrared) NIR and Blue band to differentiate water bodies from other features [35].
Using the NDWI and the Landsat images, we digitized the glacier lakes in the NPI and in the topographic map.We observed a significant increase in the area of the glacier lakes between 1979 and 2013.
For the mass balance estimation, we used two ASTER images, one from April 4, 2007 and one from March 18, 2012.We used ASTER satellite data because it has the stereoscopic capability to generate Digital Elevation Models (DEMs) [40].ASTER bands 3N and 3B can be used in a DEM creation.Bands 3N and 3B are the Nadir and the backward-looking telescopes on the NIR wavelength [40].These bands have a wavelength range of 0.78 to 0.86 µm.For this study, we geo referenced each of the ASTER images to WGS84 datum, UTM Zone 18S.

DEM Creation
We used the Ortho Engine extension of PCI Geomatica 2012 software to create DEMs.To create DEMs of 15 m resolution, we collected 50 well-distributed Ground Control Points (GCPs) and 100 well-distributed Tie Points (TPs) from the 2007 and 2012 images.During the GCPs and TPs collections, we gave priority to non-glaciated regions, because, in the absence of a natural disaster, such regions hardly change.Glaciated areas, shaded areas, and cloud-covered areas were excluded from GCPs, as such areas are highly variable and could introduce inaccuracies to the DEM.The root mean square error (RMSE) after the GCPs collection was 2.14 m for the 2012 image, with an X RMSE of 1.81 m and a Y RMSE of 1.13 m.The RMSE for the 2007 image was 1.37 m, with an X RMSE of 1.03 m and a Y RMSE of 0.90 m.We computed the change in elevation values by finding the difference in pixel-by-pixel DEM values between the 2012 image and the 2007 image.We converted the elevation difference image to points.We excluded shaded areas, steep slopes, and clouded areas from the analysis.The density of ice used to calculate mass balance was 900 kg•m −3 [9].The mass balance estimate was presented as meter water equivalent per year.The elevation of the master DEM (2007) was compared with the elevation of the 1:50,000 topographic maps.The correlation of over 70 well-distributed random points from topographic map elevations and the DEM elevations was very high (b = 0.91).Considering the lower resolution of the topographic map, DEM values seem to be highly accurate.
We expected the non-glaciated region to have a zero-elevation difference value, because there should be no change in elevation unless there is a natural disaster or a vertical error in the DEM.We used elevation differences between 2007 and 2012 DEMs in non-glaciated areas to calculate the uncertainty in the mass balance estimation in the glaciated areas.
We computed uncertainties in the elevation difference in glacier and non-glacier areas using individual Standard Error (SE) [41].We calculated SE in the non-glaciated area using the equation 2, where n is the number of included pixels.

STDV non glac
In accordance with the law of error propagation, we used the SE and the mean difference of elevation (MED) of the non-glaciated area to calculate uncertainty (equation 3) in glaciated areas [41].
We measured retreat by the change in position of the tongue of the glacier during the study period [21], as seen in Figure 2.
Table 1 shows the extent of retreat of all nine glaciers.San Quintin retreated the most, 4.5 km.This glacier tongue was intact in 1979 but by 2013 had retreated to form a new lake.Another glacier tongue in San Quintin retreated by 3.03 km.San Rafael retreated about 3.3 km between 1979 and 1987 [21].San Quintin, conversely, retreated by only approximately 0.5 km within the same time period.In general, all glaciers in the study area retreated, which is attributable to atmospheric warming in the region [21].The rate at which these glaciers retreat is a direct response to climate, because these are temperate glaciers and they respond rapidly to change in temperature [23].Due to the lack of weather data in the region, however, we could not analyze long-term temperature changes.
The number of glacier lakes in the study area increased significantly between 1979 and 2013.In 1979 there were 9 glacier lakes; in 1985 there were 10; in 2000 there were 11; and by 2013, there were 13 glacier lakes.The glacier lake area increased from 13.49 km 2 in 1979 to 65.06 km 2 in 2013 (Figure 3).
San Quintin Glacier Lakehad the maximum expansion of area, from 5.42 km 2 to 42.37 km 2 .This expansion is the result of the creation of a new glacier lake at the tongue of San Quintin.A study [42] observed a similar trend between 1945 and 2011.They also found that this drastic increase in San Quintin Lake compared to other glaciers.Table 2 shows the number of lakes and their total areas.
In the non-glaciated area, we selected 522,513 random points to estimate the uncertainty in the mass balance estimation.The average elevation difference in the non-glaciated area was 2.82 m, with a standard deviation of 8.99 m and a standard error of mean of 0.012 m.Based on these parameters, we estimated uncertainty in mass balance of the glaciated area to be ±0.51m.w.e.a −1 .In glaciers, we estimated mass balance for accumulation and ablation zones separately.Table 3 presents mass balance values from 2007 to 2012.
To estimate mass balance, we used only the areas of glaciers that fall within the boundaries of the satellite images.In general, ablation zones had higher loss of mass than accumulation zones (Table 3).The HPN1 glacier ablation zone showed the highest loss of mass, with −4.76 ± 0.51 m.w.e.a

Discussion and Conclusions
The results of this study clearly indicated general retreat, expansion of glacier lakes, and loss of mass in NPI.However, other studies have found the advancement of glacial tongues in some glaciers in the previous decade.For example, San Quintin glacier advanced in 1993 and then retreated between 1996 and 2000 [43].Our study showed that the retreat is continuing in San Quintin.Overall, there was a 4.5 km retreat in the San Quintin glacier tongue between 1979 and 2013.A prior study observed the retreat of San Rafael glacier stopped in 1991 [44].However, that study found a 3.92 km retreat in San Rafael glacier during the study period.
A study [29] observed a mass balance loss in all the glaciers in the NPI between 1975 and 2001.These average ice-loss rates have more than doubled to an equivalent sea level rise in the North and South Patagonian Ice-fields between 1995 and 2000 [19].Similarly, this study observed negative mass balance in most glaciers.This study also found that accumulation zones are losing mass at a higher rate than ablation zones.
Glaciers on the western side of the nunataks experienced more retreat than glaciers on the eastern side of the nunataks.Another study [28] referred to orographic factors as the cause of this pattern.He further explained that the eastern side of the NPI receives less snowfall compared to the western side, because the NPI lies in the zone of the westerlies.The eastern side thus receives more snow from the western side through drift snow that is blown by wind.Due to this drift snow, the eastern side accumulates more snow and experiences less retreat in the NPI.A study [45] estimated the thinning rate of HPN-1 glacier as −0.7 m.w.e.a −1 between 1945 and 1975; another study [29] estimated the thinning to be −4.0 m.w.e.a −1 between 1975 and 2001; and we recorded a negative mass loss of HPN-1 to be −4.76 ± 0.51 m.w.e.a −1 in the ablation zone and −0.80 ± 0.51 m.w.e.a −1 in the accumulation zone.
Although these results were not validated with in-situ data for either mass balance or retreat, similar results were observed by [18] [19] [28] [45] for mass balance studies and [21] [29] [41] for glacier retreat and lakes.The results of this study and previous studies observed overwhelming retreat of glaciers and loss of mass in the NPI.These results must be taken seriously, as these conditions are linked to climate fluctuations in the region.

Figure 2 .
Figure 2. Retreat of nine major glaciers in the NPI.

Table 1 .
Retreat of glaciers between 1979 and 2013.

Table 2 .
Glacier lakes area for 1979 and 2013.

Table 3 .
Mass balance of glaciers between 2007 and 2012.
Note: Area is not the actual area of the glaciers, but it represents only the area that was covered in the image.