60 Years of Glacier Elevation and Mass Changes in the Maipo River Basin, Central Andes of Chile

: Glaciers in the central Andes of Chile are fundamental freshwater sources for ecosystems and communities. Overall, glaciers in this region have shown continuous recession and down-wasting, but long-term glacier mass balance studies providing precise estimates of these changes are scarce. Here, we present the ﬁrst long-term (1955–2013 / 2015), region-speciﬁc glacier elevation and mass change estimates for the Maipo River Basin, from which the densely populated metropolitan region of Chile obtains most of its freshwater supply. We calculated glacier elevation and mass changes using historical topographic maps, Shuttle Radar Topography Mission (SRTM), TerraSAR-X add-on for Digital Elevation Measurements (TanDEM-X), and airborne Light Detection and Ranging (LiDAR) digital elevation models. The results indicated a mean regional glacier mass balance of − 0.12 ± 0.06 m w.e.a − 1 , with a total mass loss of 2.43 ± 0.26 Gt for the Maipo River Basin between 1955–2013. The most negative glacier mass balance was the Olivares sub-basin, with a mean value of − 0.29 ± 0.07 m w.e.a − 1 . We observed spatially heterogeneous glacier elevation and mass changes between 1955 and 2000, and more negative values between 2000 and 2013, with an acceleration in ice thinning rates starting in 2010, which coincides with the severe drought. Our results provide key information to improve glaciological and hydrological projections in a region where water resources are under pressure. Our results represent the first estimation of glacier mass changes for the second half of the 20th century in this region. Between 1955 and 2000, the highest thinning rates in the Maipo River Basin were observed in the Olivares sub-basin. possible for this is related to the location of this sub-basin, which is adjacent to two open-pit mines. mines generate that Olivares. leads


Introduction
Glacier mass balance is a key variable for understanding the response of glaciers to climate change and their contribution to sea level rise [1][2][3]. Traditionally, glacier mass balance has been obtained

Glacier Outlines
Two complete glacier inventories for the Maipo River basin are available. The earliest was compiled by Marangunic [34] for the year 1955. The second is the official Chilean glacier inventory based mostly on satellite imagery acquired between 2000 and 2003 [33]. In order to calculate glacier elevation and mass changes for the period 1955-2000, we used the glacier outlines from Maranguinc [34], which are mainly based on a stereoscopic analysis of HYCON aerial photographs from 1955 [35]. These aerial photographs were obtained as part of an agreement between Chile and the Inter-American Geodetic Survey of the United States. As a small fraction (3%) of the glacier area was not covered by the HYCON photogrammetric flight, these missing areas were completed by Marangunic [34] using the maps drawn by Lliboutry in 1956 [36].
For the period 2000-2013, we used the outlines from the Chilean glacier inventory [33]. The Chilean glacier inventory is mainly based on multispectral data from the Landsat satellite mission. However, some inconsistencies were detected over debris-covered areas. In order to improve the glacier outlines, we used a supervised semi-automatic method based on band ratio segmentation of TM4 and TM5 Landsat bands using the original Landsat [37]. Glacier elevation and mass change calculations were estimated over debris-free and debris-covered glaciers. Although this is a region with an abundant number of rock glaciers, these are not considered and should be thoroughly investigated.

Glacier Outlines
Two complete glacier inventories for the Maipo River basin are available. The earliest was compiled by Marangunic [34] for the year 1955. The second is the official Chilean glacier inventory based mostly on satellite imagery acquired between 2000 and 2003 [33]. In order to calculate glacier elevation and mass changes for the period 1955-2000, we used the glacier outlines from Maranguinc [34], which are mainly based on a stereoscopic analysis of HYCON aerial photographs from 1955 [35]. These aerial photographs were obtained as part of an agreement between Chile and the Inter-American Geodetic Survey of the United States. As a small fraction (3%) of the glacier area was not covered by the HYCON photogrammetric flight, these missing areas were completed by Marangunic [34] using the maps drawn by Lliboutry in 1956 [36].
For the period 2000-2013, we used the outlines from the Chilean glacier inventory [33]. The Chilean glacier inventory is mainly based on multispectral data from the Landsat satellite mission. However, some inconsistencies were detected over debris-covered areas. In order to improve the glacier outlines, we used a supervised semi-automatic method based on band ratio segmentation of TM4 and TM5 Landsat bands using the original Landsat [37]. Glacier elevation and mass change calculations were estimated over debris-free and debris-covered glaciers. Although this is a region with an abundant number of rock glaciers, these are not considered and should be thoroughly investigated.

