Glacier changes on Sierra Velluda massif , Chile ( 37 ◦ S ) : mountain glaciers of an intensively-used mid-latitude landscape

Abstract. The central-southern section of Chile is defined as one of the Latin American hot spots in the last IPCC Report due to the impact of glacier retreat on water resources, the transitional character of the climate, and its importance in terms of agricultural and forestry activities. In order to provide a better understanding of glacier behavior in this zone, this paper analyzes the volumetric changes of glaciers in the Sierra Velluda, located in the upper Bio Bio River Basin. Bibliographic sources, satellite images, and DEMs were used to estimate frontal, areal, and volumetric changes. An analysis of significance was performed in order to provide accurate estimations of the fluctuations. The results indicate that Sierra Velluda glaciers have suffered a significant reduction since the 1960s, despite some short periods of positive fluctuations. A maximum position of a glacier for the year 1828 was identified, which is in concordance with other proxies registered elsewhere in Chile. These changes agree with measurements of glacier fluctuation elsewhere in Chile. While short-term fluctuations are consistent with the inter-annual precipitation variability, lake levels records, and a warm phase of the El Nino Southern Oscillation (ENSO), the general shrinkage agrees with the shift of the ENSO (PDO) in 1976. Therefore, it is proposed that Sierra Velluda's glaciers are highly sensitive to high frequency climatic fluctuations and even to inter-annual variability. Considering that models project a reduction in Andean precipitation and an altitudinal increase in the 0 °C isotherm, these ice bodies are expected to continue to shrink.