Glacier Elevation and Mass Changes
Topography maps (1:50,000) of the Chilean central Andes were published by the Instituto Geografico Militar (IGM) in 1975, based upon HYCON (1955) aerial photographs. The original 1955 maps were horizontally referenced to the PSAD56 datum, UTM projection zone 19S. The vertical reference is based on orthometric height (m a.s.l.) obtained from precise trigonometric levelling with theodolites from tide gauge benchmarks located on the coast.
In this study, we used the 1955 topographic maps scanned by the IGM at 1200 DPI. We georeferenced these scanned maps to a common base using standard GIS procedures. We converted the original PSAD56 datum to WGS84 datum ( Figure 2) following the processing steps described by Farías-Barahona et al. [23]. This process consists of digitizing the 50 m contour lines plotted on the 1955 maps. Using kriging interpolation over these vectorised contours, a DEM was generated with a spatial resolution of 30 m.

Glacier Elevation and Mass Changes.
Topography maps (1:50,000) of the Chilean central Andes were published by the Instituto Geografico Militar (IGM) in 1975, based upon HYCON (1955) aerial photographs. The original 1955 maps were horizontally referenced to the PSAD56 datum, UTM projection zone 19S. The vertical reference is based on orthometric height (m a.s.l.) obtained from precise trigonometric levelling with theodolites from tide gauge benchmarks located on the coast.
In this study, we used the 1955 topographic maps scanned by the IGM at 1200 DPI. We georeferenced these scanned maps to a common base using standard GIS procedures. We converted the original PSAD56 datum to WGS84 datum ( Figure 2) following the processing steps described by Farías-Barahona et al. [23]. This process consists of digitizing the 50 m contour lines plotted on the 1955 maps. Using kriging interpolation over these vectorised contours, a DEM was generated with a spatial resolution of 30 m.  The SRTM-DEMs were generated using bi-static radar interferometry. These DEMs are the final products of SRTM, which was carried out in February 2000 on a near-global scale between 56 • S and 60 • N in both C-band and X-band radar frequencies by the National Aeronautics and Space Administration (NASA), and German Aerospace Centre (DLR) [38]. We used the non-gap filled SRTM C-band, available at the US Geological Survey, which is vertically referenced to the Earth Gravitational Model 1996 (EGM96). The non-gap filled DEM shows a very good coverage for the Maipo River Basin, except for a few small glaciers, which are not included in the analysis (e.g., Morado Glacier).
The TanDEM-X mission is a German bistatic satellite configuration jointly operated by the DLR and Airbus Defence and Space. It was launched in 2010, and since then, it has been acquiring single-pass InSAR image pairs [39]. The DLR provides InSAR image pairs as co-registered single-look complex images in HH polarization. We processed multiple along-track TanDEM-X scenes acquired on the same day, which are concatenated and processed using differential SAR interferometry (DInSAR) with SRTM C-band DEM as a reference (e.g., [4,40]) ( Figure 2). The topographic phase is simulated from SRTM DEM using the TanDEM-X interferometric parameter to obtain a differential interferogram. This differential approach is generally applied to reduce difficulties in subsequent phase unwrapping, especially in high-mountain regions [41].
To remove the noise from the differential interferogram, a Goldstein filter was applied [42]. Areas with low coherence (<0.2) were masked out [41] before unwrapping the phase. The Minimum Cost Flow (MCF) algorithm was applied for phase unwrapping [43], keeping the reference point on off-glacier flat terrain. Phase jumps and residual phase ramps were manually checked and removed from the unwrapped phase, which is finally converted into absolute differential heights. These differential heights were added to the topographic heights from the SRTM-C DEM to generate the TanDEM-X DEM. The resulting TanDEM-X DEM was geocoded with the SRTM-C DEM to keep planimetric consistency.
In order to ensure three-dimensional (3D) alignment between SRTM and TanDEM-X DEMs, we further corrected for any horizontal and vertical biases before mosaicking the TanDEM-X DEMs for our study area. This is based on an off-glacier mask derived from a slope threshold of 15 • on SRTM DEM (close to the mean slope of the glacierised area) and Landsat-based Normalised Difference Vegetation Index (NDVI) and Normalised Difference Water Index (NDWI) masks. These masks were used to eliminate areas of dense vegetation and water bodies from the analysis [4]. The stable mask was used during the co-registration process of the TanDEM-X and SRTM-C DEMs, and it was performed according to Nuth and Kääb [44]. The main objective is that the difference (dh) between two DEMs on stable terrain should be as close as possible to zero. A bilinear plane fit was used to remove any remaining vertical offset, where the elevation difference is a function of projected x-and y-coordinates in the mapped DEM on stable ground [45]. Finally, the corrected DEMs were mosaicked [4] (Figure 2).
In order to derive the geodetic mass balances, we subtracted the difference of the topographic DEM map and STRM DEM after horizontally co-registering them in order to estimate elevation change for the period 1955-2000 and between the SRTM DEM and TanDEM-X for the period 2000-2013. For the 1955-2000 period, we discarded unrealistic positive values in the retreated part of the glacier present in the period 2000-2013. We removed strongly deviating values by applying a quantile filter (1%-99%) (e.g., [46]) and filled them in using the hypsometry function (e.g., [4]). Since the SRTM and TanDEM-X DEMs were acquired at the end of the austral summer (February-March), we converted the elevation changes to mass change using a density conversion factor of 850 ± 60 kg m −3 [47]. For the specific elevation and mass change, we used the mentioned inventories (Section 3.1) as a reference and not a temporal mean, as suggested by [1] due to the lack of third inventory.  (Table 1 and Figure 2). Results for four of the fifteen glaciers were presented and discussed by Burger et al. [20] and Farías-Barahona et al. [23].

LiDAR Measurements
LiDAR is a surveying method consisting of an active system, which measures the time taken for a laser signal to travel from the source to a ground surface object and its return to the sensor. It allows the creation of a three-dimensional set of elevation points when the platform location is known. During the four field campaigns, two different types of LiDAR equipment were used: Trimble Harrier 68i and Terra Remote Sensing Mark II [48][49][50]. Both airborne system devices were equipped with Inertial Measurement Unit (IMU) and GNSS (Global Navigation Satellite System). Both sensor systems allow multiple laser pulses and echoes in the air, increasing the number of points per area. During the field measurements, a self-calibration method was implemented in order to minimise systematic errors, which consist of a GNSS kinematic track parallel to the airborne flight.
The LiDAR survey achieved a nominal point density between 2 to 4 points per m 2 . In Table 1, we present the LiDAR coverage for the measured glaciers. Due to cost and time constraints, some glaciers were surveyed along their centre lines instead of their entire area; hence, for some glaciers, the LiDAR coverage reaches only~30% of the total glacier area ( Table 1). In order to estimate the elevation changes of the measured glaciers, we interpolated the LiDAR's point clouds using a Triangulated Irregular Network (TIN) to obtain DEMs of 1 m spatial resolution. The reference system is WGS84 (UTM projection 19S) and ellipsoidal height. Finally, we calculated the elevation changes using the local hypsometry method with 100 m elevation bands.