Introduction
South-central Chile is identified as one of Latin American's area of particular concern in the most-recent IPCC Report (Magrin et al., 2007) due to the impact of glacier retreat on water resources.Due to its hydrologic potential (Mardones, 2001), this zone has been the principal source of hydroelectric energy for most of the Central Interconnected System in Chile.In fact, a series of 10 hydroelectric dams presently feed the system  (CDEC-SIC, 2009).Additionally, due to the transitional character of the climate (Devynck, 1970), this zone is crucial for the development of studies on climate change in mid-latitudes and its influence on human activities.This zone is fundamentally dedicated to agricultural (54% of the nation's total) and forestry (58% of the nation's total) activities (INE, 2007).
Although there are detailed glacier inventories for some river basins in the zone between 32 • to 41 • S (Valdivia, 1984;Rivera, 1989), they have not been updated and the scarce extant data on volumetric fluctuations is concentrated between 33 • and 34 • S (80%).Consequently, their response to current changes in temperature and precipitation is unknown.
In order to provide a better understanding of glacier behavior in this zone, this paper analyzes the volumetric changes of glaciers in the Sierra Velluda (37 • 27 S and 71 • 24 W), which are located in the upper Bío Bío River Basin and inside the Laguna del Laja National Park in central-southern Chile (Fig. 1).With a maximum altitude of 3585 m, it is the highest peak in the basin.Its genesis has been dated at 495 000 y BP (Vergara et al., 1985).In addition, evidences of both Last Glacial Maximum and Holocene fluctuations have been documented (Mardones and Jaque, 1991).

Frontal and areal changes
Both bibliographic sources and satellite images were used to measure changes to the glacier's surface and frontal positions.Two bibliographic sources were analyzed.The first one corresponds to the description and drawings of the exploration led by Eduard Poeppig, which is an extract of a more extensive trip made in Peru and Chile between 1826 and 1829.The direct translation from German by Keller (1960) was used.Since this document includes a series of contemporary commentaries on the context in which the exploration was made, an adjustment of "epoch-specific" meaning, as has been 687 Discussion done elsewhere (e.g.Araneda et al., 2009), was considered unnecessary.To estimate the position of the glacier tongue, the functional landscape analysis method was used (de Bol ós, 1992).
The second source corresponds to the glacier inventory performed by Rivera (1989), who mapped and described 11 glaciers corresponding to the most glaciated area of the Bío Bío River Basin.The map of this inventory was based on the interpretation of OEA (Organizaci ón de Estados Americanos) flight aerial photographs (1961), with a nominal scale of 1:60 000, as well as the 1:50 000 topographic maps of the Chilean Instituto Geogr áfico Militar.Rivera's map, originally in paper format, was created using the Datum SAD 1969.In this study, that map was digitalized in a Geographic Information System (ArcGIS 9.2), and subsequently transformed into Datum WGS 84 using standard parameters (NIMA, 1997).The nomenclature used to identify each glacier follows Rivera (1989), which was based on the specifications of M üller et al. (1977) and the national coding system for river basins (e.g.DGA, 2004).However, since the Sierra Velluda presents four watersheds which contribute to the Laja River (a tributary of the Bio Bío River) the nomenclature was expanded by 2 digits in addition to the code for the Bío Bío River Basin (083, DGA, 2004).Additionally, each glacier inside a subbasin of the Laja River was labeled accordingly to its location from north to south.This procedure gave a label of 2 letters and 7 numbers.
Satellite images were also used.First, a scene search was performed using two internet servers.The main criterion for image selection was their availability for the end of the austral summer (between late January and late March), at the end of the annual ablation period (Table 1).The Landsat ETM+ images downloaded from GLCF were already ortho-rectified to WGS 84 UTM zone 19 south and were used to co-register the other Landsat MSS, Landsat TM and Aster scenes used in this study.Nevertheless, these were checked against the topographic maps and found to be precise within one pixel (30 m).To correct for reflectance of these scenes, processing software (ENVI 4.2) was used to convert the digital numbers into reflectance values using the header file.
Additionally, the Terra-Aster sensor images had two co-register procedures.The first one used the orbital parameters of the header file for geometrical correction using ENVI 4.2.Secondly, more than ten control points in the Landsat images and their corresponding points in each Aster scene were identified.The Aster scenes were resampled using the cubic convolution method.Since these images were acquired in Level 1A, atmospheric and reflectance corrections were needed in order to convert the raw data to radiance.This task was performed using the metadata of each scene and the FLAASH tool (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) in ENVI 4.2 (ITT, 2009).
Various remote-sensing techniques are used to determine glacier cover depending upon factors such as the type of sensor used, the spatial resolution, and especially the spectral resolution (for more details see Rees, 2006).Landsat MSS images can be used to delimit glaciers using false color composites (Ye et al., 2006).In this case, an image with the bands MSS4, MSS3 and MSS2 were composed and superimposed on a Digital Land Model (SRTM90V.4)to generate three-dimensional views.The glaciers were delimited manually on the screen, using the software ArcGIS 9.2, as done in other studies (Ye et al., 2007;Fern ández et al., 2010).Since the study area was visited prior to the delimitation in order to pre-classify these images, the procedure can be considered as supervised classification.
Band ratios were used for the Landsat and Terra-Aster images.In the contemporary literature, a series of ratios have to be used, and each author argues that the one used worked optimally for the purpose of each study.Paul et al. (2002) argues that the optimal solution is built from the bands TM4 (0.76-0.90 µm) and TM5 (1.55-1.75µm) (Hall  , 1988), equivalent to ETM+4 (0.76-0.90 µm) and ETM+5 (1.55-1.75µm) and to Aster 3 n (0.78-0.86 µm) and Aster 4 (1.6-1.7 µm), while Andreassen et al. (2008) used the bands TM3 (0.63-0.69 µm) and TM5 (1.55-1.75µm) (Rott, 1994), equivalent to ETM+3 (0.63-0.69 µm) and ETM+5 (1.55-1.75µm) and to Aster 2 (0.63-0.69 µm) and Aster 4 (1.6-1.7µm).Here, both ratios were compared using a Landsat ETM+ scene of the year 2001, the field classification, and the information existing in the inventory of Rivera (1989), who calculated a total surface area of 20.32 km 2 using aerial photos taken in 1961.The application of the ratio of Hall et al. (1988) provides a value of 46.62 km 2 , while 25.40 km 2 is calculated using the ratio of Rott (1994).This second method agrees more closely with the supervised classification done during the field trip.Thus, for this study, the ratio of Rott (1994) is considered optimal for application in this location.
In this work, the criterion of Andreassen et al. (2008), which defined glaciers with a surface area greater than 0.01 km 2 , was used.The uncertainty in the delimitation of margins and surfaces is determined by the image resolution, the precision of the co-registration (Ye et al., 2006) and the analyst's subjectivity, especially in the visual interpretation of false color image composites.In this case, since the interpretation of the MSS composites were supported by field observation, this error is considered insignificant.In the case of the other components, Ye et al. (2006) use an evaluation method based on Williams et al. (1997), Hall et al. (2003) and Silverio and Jaquet (2005).This method considers the spatial resolution and the co-registration error of each image to the base image or map from where the controls points were extracted.As mentioned before, in this study the base corresponds to the Landsat ETM+ scene.In all the images, the co-registration error was calculated at smaller than one pixel.Thus, a halfpixel size was assumed as the co-registration value for the application of the Eq.(1).
Following Ye et al. (2006), in the case of changes in the glacier's margin, the formula is expressed as: where, U T = Uncertainty in the measurement of the frontal position.λ = Spatial Resolution of each image considered on the evaluation.ε =The co-registration error of each image to the Landsat ETM+.
The measurements of changes in area are based on Hall et al. (2003), who recognize the importance of the resolution in this estimate.Ye et al. (2006) indicates that the co-register error between images is important, and consequently the uncertainty measurement is established as (Ye et al., 2006): where, U A = Uncertainty in the measurement of the change in area.This analysis showed that the greatest uncertainty in the frontal measurement is produced when comparing the margins resulting from the analysis of the MSS images (± 202 m).This also occurs with the estimation of the change in areas, where the maximum value corresponds to 5 hectares (Table 2).The present paper used these measurements as a criterion to select the periods in which the fronts and the areas register a significant change.When a comparison between frontal and areal changes showed a signal-to-noise ratio below 1, it was eliminated from the subsequent analysis.In applying this procedure, it is possible that the dates of frontal and areal changes recorded as significant do not match; this is the case for some glaciers in this study (Fig. 4).The implications of this procedure are discussed in Sect. 4.

Thickness changes
To estimate the changes in volume, map algebra was used to compare the altimetry of the regular cartography of the Instituto Geogr áfico Militar (IGM) with the DEM of the Shuttle Radar Topographic Mission (SRTM) project.

The IGM DEM
The contour line elevations and spot heights of regular IGM cartography were used to build a DEM (referred to here as the DEMIGM).These data present voids in the high zones (attributable to glacier accumulation zones) owing to an insufficient matching in the stereoscopic models due to the weak contrast for the snow.Consequently, these voids needed to be defined as Boolean masks to avoid artifacts in map algebra.Two interpolation methods were tested: the Inverse Distance Weight (IDW) and the Triangulated Irregular Network (TIN).Both procedures were applied in IDRISI ANDES (Eastman, 2006).The spatial resolution used was 30 m in order to obtain a comparison between the DEMIGM and the SRTM (see below).
While many studies do not incorporate an analysis of the interpolation method in order to determine if it is adequate for the study area (e.g.Schneider et al., 2007;Bown and Rivera, 2007), others have used qualitative (Cogley and Jung-Rothenh äusler, 2004) and quantitative (Rivera and Casassa, 2004) procedures to do so.Doing so appears necessary, since multiple studies that identify a source of error in the interpolation method (for details, see Fisher and Tate, 2006).
To select the appropriate interpolation method, a quantitative comparison was performed with the construction of DEMs using distinct quantities of source data and map algebra to compare them.
This comparison considered three stages.First, a testing area that had areas with and without ice was selected.Second, the sector was interpolated with both methods and a series of new DEMs were produced by randomly eliminating 10% of the original data from each new iteration.The quantity of iterations was established as n−1, where "n" is the number of pixels to eliminate in the testing area.Third, each new DEM was subtracted to the DEM built with 100% of the data, calculating the RMSE of each difference.Finally, the spatial variability of each interpolation method was analyzed by the application of a standardized Principal Component Analysis (Eastman and Fulk, 1993) to the series of iterations.This procedure seeks to improve the robustness of the interpolation method when using a procedure similar to jack-knifing (Rivera et al., 2007), while recognizing that the DEM produced can be one of several possibilities, as has been exemplified in studies in which the validation was performed with Monte Carlo simulations (Wechsler and Kroll, 2006).The sector analyzed was approximately 5 km 2 , corresponding to 5518 pixels.An "n" equivalent to 552 iterations was used (Figs. 2 and     6).As can be observed in Fig. 2a, the spatial distribution of the errors is observed to be more stochastic without an apparent pattern for the first component of interpolation iterations with IDW (36.1%), although there is a weak tendency towards a concentration of negative values to the west where the highest altitudes of the sampling area are located.This variability disappears in the second component (0.62%, not shown here), showing the values closest to zero in the area.In contrast, the first TIN component (55.7%, Fig. 2c) shows high variability located in the corners and borders.When comparing the difference for the original DEM graphs and the successive iterations (Fig. 2b and d), the IDW is observed to not surpass ± 2 m with a calculated RMSE of ± 0.53 m and TIN is observed to present differences greater than ± 20 m with a RMSE of ± 5.01 m.The distribution of differences between each iteration and the DEM original for IDW is adequately modeled by a Gaussian distribution (Kolmogorov-Smirnov test p > 0.1) with an asymmetry of −1.23 and kurtosis of 1.18.This result implies a slight tendency to interpolate values higher than the original ones, although with a leptokurtic normality.On the other hand, the TIN presents a massive underestimate of elevations, although they are impossible to model with a normal distribution (Kolmogorov-Smirnov test p < 0.01).
These results indicate that, even when there is reduced possibility of detecting a spatial pattern for the uncertainty of the interpolation method, the IDW does not present a bias towards the corners and borders.For the present type of comparison, this characteristic is key when there are void zones that can distort the like the ones that exist in the topographic map.Additionally, the IDW presents two more strengths: a RMSE almost 10 times smaller than the TIN and the possibility that the error can be The final calculation of uncertainty for this method is defined as the method applied by Bown and Rivera (2007), who utilize the formulation of Falkner (1995) and the RMSE of the interpolation.Using this method, the RMSE of DEMIGM corresponds to ± 17.2 m.

The SRTM DEM
The second model used to calculate the volumetric change corresponds to the DEM of the SRTM in its version of 1 arc second (approximately 30 m, referred hereafter as SRTM30).This DEM was obtained in format DTED 2 with absolute precision (nominal) of 23 m and 18 m for the horizontal and vertical, respectively (Advis and Andrade, 2007).
Recently, there has been a discussion on the validity of the results when using SRTM to calculate volumetric changes in glaciers.Indeed, Berthier et al. (2006) has documented a bias towards the underestimation of altitudes in mountain areas in the version of 90 m (the re-sampled version of the same SRTM30), which has resulted in the application of a correction in some subsequent studies (e.g.Berthier et al., 2007;M öller et al., 2007).However, Paul (2008) has questioned the existence of this bias in rough terrain where the combination of changes in slopes and curvatures could be an important factor for underestimation of elevations in coarse resolution DEMs.Thus, it is not possible to assign a unique bias in glacierized (normally smoother) and non-glacierized (normally rougher) terrain (M öller and Schneider, 2010).
The number of GPS points with geodetic quality in the study area is insufficient for local analysis since only one is located inside the limits of the national park (Fig. 1).As a result, the errors associated with the use of SRTM30 were quantified via three comparisons between the DEM and other data.The first comparison was performed by map algebra, comparing the non-glacier areas between the DEMIGM and SRTM30.The second comparison was between the non-glacier zones of the same sources 694 but above 2000 m, which corresponds to the lower limit of the glaciers inventoried by Rivera (1989).Finally, a larger comparison was made, in which the SRTM30 tiles which cover the administrative region was compared with the altimetry geodetic quality of 32 GPS points located within the same region (Figs. 1 and 4).In all the cases, a linear trend inversely proportional to the increase of altitude was observed.It means there is an underestimation of the SRTM30 when compared with the other data.The highest value corresponds to the differences between the maximum and minimum values of the SRTM-IGM pair (Table 3).In the case of the tandem SRTM-GPS, the values tended to be more negative above 1500 m, even when the comparison involves only three points (Fig. 3).This result is probably related to the null significance for the bias calculated with the linear fit.Above 2000 m, the differences are +7.8 and −9.9 m (RMSE = ± 4.31).
When both DEMs are compared, the bias derived from the linear fit tends to increase and becomes significant.The Durbin-Watson test indicates a possible serial correlation, which can indicate a certain robustness of this bias.However, since the R 2 is low, the variability of the data is not well modeled by this linear fit (Table 3).Additionally, it is notable that the RMSEs are similar and are twice the value calculated from the SRTM-GPS pair.Above 2000 m, two important altitude bins are observed.The first is located between 2250-2500 m, and the positive tendency is observed.The second, located above 2700 m, has a more pronounced tendency.According to these results, it can be concluded that there is a tendency of SRTM30 to underestimate altitude, although it tends towards a higher value than presented in previous studies at altitudes above 2000 m (see Berthier et al., 2006).
The impossibility of statistically validating the linear fit indicates that the significant bias does not represent a definitive tendency for the data.As a result, no correction was made.Instead, the highest RMSE value was preferred in order to constrain the uncertainty of the comparison.This implies that the uncertainty corresponds to ± 26.24 m.

Poeppig's chronicle
In 1835, German's explorer Eduard Poeppig published, in two volumes, the book "Reise in Chile, Peru und auf dem Amazonenstrome wahrend der jahre 1826-1829"; Keller (1960) translated volume I, which is completely dedicated to the exploration of Chile.The prologue of the translation (Keller, 1960:9) indicates that Poeppig initiated his exploration of Chile in the austral summer of 1827 after a voyage that began in the boreal autumn of 1826 in Baltimore (United States), arriving in Valparaíso (central Chile) 115 days later.Poeppig's access to the area of Sierra Velluda took place in spring 1828 when his notes indicate that he visited Concepci ón and Talcahuano on 30 October (Keller, 1960:347); both cities are part of the littoral landscape of the Bío Bío River basin.Sierra Velluda was accessed by the road that had been damaged in 1820 by the Lahar of the last eruption of the Antuco Volcano; he followed the course of the Laja River, calling the Volcano the "Silla Velluda" (Keller, 1960:375) because it looked like a riding saddle.The visit's record consists in a drawing that characterizes the landscape of the Trubunleo River valley, a southern tributary of the Laja River, located at 37 • 26 05 S and 71 • 27 03 W, with a NNW orientation.Here, Poeppig describes an oval-shaped valley with a spillway ending in a waterfall which drains into Laja River (Keller, 1960:377).Additionally, the saddle wall in the valley's interior is described as having a vertical section, which connects the valley's end with the glaciers (referred to by Poeppig as "ventisqueros"), with a relief equivalent to 3000 feet (900 m), which agrees quite well with the DEMIGM.
The drawing presented by Poeppig (Fig. 4a) was made "in the low mountains that close the valley (of the Trubunleo River) towards the north" (Keller, 1960:382), which means that it was drawn from the W side of the valley.The contemporary Fig. 4b seeks to emulate the position and field of view at the moment of the drawing.The glacier tongue that descends on the right side is notable, as described at the foot of the image (Keller, 1960:382).Because Keller was unable to exactly match Poeppig's viewpoint, the figure presented here gives a better understanding of the landscape features which permit the plotting of the glacier front in this year.The description and the drawn made by Poeppig allows the identification of the frontal position of glacier RC108371/2 in 1828 (Fig. 4c).From this position, a retreat rate of 3.5 m y −1 until 1961, was calculated.

Imagery analysis
Table 4 presents the frontal changes detected since 1961.However, the analysis of glaciers RC108376/1 and RC108376/2 begins in 1975 because they were not registered in the inventory of Rivera (1989).It is important to emphasize that the complex morphology of some glaciers (usually ice aprons and cirque glaciers) makes it difficult to extract a central line to measure length changes.In this sense, there is the potential for identifying negative length changes with positive areal changes.This was the situation in certain periods with some glaciers in the Sierra Velluda (Fig. 4).The implications of this are discussed in Sect. 4.
In general, the glaciers show a frontal retreat with only three exceptions.The first one corresponds to the section SW of the glacier RC108371/1 (referred to as Front 1) because it was difficult to discriminate the ice divides in this glacier.Indeed, this glacier can be morphologically classified as an ice apron (M üller et al., 1977), presenting three fronts of which Front 3 (the most northern face) did not present significant changes.The second exception is glacier RC108371/3, which has a NE aspect and presents a change equivalent to +1.9 m y −1 until year 2001.The third exception is the glacier RC108376/1, located on the W slope of the Sierra Velluda but with S aspect.This glacier registers the highest advance rate equivalent to +6.8 m y −1 until the year 2001.
The greatest retreat is observed in the majority of the glaciers of the SE slope, a pattern sustained since mid 1970.Glaciers RC108370/3 and RC108370/4 present the highest retreat, with practically twice the rates of the rest of the retreating glaciers.In 697 As can be appreciated from the graphs of frontal change (Fig. 5), two noticeable advance periods were detected in the time series.The first one corresponds to the 1970s, specifically to the period 1976-1979, when 64% of the glaciers show frontal advance.The second one occurs during the period 2001-2007, where 35% of the glaciers advanced.In addition, the situation of the glacier RC108370/8 is remarkable, since the only two periods in which significant change was observed (1961-1975 and 1975-2001) tend to compensate each other, making the overall change small.That is the reason why in Table 4 the complete change seems to be less than the error estimation.
These changes have modified the morphology of some glaciers.In effect, the change that glacier RC108371/2 has experienced is notable: it was first classified as a mountain glacier with a simple basin and a detached front, while it presently (position 2006) is much closer to a cirque glacier.On the other hand, the dynamics of the glaciers RC108370/7 and RC108370/8, which in 2006 did not exist or had a smaller dimension than the threshold used in this work, indicates that they were probably glacierets.
The changes in area for all the glaciers (Table 5) are generally for the period 1961-2007.Glacier RC108371/2, which presents values from 1828, together with glaciers RC108376/1 and RC108376/2 are an exception because they were not included in the 1961 inventory (Rivera, 1989).As in frontal changes, area changes are mainly positive in the 1970s and in some years between 2000 until 2007.Indeed, while 76.9% of the glaciers experience significant areal increase between 1975 and 1979, 84.6% increased between 2006 and 2007.Of the 13 glaciers, 9 have lost surface area in the period considered, although there is no spatial pattern (e.g.glacier aspect) associated to the retreat.The glaciers RC108370/3, RC108370/8, RC108376/1 and RC108376/2 increased their surface area in 19.6% (RC108376/1) to 66.7% (RC108376/2); all the advancing glaciers have S to SE aspect.
On the other hand, the largest absolute losses in surface area were observed in the ice apron RC108320/1, which with 3.35 km 2 have lost an area equivalent to 38.4% since 1961.With respect to their area in 1961, RC108370/4 (89.5%),RC108371/1 (69.4%) and RC108371/3 (51.4%) are the glaciers with the greatest losses.Indeed, glacier RC108371/2 has lost practically the same quantity of surface area (−0.53 km 2 ) since 1961 as it did in the period 1828-1961 (−0.61 km 2 ).In total, the Sierra Velluda has lost 43.8% of its surface area between 1961 and 2007.

Thickness change
The boundary of the glaciered area from 2007 was used as a Boolean mask for map algebra, extracting the areas without data from SRTM30 and those without stereoscopic vision of the regular IGM cartography.This means that not all the glaciers are equally represented in the calculation of the thickness change.Indeed, ice thickness changes in practically 90% of the glacier RC108370/2, 40% of RC108320/1 and 10% of RC108371/1 could not be computed.Additionally, since the values of change cannot be modeled by a normal distribution (Kolmogorov-Smirnov test p < 0.01), we did not eliminate extreme values.Considering these restrictions, a mean ice thickness change for the 13 glaciers equivalent to −20.06 ± 26.24 m with a rate of −0.51 ± 0.67 m/yr was registered.Even though this change is not significant, it is notable that 35% of the compared pixels have values superior to RMSE with an average equivalent to −38.63 m and a rate of −0.99 m/yr.If the change is observed within an altimetric range (Fig. 6b), a tendency toward negative values is observed as elevation increases.In effect, the linearization of the tendency found a rate of −0.011 m m −1 (p < 0.01) considering all the data, and −0.009 m m −1 (p < 0.05) considering only the values that are outside the range of RMSE.For the entire massif, there are two hypsometric concentrations of the changes.
The first and most important is observed around 2600 m, where the highest losses are found.The second one corresponds to altitudes higher than 3500 m.In the case of the positive values, the bin with the highest values, and coincidentally the only significant ones, is found between 2100 and 2300 m.The glaciers that show the most significant changes are principally located on the NE side of the massif, where recorded thinning rates were almost double the RMSE.On the other hand, it is notable that practically all the N and NW sides present insignificant, negative changes.Consequently, 92% of the glaciers thin with 46% of them presenting significant rates (Table 6).Additionally, the change of the glaciers RC108371/3 and RC108376/1 are noticeable, with more than 90% of the surface area experiencing significant thinning (Fig. 6a).

Discussion and conclusions
The Sierra Velluda glaciers have suffered a significant reduction in the last decades with up to 90% of area loss for some glaciers since 1961.These changes agree with the evolution that the Chilean glaciers have shown in the last few decades (Rivera et al., 2000) and also at a global level (Oerlemans, 2005).Zenteno et al. (2004) found losses of 78% for the glacier area of the Nevados de Chill án peaks (80 km north of Sierra Velluda) since 1862 with 46% of this loss registered since 1975.
This study has identified a maximum position of the glacier RC108371/2 for the year 1828.This position has certain concordance with those registered for the glacier Cipreses ( 34• S).Although with certain interpretation differences of the subsequent tendency, Le Quesne et al. ( 2009) and Araneda et al. (2009), agree that a maximum position of Cipreses Glacier was attained in 1842, associated to a frontal moraine.Araneda et al. (2009) attribute that position to the Little Ice Age.Espiz úa and Pitte ( 2009) register an advance of the Río Grande glaciers in Argentina ( 35• S) around 1830.On the other hand, Oerlemans (2005) suggests that glacier retreat on a global scale began around 1850.The evidence indicates that a wet period was produced between 1300 and 1800 as a result of persistent migration of the westerlies (Veit, 1996).Specifically, in the Laja Lagoon a cold period was detected between 1500 and 1900 by analyzing a sediment core for the last 2000 yr (Urrutia et al., 2010).Additionally, at 40 • S, a wet period between 1780 and 1820 was observed (B öes and Fagel, 2008).Consequently, even when this glacier can present a certain tendency of iceblock detachment, the maximum position defined here corresponds well with the other observations of central-southern Chile.
For the 20th Century, Le Quesne et al. ( 2009) found negative frontal changes in glaciers between 33 • S and 36 • S, together with a significant thinning determined at the glaciers Cipreses, Universidad and Palomo.South of the study area, Bown and Rivera (2007) identified significant thinning for the Casa Pangue glacier (41 • S) in the order of 85.1 m between 1961 and 1998.Even when Sierra Velluda has not registered values as high as Casa Pangue, the existence of glaciers with significant changes implies that the change process recorded in the last decades is a regional trend.
A notable fact of these changes registered is that produced during the 1970s, when more than 60% of the glaciers studied here registered frontal advance and surface area gains.Previous studies report advances in this decade (Rivera et al., 1997;Le Quesne et al., 2009;Masiokas et al., 2009).In fact, Le Quesne et al. (2009) found frontal advance for three of the five glaciers analyzed in this period.Coincidently, the Echaurren glacier (33 • S) registered a positive mass balance for the period 1977-1981 (Escobar et al., 1995).These positive fluctuations are concurrent with the record of precipitations of 1971 and 1972 for the meteorological station of the Laja Lagoon, which registered amounts greater than 3000 mm, indicating that these were the rainiest years of the decade (Mardones and Vargas, 2005) as observed in other stations in the region (Carrasco et al., 2005) and in the recent modeling (CONAMA, 2006).Additionally, these positive anomalies resulted in an increase in the lake's level in 1973 (Mardones and Vargas, 2005).This period is defined as a year of the warming phase of the El Ni ño Southern Oscillation (ENSO) due to the occurrence of two Sea Surface Temperature (SST) increase events around 2.5 • C: one between February and March and another between September and November (Reed, 1986).Since more than 62% of the lake's feed comes from runoff (Gonz ález et al., 2004), it is likely that a part of this fluctuation has been channeled through the accumulation of glaciers in the area.Indeed, if the formula of Oerlemans (2005) is applied to estimate the glaciers' reaction-times using the mean precipitation of 2000 m y −1 (Mardones, 2001) and the slopes of the glaciers derived from the inventory data of Rivera (1989), the response scale goes from one to five years for the smallest and largest, respectively.These reaction-times coincide with the generality of the behavior in the 1970s.
On the other hand, the consequent reduction of glaciers from the 1970s to the 2000s agrees with the shift of the ENSO (PDO) in 1976 since a negative correlation between precipitation post-1980 and the occurrence of the warm phase of the ENSO was found at 40  (Carrasco et al., 2008).However, this paper shows that there has been warming in the lower troposphere between 32

S to 41
• S, statistically significant for winter.Since winter is the rain season, a change in temperature can be changing the proportion of liquid water versus snow, decreasing accumulation.Indeed, these changes have been the base to estimate an increase in the isotherm of 0 • C and a consequent Equilibrium Line Altitude (ELA) rise of 127 m (Carrasco et al., 2005(Carrasco et al., , 2008)).The consequences of these changes are an increase of the ablation and a diminishment of the glacier accumulation area.
These results clearly show that even when the Sierra Velluda glaciers register behavior similar to that of most glaciers in central-southern de Chile, they are highly sensitive to high frequency climatic fluctuations and even to inter-annual variability.In this case, the most important inter-annual change in the climate elements is the in the precipitation registered.In this way, the Sierra Velluda glaciers have been more sensitive to changes in the precipitation than temperature in the last few decades.Indeed, the negative gradient of thickness change indicates that the most important changes are produced at the highest altitudes, which indicates that the accumulation areas are the ones being reduced.Considering that the models project a diminishment in Andean precipitation and an ascent in the 0 • C isotherm (CONAMA, 2006), these ice bodies are expected to reduce even more.
If we consider that a high correlation between the fluctuation of the Laja Lagoon level and the use of hydroelectric energy since 1970 as well as the descent of 27 m on the side and a low correlation between this level and precipitation (Mardones and Vargas, 2005), the region's water situation is quite fragile.This situation is important, especially since hydroelectric regulation is naturally limited by resource availability.Indeed, considering the results reported here, the limited correlation between the changes in lake level with precipitation can be reinterpreted as a growing tendency to use the lake reserves rather than an effective control.
In this paper, is remarkable the importance of estimating how significant the changes detected are, which is particularly important for mountain glaciers such as the ones analyzed.In contrast with larger glaciers, such as in the Patagonia, the high interannual variability and the relative small size of the glaciers require robust time series in order to analyze the changes when significant.In this case, the application of a significance assessment procedure gave more significant periods of change for area changes than for frontal changes (35% less records).It means that frontal changes are more sensible to image resolution.However, it must be considered that in this kind of glaciers can be more difficult to determine a central line for glacier length calculations.In this study, this difficulty can be an explanation for apparently anomalous behavior of frontal and areal changes for certain glaciers.Indeed, 9 glaciers showed one year with this discrepancy (Fig. 4).Nevertheless, since on average each glacier had 5 periods of significant frontal change and 8 of significant area change, this discrepancy is considered insignificant and the assessment method seems to be robust.
Since new change detection techniques based in images with better spatial resolution and methods that must be used to compare past measurements with future ones are required, these findings point out to the need for making significance analysis as Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

689
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | et al.

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper |

691
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

693
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | correctly modeled by a normal distribution.Consequently, IDW was selected as the interpolation method.

695
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 3 Results Paper | Discussion Paper | Discussion Paper | Discussion Paper | contrast, the lowest retreat values do not present an apparent spatial pattern since they occur in the glaciers RC108371/1 (N slope) and RC108370/8 (SE slope).

699
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper |

701
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

703
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | a previous step for monitoring purposes.Thus, an analysis of significance should be performed when comparing information from data sources with different resolutions.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .
Fig.1.Location map showing the study area inside its administrative region (Bío Bío).In cyan, the divide of the Bío Bío River Basin is highlighted.Red dots correspond to the GPS points used for the DEM validation (Fig.3b).

Table 1 .
Description of the satellite images used in this study.

Table 2 .
Rivera (1989)ed for each scene pair from the images used in this study.Aerial photo denotes the source of the information inRivera (1989).