Uncertainties and Error Assessment
To calculate the uncertainties of our geodetic mass balance (dM), we followed the procedure described in Braun et al. [4] and Malz et al. [45], who use Equation (1), which includes: (1) uncertainties in the DEM differencing (δ dh dt ) (including spatial autocorrelation), (2) error from the glacier outlines (δ A ), (3) uncertainty from the volume-to-mass conversion using a fixed density (δ ρ ), and (4) uncertainty from radar signal penetration (V pen ).
where (δ dh dt ) corresponds to the relative elevation differences ( dh dt total) between SRTM-C and TanDEM-X estimated over areas of stable ground. (A) and (δ A ) correspond to the glacier area and the uncertainty in extent mapping, respectively. For the manual digitisation of glacier outlines, an uncertainty of 5% Remote Sens. 2020, 12, 1658 7 of 19 of the total area was assumed, based on Paul et al. [51]. To calculate mass conversion uncertainties, we considered the uncertainty of the volume to mass conversion (δ ρ ) as an extra 7% of the elevation change uncertainty, which corresponds to ±60 kg m −3 [47].
To obtain (δ dh dt ), we used standard deviations (σ dh/dt ) derived from 5 • slope bins on stable ground and calculated the slope-binned area-weighted standard deviation (σ dh dt AW ) on glacierised areas (e.g., 4, 40, 45) (Equation (2)). Additionally, we calculated the spatial autocorrelation to obtain the error budget [52]. We generated semi-variograms of the elevation change values on stable ground to determine the mean lag distance (dc) to determine the correlation area S e = dc 2 *π, and S f corresponds to the glacier area.
Other sources of uncertainties come from the radar signal penetration of SRTM-C and TanDEM-X. Radar signal penetration depends on several factors, such as the type of surface (ice/snow), water content, and structure of the medium, alongside the radar frequency. Signal penetration was anticipated to be high in the Northern Hemisphere as SRTM data was acquired in February (mid-winter), and therefore, must be accounted for (e.g., [41]). In the hyper-humid region of Patagonia (Northern Patagonia Icefield), studies reported negligible radar signal penetration into the firn/snow below 2900 m a.s.l. (e.g., [53]), and only a few centimetres of band penetration have been found between the SRTM-X and SRTM-C bands (e.g., [54]). From the evaluation of 20 glaciers located between 3500 and~5000 m a.s.l, the Maipo River Basin also shows no significant band penetration between X and C bands [23]. Photographs taken four days before the TanDEM-X acquisition show the glacier surface comprising bare ice~4600 m a.s.l (Supplementary Figure S1). However, for consistency with other studies regarding uncertainty assessment using TanDEM-X over the Andes [4,40], we applied a similar procedure as Braun et al. [4], who assumed a linear increase of radar signal penetration (V pen ) difference from 0 m up to 5 m at the maximum elevation of the study region, which is integrated over the area above the Equilibrium Line Altitude (ELA) to obtain (V pen ), in this case,~4000 m a.s.l [15].
We identified other sources of uncertainties as: (1) The classification of very small glaciers in the first glacier inventory, (2) potential bias in the vertical reference (m a.s.l) between the 1955 topographic maps and SRTM-C, and (3) LiDAR accuracy. There are some small glaciers in the second inventory that were not mapped in the first one , this led to uncertainties related to whether the first inventory interpreted these as snow patches or glacierets (<0.05 km 2 ). In any case, these were not considered in the glacier mass balance, because they represent a very small fraction of the total area (~1%). Although the topographic maps have been previously used in the Chilean Andes (e.g., [55][56][57]), another potential bias can be caused by different vertical reference used in the 1955 topographic maps and the SRTM DEMs. The elevation of the IGM topographic maps is referenced to orthometric heights (m a.s.l.) obtained from trigonometric levelling, whereas SRTM-C is converted to elevation (m a.s.l) with the EGM96 system. However, the co-registration results, in which a good agreement between DEMs of each sub-basin in stable areas was observed (±0.05 m a −1 ), may also cover any geoid undulation residual.
The LiDAR measurements did not cover areas of stable ground. Therefore, we assumed the reported accuracy corresponding to ±0.30 m between GNSS and LiDAR. This accuracy is based on simultaneous GNSS and LiDAR airborne measurements on stable ground and glacier surfaces during the field campaigns [49,50].

Glacier Elevation and Mass Changes from Satellite Products
Spatially distributed glacier elevation changes in the Maipo River Basin for the 1955-2000 and 2000-2013 periods are shown in Figure 3. The regional average glacier elevation change rates for the Maipo River Basin are −0.13 ± 0.05 m a −1 and −0.18 ± 0.08 m a −1 for the 1955-2000 and 2000-2013 periods, respectively. While the first period evidenced spatially heterogeneous patterns between sub-basins (Table 2), we observed a consistent mass loss during our second period, except for the Colorado sub-basin. We computed a total mass loss of 2.43 ± 0.26 Gt for glaciers in the Maipo River Basin between 1955 and 2013. Table 2 shows the mean geodetic mass balance rates for each sub-basin of the Maipo River Basin. Glaciers in the Olivares sub-basin show the most negative glacier mass balance rates with −0.29 ± 0.07 m w.e.a −1 between 1955 and 2013 (Figure 4), and the Olivares Alfa Glacier shows the highest negative glacier mass balance of −0.62 ± 0.08 m w.e.a −1 (1955 and 2013). Figures 5 and 6 show the correlation between rates of elevation change and latitude with topographic and glaciological parameters (aspect, elevation, and area) for glaciers with areas larger than 0.05 km 2 in both study periods. In the first interval, the largest negative rates of elevation change correspond to glaciers with easterly (E) and northwest (NW) aspects, whereas glaciers with south (S) and southwest aspects (SW) show positive elevation changes ( Figure 5), specifically in the latitude between 33.5 • and 34 • S. However, in glaciers located at~33 • S, high rates were observed (Olivares sub-basin) independent of their aspect ( Figure 5). During this period, the thinning rates were concentrated in glaciers with an area larger than 1 km 2 ( Figure 5 and Supplementary Figure S2).
In the second interval (between 2000 and 2013), many glaciers have changed their aspect due to their fragmentation. Overall, we observed the highest percentage of thinning rates in the glaciers with E aspect. While a reduction in the thinning rates is observed on glaciers with W and NW aspects, this is related to an increase of the mean elevation of the glaciers (Figure 6), especially for glaciers located at 33.5 • S. Unlike the first period, between 2000 and 2013, we observed an increase in the thinning rates in glaciers with S and SW aspects ( Figure 6). During this period, the thinning rates were concentrated in the largest glaciers (>5 km 2 ) ( Figure 6 and Supplementary Figure S2).
It is interesting that some of the sub-basins exhibited opposite temporal behaviour. While negative mass balance results in Colorado were obtained in the first period, a slight positive or neutral glacier mass balance was observed between 2000 and 2013. The opposite behaviour was found in the Volcan and the Upper Maipo sub-basins: positive/neutral mass balance between 1955 and 2000 and negative mass balance between 2000 and 2013. Our observations also revealed that during the second analysis period, the debris-covered glaciers have also been thinning, in some cases at similar rates as the debris-free glaciers (Supplementary Figures S3 and S4). We also confirmed some surge events inferred from positive elevation changes in the ablation area of some glaciers during both periods (Supplementary Figure S5).

Glacier Elevation Changes with LiDAR
Strong acceleration in the thinning rates after 2010 were detected for the glaciers measured with airborne LiDAR. The elevation change maps show high thinning rates for most glaciers (Supplementary Figure S6). Overall, we observed the highest negative glacier elevation changes in the Olivares and Molina sub-basins. For glaciers of the Olivares sub-basin (Juncal Sur, Olivares Alfa, Olivares Beta, and Olivares Gama) an average elevation change of -1.51 ± 0.10 m a -1 was found between 2011 and 2015. Both glaciers located at the Molina sub-basin, the Paloma, and Rincon glaciers, together showed an average of -1.81 ± 0.14 m a -1 between 2012 and 2015.

Glacier Elevation Changes with LiDAR
Strong acceleration in the thinning rates after 2010 were detected for the glaciers measured with airborne LiDAR. The elevation change maps show high thinning rates for most glaciers (Supplementary Figure S6). Overall, we observed the highest negative glacier elevation changes in the Olivares and In Figure 7, we show a comparison of elevation change rates for all observation periods for the fifteen glaciers measured by LiDAR in comparison with the InSAR and topographic maps dataset. Overall, the differences between glaciers are related to mean altitude and debris-covered factors. In terms of individual glaciers, the most negative glacier elevation change was the San Francisco Glacier located in the Volcan sub-basin, with −2.12 m a −1 on average (2009)(2010)(2011)(2012)(2013)(2014)(2015). This glacier presents the lowest mean altitude (~3300 m a.s.l.) of the glaciers measured with airborne LiDAR.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 20 In Figure 7, we show a comparison of elevation change rates for all observation periods for the fifteen glaciers measured by LiDAR in comparison with the InSAR and topographic maps dataset. Overall, the differences between glaciers are related to mean altitude and debris-covered factors. In terms of individual glaciers, the most negative glacier elevation change was the San Francisco Glacier located in the Volcan sub-basin, with -2.12 m a -1 on average (2009)(2010)(2011)(2012)(2013)(2014)(2015). This glacier presents the lowest mean altitude (⁓3300 m a.s.l.) of the glaciers measured with airborne LiDAR.

Discussion
Our results revealed that negative values of glacier mass balance dominated between 1955 and 2013, but with some spatial and temporal differences in the responses of sub-basins and glaciers. While glacier mass balance showed negative and positive values in the different sub-basins and glaciers for the first study period , there was a generalised negative trend in the second one (2000-2013), with the exception of the Colorado sub-basin. After 2010, we observed that mass losses accelerate everywhere, with nearly all analysed glaciers showing high negative mass balance. The differences between these periods can be explained by a combination of precipitation variability and a sustained increase of air temperature associated with global warming. Glacier mass balance in the central Andes is very sensitive to precipitation, which is mainly related to El Niño Southern Oscillation (ENSO) events, where positive mass balances have been associated with El Niño events and vice-versa with La Niña events [21,22,23,58]. However, at longer time scales, the Pacific Decadal Oscillation (PDO) and the Interdecadal Pacific Oscillation (IPO) play an important role in the climatic conditions over the Central Andes [59,60]. It is noteworthy that at present, this region is being affected by a severe drought that started in 2010, with an annual rainfall deficit that ranges between 25% and 45% [19,25]. Besides the precipitation deficit, recent studies have indicated an increase in the autumn and spring temperature of inland Maipo [61]. Both factors may affect glaciers at a local and regional spatial scale, but also depend on the response time of the individual glaciers to those forcings, latitude, and the glacier characteristics, such as elevation and aspect, and the variability of ice melt [20].

Discussion
Our results revealed that negative values of glacier mass balance dominated between 1955 and 2013, but with some spatial and temporal differences in the responses of sub-basins and glaciers. While glacier mass balance showed negative and positive values in the different sub-basins and glaciers for the first study period , there was a generalised negative trend in the second one (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013), with the exception of the Colorado sub-basin. After 2010, we observed that mass losses accelerate everywhere, with nearly all analysed glaciers showing high negative mass balance. The differences between these periods can be explained by a combination of precipitation variability and a sustained increase of air temperature associated with global warming. Glacier mass balance in the central Andes is very sensitive to precipitation, which is mainly related to El Niño Southern Oscillation (ENSO) events, where positive mass balances have been associated with El Niño events and vice-versa with La Niña events [21][22][23]58]. However, at longer time scales, the Pacific Decadal Oscillation (PDO) and the Interdecadal Pacific Oscillation (IPO) play an important role in the climatic conditions over the Central Andes [59,60]. It is noteworthy that at present, this region is being affected by a severe drought that started in 2010, with an annual rainfall deficit that ranges between 25% and 45% [19,25]. Besides the precipitation deficit, recent studies have indicated an increase in the autumn and spring temperature of inland Maipo [61]. Both factors may affect glaciers at a local and regional spatial scale, but also depend on the response time of the individual glaciers to those forcings, latitude, and the glacier characteristics, such as elevation and aspect, and the variability of ice melt [20].   Our results represent the first estimation of glacier mass changes for the second half of the 20th century in this region. Between 1955 and 2000, the highest thinning rates in the Maipo River Basin were observed in the Olivares sub-basin. One possible reason for this is related to the location of this sub-basin, which is adjacent to two open-pit mines. These mines generate dust that can be transported and deposited over the surface of glaciers in Olivares. This process potentially leads to a lowering of the surface albedo, which in turn increases the absorbed incoming shortwave radiation and energy Our results represent the first estimation of glacier mass changes for the second half of the 20th century in this region. Between 1955 and 2000, the highest thinning rates in the Maipo River Basin were observed in the Olivares sub-basin. One possible reason for this is related to the location of this sub-basin, which is adjacent to two open-pit mines. These mines generate dust that can be transported and deposited over the surface of glaciers in Olivares. This process potentially leads to a lowering of the surface albedo, which in turn increases the absorbed incoming shortwave radiation and energy available for melt. It is beyond the scope of this work to assess and quantify the impact of the deposition of dust originated from mining activity; hence, additional studies are required. Another possible explanation for these differences between sub-basins is that the Maipo River Basin presents a shift between two hydroclimatic regimes, where an increment in the precipitation amounts towards south of the 34 • S have been observed [60]. This could also explain the observed north to south glacier mass balance gradient observed during this period ( Table 2 and Figure 5), where positive and nearly neutral mass balance was observed for Volcan and Upper Maipo for the period 1955-2000, respectively ( Figure 5). The difference between both sub-basins is related to elevation of the glaciers, where the Volcan sub-basin contains big glaciers that are located at higher elevations on the western and southern flanks of active volcanoes (Marmolejo and San Jose Volcano).
Throughout the second period (2000-2013), our results agree well with recent findings in the central Andes of Chile and Argentina [4,10,24], including those focusing on individual glaciers (e.g., [20][21][22][23]). Again, in the Olivares sub-basin, we obtained the highest negative glacier mass balances during this period. However, in this period, the latitudinal gradient was not observed, instead, we obtained similar negative glacier mass balance, with the exception of glaciers located in the Colorado sub-basin, where it seems that these glaciers are in balance with the present climate, likely due to their high elevation. This was also confirmed by the LiDAR measurements.
Comparing our results with those for the Southern Andes macro-region (RGI 17) [5,62] can be misleading, since this macro-region encompasses the Patagonian Icefields which provide nearly all of the regional mass loss [4,10,62]. In Figure 8, we illustrate the glacier elevation changes of these previous studies that include glaciers of the Maipo River Basin from 2000 onwards. On a local scale, our estimate of glacier mass balance for the Maipo River Basin was more positive than that derived by the glaciological method for the Echaurren Norte Glacier. This glacier shows largely negative mass balance rates [21][22][23], but the representativeness of the Echaurren Norte Glacier is still a matter of further analysis, since there are differences between the glaciological and geodetic methods. The geodetic mass balance of the Echaurren Norte Glacier [23] was positive between 2000 and 2009, and this trend has been confirmed by other authors [10,20]. The acceleration thinning rates based on LiDAR data indicates the role of the extensive drought that has been affecting the region since 2010 [19,25]. It is likely that after 2013, the drought exacerbated thinning rates, which may also explain the differences between the two recent region-wide glacier mass balance records for the central Andes [4,10]. In fact, observing other studies, it seems that the effects of the drought were not prominent before 2013 for glaciers located in the Yeso sub-basin [20,23,63,64]. For instance, GNSS measurements on Bello and Yeso glaciers showed −0.24 and −0.17 m a −1 between April 2012 and March 2014 [63], and measurements during the ablation period of 2012 and 2013 showed 0.67 and 0.44 m, respectively [64]. The acceleration of thinning rates has also been observed by other authors [10,20,65]. For instance, Hernandez et al. [65] estimated thinning rates of −1.64 ± 0.08 m a −1 and −0.94 ± 0.10 m a −1 for the Olivares Alfa and Olivares Beta glaciers between 2013 and 2016 respectively, using airborne LiDAR. Our results also showed similar rates: Olivares Alfa Glacier: −1.61 ± 0.08 m a −1 and Olivares Beta Glacier: −1.39 ± 0.08 m a −1 . The slight difference may be related to the time span and the dataset coverage (Table 1 and Supplementary Figure S6). Conversely, glaciers located in the Colorado sub-basin seem to be in equilibrium after an important reduction of their area size (Figure 3), as we observed in our elevation change results (2000-2015) (Figure 7 and Table 2). This may be related to the new mean altitude of these glaciers, which in the case of Yeso 1 and 2, are nearly located at Remote Sens. 2020, 12, 1658 14 of 19 5000 m a.s.l., whereas in a prior period  both glaciers were part of a larger, lower-lying glacier ( Figure 3).
Our results also showed similar rates: Olivares Alfa Glacier: -1.61 ± 0.08 m a -1 and Olivares Beta Glacier: -1.39 ± 0.08 m a -1 . The slight difference may be related to the time span and the dataset coverage (Table 1 and Supplementary Figure S6). Conversely, glaciers located in the Colorado subbasin seem to be in equilibrium after an important reduction of their area size (Figure 3), as we observed in our elevation change results (2000-2015) (Figure 7 and Table 2). This may be related to the new mean altitude of these glaciers, which in the case of Yeso 1 and 2, are nearly located at 5000 m a.s.l., whereas in a prior period  both glaciers were part of a larger, lower-lying glacier ( Figure 3).  We also calculated thinning rates for several debris-covered glaciers during the second study period and we estimated that some of them are thinning at similar rates (Supplementary Figure S2). Ayala et al. [26] showed that the glacier runoff contribution of the low-lying, debris-covered Piramide Glacier is similar to that of the high-elevation, debris-free Bello and Yeso of the same sub-basin, suggesting that elevation differences and debris cover can explain the similar mass balance of these two types of glaciers [20]. In fact, when comparing the Echaurren Norte Glacier and Piramide Glacier (debris-covered), which present a similar mean elevation (~3700 m a.s.l.), thinning rates at the Echaurren Norte Glacier are twice of those at the Piramide Glacier ( Figure 7). Hence, although the Piramide Glacier presents negative rates, the supraglacial debris clearly reduces the thinning compared with debris-free surfaces at the same elevation. Nevertheless, it is likely that ice cliffs and supraglacial lakes are playing a key role by enhancing the total ablation [66,67]. Ice cliffs on debris-covered glaciers have been identified as "hot-spots" and hence as major contributors of mass loss, which outplay the insulating role of debris cover. This phenomenon was confirmed by the LiDAR measurements, which show strong ablation at several ice cliffs and supraglacial lakes with depths of up to 40 m on the Piramide Glacier (Supplementary Figures S3, S4, and S6). Other debris-covered glaciers in the region were also displaying ice cliffs and supraglacial lakes in the SRTM and TanDEM-X elevation changes (Supplementary Figure S4), but further monitoring studies with high-resolution data (e.g., LiDAR) are required to understand the role of ice cliffs in the mass balance of debris-covered glaciers at the basin or region scale in the central Andes.
We observed signals of surge in some glaciers, although not with the same magnitude as in other regions, e.g., Himalaya [68]. Overall, only a few surging glaciers have been the focus of historically detailed studies in the central Andes, such as the Horcones Inferior [69,70] and Grande del Nevado [71,72] in Argentina, and Cachapoal Glacier in Chile [73]. Other surge events were reported in the past by Lliboutry [35], such as those of Juncal Sur (1947), Nieves Negras (1927), and Glaciar del Río (1935). However, we did not observe clear surge events in those glaciers between 1955 and 2000, with the exception of Nieves Negras Glacier. Instead, there were small signals of surge at Picos del Barroso, Loma Larga, Sierra Bella, Azufre (Tupungatito volcanic complex), and Oeste del Cerro Alto Glaciers (Supplementary Figure S5). For the second period, our results also confirmed small signals of surge-type behaviour for a few glacier systems, as noted by Falaschi et al. [74], in glaciers such as Loma Larga, Oeste del Cerro Alto, and Sierra Bella.

Conclusions
In this study, we have estimated the long-term region-wide glacier elevation and mass changes for the Maipo River Basin in the central Chilean Andes, which was previously sparsely investigated. We estimated a glacier mass balance of −0.12 ± 0.06 m w.e.a −1 with a total mass loss of 2.43 ± 0.26 Gt for the glaciers throughout the Maipo River Basin between 1955 and 2013. Heterogeneous spatial patterns in the elevation changes were observed between sub-basins in the first period , and an increase in the thinning rates in the following analysed period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013). A strong acceleration in the thinning rates was found until 2015, where fifteen glaciers were observed using in-situ techniques. The acceleration of glacier thinning after 2010 is coincident with the current severe dry period (2010 onwards) in the central Andes of Chile and Argentina.
Over the entire study period, the most negative glacier mass balance was found in glaciers located in the Olivares sub-basin. This large mass loss may be attributed to a decreasing surface albedo and less precipitation amounts in comparison with the southern sub-basins. However, additional studies are required to elucidate the relationship between dust on the glacier surface and anthropogenic activities (e.g., mining and transport of city pollution to mountain areas).
Our results provide the first long-term spatio-temporal glacier elevation and mass change analysis, and again prove the feasibility of using remote sensing techniques to monitor glaciers. However, it is critical to extend the field glacier monitoring over the basin, due the relevance of glaciers in the Maipo River Basin. Here, we also provide a key dataset to further validate hydrological projections as well as for the implementation of water management plans, as the demand for water resources has been considerably increasing during the last decades.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/10/1658/s1: Figure S1: Melt conditions and bare ice upper sections of (a) Bello Glacier and (b) Yeso Glacier at 4600 m a.s.l. Photographs taken by Geoestudios© on 27 March 2013, which is approximately 4 days after the TanDEM-X acquisition. Hence, we discard any radar signal penetration at this altitude. Figure S2: Basic glaciological parameters plot. Elevation change in percentage (%) per number of glaciers in relation to size interval (km 2 ). * corresponds to very small glaciers (<0.05 km 2 ) in the first glacier inventory , which were not taken into account. Figure

Funding:
We thank the financial support by the German Research Foundation DFG under grant BR 2105/14-1, the German Space Agency (DLR), and the Ministry of Economy (BMWi) under contract 50EE1544, and the support by the Chilean National Science Foundation (Fondecyt) under grant agreement Nº 1161130 "Present and future projections of glacier extent and water yield at Maipo basin in Central Chile". The University of Erlangen-Nuremberg supported this study from its Open Access Publishing fund.