Mass balance and area changes of glaciers in the Cordillera Real and Tres Cruces, Bolivia, between 2000 and 2016

Abstract

Climate change has led to a significant shrinkage of glaciers in the Tropical Andes during the last decades. Recent multi-temporal quantifications of ice mass loss at mountain range to regional scale are missing. However, this is fundamental information for future water resource planning and glacier change projections. In this study, we measure temporally consistent glacier area changes and geodetic mass balances throughout the Bolivian Cordillera Real and Tres Cruces based on multi-sensor remote-sensing data. By analyzing multi-spectral satellite images and interferometric SAR data, a glacier recession of 81 ± 18 km² (29%; 5.1 ± 1.1 km² a⁻¹), a geodetic mass balance of −403 ± 78 kg m⁻² a⁻¹ and a total ice mass loss of 1.8 ± 0.5 Gt is derived for 2000–2016. In the period 2013–2016, ice mass loss was 21% above the average rate. A retreat rate of 15 ± 5 km² a⁻¹ and a mass budget of −487 ± 349 kg m⁻² a⁻¹ are found in this more recent period. These higher change rates can be attributed to the strong El Niño event in 2015/16. The analyses of individual glacier changes and topographic variables confirmed the dependency of the mass budget and glacier recession on glacier aspect and median elevation.

1. Introduction

Glaciers and ice caps outside of the polar ice sheets are defined as essential climate variables by the Global Climate Observing System (GCOS) and key indicators for climate change by the Intergovernmental Panel on Climate Change (IPCC, 2014). In particular, glaciers in low latitudes, also known as tropical glaciers, are very sensitive to changes in climatic conditions (e.g. Kaser and Osmaston, 2002; Rabatel and others, 2006, 2013). Small changes in temperature and precipitation strongly influence the mass balance of tropical glaciers since they are close to melting conditions throughout the year (Francou and others, 2003; Vuille and others, 2008).

The El Niño Southern Oscillation (ENSO) is known to affect the precipitation pattern in the Tropical Andes. The decadal precipitation variability appears larger than the long-term trend, which is additionally regionally contrasted. Indeed, in Ecuador and northern Peru, a long-term increase in precipitation since the mid-20th century is documented (Vuille and others, 2003). In contrast, a decreasing trend is reported for the Bolivian Altiplano region (Vuille and others, 2003; Haylock and others, 2006). Moreover, between 1939 and 2006, Vuille and others (2008) documented a near-surface air temperature increase of 0.68 °C for the Tropical Andes (1°N to 23°S), which agrees with local observations by other studies (e.g. Vuille and Bradley, 2000; Mark and Seltzer, 2005). Gilbert and others (2010) analyzed englacial temperature measurements in a 138 m deep borehole close to the summit of Illimani (6340 m a.s.l., Bolivia). The authors derived a general warming trend of 1.1 ± 0.2 °C for the 20th century, with a higher rate in recent years (4.3 ± 1.4 °C per century for 1985–1999). As a consequence of the changes in climatic conditions, the equilibrium line altitude (ELA) of the glaciers in this region rose (Rabatel and others, 2013; Réveillet and others, 2015). At the Quelccaya Ice Cap, Peru, an ELA shift from 5275 to 5414 m a.s.l. between 1975 and 2012 was observed by Hanshaw and Bookhagen (2014), and at Zongo Glacier, Bolivia, a temperature sensitivity of the ELA of 150 ± 30 m °C⁻¹ was estimated (Lejeune, 2009).

In the Tropical Andes, glaciers and ice caps are an important water resource and regulate the seasonal runoff, particularly in central and southern Peru and Bolivia (Vuille and others, 2018). No seasonal snow cover exists outside the glaciers, which could act as a buffer for the dry season runoff (Lejeune and others, 2007). During the dry season, from April until September, the glacier meltwater significantly contributes to the drinking water supply of Bolivia's capital La Paz (e.g. Soruco and others, 2015). Additionally, glacier runoff is highly important for irrigation, hydropower generation and water consumption of local communities (Vuille and others, 2018). The extreme drought in Bolivia in 2016 forced the administration to declare an emergency and to ration water (Schoolmeester and others, 2018). The authors also reported a very high contribution of 61% of the glacier meltwater to the available water supply in La Paz during this extreme event. Buytaert and others (2017) estimated an even higher monthly maximum of 85.7% during drought years. On a multi-decadal timescale (1963–2006), Soruco and others (2015) reported an average glacier meltwater contribution to the water supply of La Paz city, the administrative capital of Bolivia, and the neighboring city El Alto (together ~2.3 million inhabitants) of 27% during the dry season. The glacier recession in the Tropical Andes leads to temporarily increased meltwater runoff. However, on long-term scales, the runoff will decrease as glaciers retreat further (Vuille and others, 2018). Estimates by Huss and Hock (2018) indicate that the Bolivian glaciers have already reached or are close to ‘peak water’, the point in time with maximum meltwater runoff. Soruco and others (2015) projected a reduction in runoff of 24% during the dry season if the glaciers disappear in the La Paz drainage basins.

The glacier area in the Bolivian Cordillera Oriental (Cordillera Apolobamba, Cordillera Real and Tres Cruces) reduced by 43% between 1986 and 2014 (Cook and others, 2016) and even a higher retreat of ~50% was found in the small mountain range of Tres Cruces between 1975 and 2009 (Albert and others, 2014). For the southern wet outer tropics, where the Bolivian Cordillera Oriental is located, a mass change rate of −0.44 ± 0.11 Gt a⁻¹ (850 kg m⁻³ density scenario) was estimated by using radar remote-sensing data for the period 2000–2013 (Braun and others, 2019). This estimate also includes largely glacierized areas in Peru and the authors do not provide information on detailed spatial scales since the study focused on a continent-wide mass-balance analysis. Detailed studies of glacier changes in Bolivia have been carried out at the Zongo, Chacaltaya and Charquini glaciers (e.g. Rabatel and others, 2006; Soruco and others, 2009b; Réveillet and others, 2015; Soruco and others, 2015). Volume changes between 1963 and 2006 of 21 glaciers in the Cordillera Real were measured by Soruco and others (2009a) based on photogrammetric measurements. The authors estimated for nearly one-third of the Cordillera Real's glaciers (376 glaciers) an area loss of 48% and a volume loss of 43% throughout the study period using an empirical relationship between mass balance, exposure and average altitude based on the 21 monitored glaciers.

At Cordillera Real and Tres Cruces, no region-wide, spatially detailed and multi-temporal information on the contemporaneous glacier area and mass changes is available. Therefore, the scope of this paper is to monitor glacier evolution throughout the Bolivian Cordillera Real and Tres Cruces for the period 2000–2016 in order to obtain enhanced information on recent glacier changes. Our primary objectives are to:

• • Generate detailed glacier inventories for the years 2000, 2013 and 2016, which are temporally consistent to the mass-balance analyses

• • Measure regional and glacier-wide geodetic mass balances for the periods 2000–2013, 2000–2016 and 2013–2016

• • Analyze area change and mass-balance results in comparison with topographic and climatic variables

The time intervals are defined by the availability of the used remote-sensing data and complement the observation periods of other studies.

2. Study site

The Cordillera Real and the smaller mountain range Tres Cruces are home of the majority of the Bolivian glaciers. In agreement with Jordan (1991), the glacierized regions of Mururata and Illimani were included in the Cordillera Real (Fig. 1). The analyzed mountain ranges extend in the NW-SE direction and have a length of ~200 km. They separate the wet Amazon Basin from the arid Altiplano Plateau. The study site is located in the outer tropics region as defined by Troll (1941), particularly in the southern wet outer tropics according to Sagredo and Lowell (2012). Glacier surface albedo and the overall radiation budget, mainly controlled by cloud cover and solid precipitation, strongly force the surface mass balance of these glaciers (Wagnon and others, 1999; Sicart and others, 2005). During the wet summer season (January–April), accumulation, as well as snow and ice melt rates, are the highest, whereas in the dry winter season (May–August), less melting occurs and sublimation becomes important. The transition period (September–December) is characterized by sporadic precipitation events and high melt rates, since the solar radiation increases (Sicart and others, 2011; Rabatel and others, 2012). ENSO events have also been identified as an important climate driver of the interannual variability of the glacier surface mass balance (e.g. Francou and others, 2003; Vuille and others, 2008; Rabatel and others, 2013). El Niño events amplify the mass loss, whereas during La Niña events, lower mass loss or even positive mass balances are found (Wagnon and others, 2001; Vuille and others, 2008). Based on the Randolph Glacier Inventory (RGI) 6.0 (RGI Consortium, 2017), the glacierized area in the study region ranges between ~3700 (too low, caused by an artifact in RGI 6.0, see Section 5.1) and 6400 m a.s.l., with an extension of ~329 km² (measured in UTM 18S, since the majority of the glaciers in RGI Region 16 are located in UTM 18S zone), corresponding to ~63% (according to RGI 6.0) of the glacierized area in Bolivia, consisting of ice caps, and mountain and valley glaciers. At three glaciers in the study region, continuous surface mass and energy-balance monitoring are (were) carried out (Zongo Glacier since 1991, Chacaltaya Glacier 1991–2008 and Charquini Sur Glacier since 2002) within the framework of a Franco-Bolivian collaboration through the French Service National d'Observation GLACIOCLIM. Zongo Glacier is a benchmark glacier of the World Glacier Monitoring Service (WGMS) and a reference Cryonet Station of the Global Cryosphere Watch led by the World Meteorological Organization.

Fig. 1. (a) Map of Bolivia. Dashed square indicates the study area, blue polygons show the glacier extent based on the RGI 6.0. Background: Natural Earth (b) Topographic map of analyzed mountain ranges and TanDEM-X (TDX) and Pléiades data coverage. Background: SRTM © NASA.

3. Data and methods

3.1 Area changes

A detailed comparison of the RGI outlines with high-resolution satellite imagery of the respective years (from Google Earth) revealed a misclassification of certain areas. As stated in the Technical Report of the RGI 6.0 (RGI Consortium, 2017), a more rigorous examination and delineation of the glacier outlines might reduce the mapped glacier area of RGI Region 16 by up to 20%. This discrepancy can bias area change computation when comparing the RGI outlines with other outlines.

Consequently, glacier outlines are mapped, using cloud-free LANDSAT images (Table S1, provided by the United States Geological Survey, USGS), for 2000, 2013 and 2016. To identify glacierized regions, the NDSI of top of atmosphere reflectance values and the band ratio (BR) (GLIMS Algorithm Working Group) of red and shortwave infrared bands of the digital number values are calculated. A combination of NDSI and BR (for deep shadows) information is used to identify the glacier outlines. By means of high-resolution satellite imagery (Google Earth) from the respective years, the threshold values for the NDSI of 0.8 and the BR of 1.7 for LANDSAT 5 TM data and 1.5 for LANDSAT 8 OLI data are carefully selected. The resulting binary masks are transferred into polygons and cropped into individual glacier basins using the ice flowshed delineations of the RGI. Afterward, the glacier inventories are then visually inspected and misclassified areas, such as proglacial lakes and snow patches, are removed by means of high-resolution satellite data (Google Earth). The surface area (S) of the final inventory and the individual glaciers is calculated in UTM 19S projection. Additionally, the topographic parameters (minimum, maximum and median elevation, slope and aspect) of the inventoried glaciers are calculated based on the void fill SRTM Version 3 digital elevation model (DEM) (NASA JPL, 2013, see also next section) according to RGI6.0 Technical Report (RGI Consortium, 2017).

3.2 Elevation changes and mass balance

In this study, multiple DEMs derived from Interferometric Synthetic Aperture Radar (InSAR) data are used to determine glacier surface elevation changes throughout 2000–2016.

The Shuttle Radar Topography Mission (SRTM) of the National Aeronautics and Space Administration (NASA) and the German Aerospace Center (DLR) acquired bi-static C-band InSAR data of landmasses between 60°N and 54°S from 11 to 22 February 2000. The goal of SRTM was to facilitate the generation of a nearly global, high-resolution and consistent DEM. Several DEM products were derived from SRTM data. In this study, the void-filled LP DAAC NASA Version 3 SRTM DEM (1 arcsec ground resolution) (NASA JPL, 2013), reprojected to UTM projection, is used. Further DEMs are generated using data sets from the TanDEM-X (TDX) mission of the DLR. The TDX mission has been acquiring bi-static X-band InSAR data since 2010 (Zink and others, 2011). A nearly complete coverage of the study site by TDX data is available for 2013 and 2016. In order to reduce the biases on the elevation change calculation caused by seasonal surface elevation changes and SAR signal penetration (see Section 5.3), TDX data within the same season as the SRTM are selected for 2013. Unfortunately, in 2016, TDX data are only available for the period September to November. The potential bias caused by this seasonal offset is discussed in Section 5.3. The footprints of the TDX acquisitions used are shown in Figure 1 and the main characteristics are given in Table S1.

The TDX DEMs are processed according to Seehaus and others (2016) and Malz and others (2018). Along track acquisitions from the same path and date are concatenated and a differential interferogram is generated using the SRTM DEM as a reference DEM. Subsequently, the differential interferogram is filtered and unwrapped using either the branch cut or minimum cost flow algorithm (the best results are manually selected) and the differential phase is converted into differential height values. Finally, the reference DEM heights are added and the resulting DEM is geocoded. In the next step, the TDX DEMs are bi-linearly vertically co-registered to the SRTM DEM on stable ground. Stable points are defined as areas with a slope of <25° (using the SRTM DEM) and by masking out vegetation (NDVI threshold of 0.3) and water (NDWI threshold of 0.1) bodies by means of LANDSAT data. Afterwards, the TDX DEMs are horizontally co-registered to the SRTM DEM using the iterative algorithm of Nuth and Kääb (2011). To remove the remaining vertical offsets, the TDX DEMs are again bi-linearly vertically adjusted to the SRTM elevations and finally mosaicked to one single DEM, including a time stamp for the individual pixels.

The SRTM DEM and the resulting TDX DEM mosaics are subtracted to estimate the elevation change rates Δh/Δt for the respective periods. As a time stamp for the SRTM elevation data, the mean date of the SRTM 2000-02-16 is taken. Voids in the SRTM DEM were filled using elevation information from different sources (NASA JPL, 2013) without providing the corresponding date information. Thus, the provided masks used are rejected the void-filling data in our Δh/Δt analyses. Since, the glaciers in Bolivia show in general a retreat pattern (Cook and others, 2016), the glacier outlines of the starting year of the respective periods are applied to measure Δh/Δt and consequently ice volume losses on the glacierized areas. Additionally, pixels with slopes >50° (3.9% of the glacierized area) are not considered in the further analyses since DEMs are less reliable on steep slopes (Toutin, 2002) and nearly no ice is aggregated there (based on field observations). Gaps in the Δh/Δt data on glaciers are filled by hypsometric extrapolation at 100 m elevation bands. In order to remove outliers, for each elevation band, the Δh/Δt values are filtered by applying three times the normalized median absolute deviation (NMAD) (Brun and others, 2017) and the resulting filtered mean Δh/Δt values are applied for the individual elevation bands. The void-filled SRTM DEM is used as the elevation reference for the hypsometric analyses.

Information on the regional glacier mass budget is revealed by calculating the geodetic mass balance (ΔM/Δt). This is defined by Fountain and others (1997) as the elevation change rate integrated over the glacierized region multiplied by an average density ρ, also called volume to mass conversion factor. Two density scenarios were applied. For the first one, an average density of 850 kg m⁻³ is considered, as suggested by Huss (2013) for alpine glaciers. The second scenario uses two densities of 600 and 900 kg m⁻³ for areas above and below the ELA (5144 m a.s.l., Rabatel and others, 2012), respectively (e.g. Gardelle and others, 2012; Kääb and others, 2012). Kääb and others (2012) suggested that the second density scenario is more suitable for glaciers where the volume changes are mainly forced by ablation and accumulation and not by glacier dynamics. Therefore, we used the results of the second scenario for the further analyses and discussion and kept the results of the first scenario for comparison with other studies, which used a similar density approach. Huss (2013) showed that the density can vary significantly. According to these findings, an uncertainty on the density of ±60 kg m⁻³ is applied for mass-balance quantification over periods longer than 10 years and mass changes larger than ±200 kg m⁻² a⁻¹. Otherwise, an uncertainty of ±300 kg m⁻³ is taken (see also Braun and others, 2019), that is, for the period 2013–2016.

The geodetic mass balance of the individual glaciers is computed according to Cogley and others (2011), and the voids in the elevation change fields are handled following the recommended method by McNabb and others (2019). For each glacier with <60% voids in the Δh/Δt fields, the hypsometric distribution of Δh/Δt is computed. Elevation intervals of 50 m are applied for glaciers with an elevation range >500 m, otherwise 10% of the glacier's elevation range is used for the binning. Information for elevation bands without Δh/Δt data is obtained by fitting a third-order polynomial to the hypsometric Δh/Δt distribution. Only glaciers with a data coverage of at least two-thirds of the elevation bands are considered in the analysis. The resulting mass budgets of the individual glaciers are analyzed in combination with topographic attributes (area, median elevation, aspect) in order to detect correlations. Mass balances of individual glaciers are computed by applying the second density scenario. Artifact in the Δh/Δt fields can bias the mass-balance results of individual glaciers. Therefore, a three times NMAD filter is applied in the hypsometric analysis (see regional analysis above) and finally a 2–98% quantile filter is used to discard possible remaining outliers.

3.3 Uncertainty analysis

The accuracy of the area measurements (δ _S) is estimated according to Malz and others (2018). The authors used the accuracy estimation of 3% for outlines of alpine glaciers, mapped by means of LANDSAT imagery from Paul and others (2013), and included a scaling factor for the different perimeter–area ratios.

The uncertainty in the geodetic mass-balance estimates (δ _ΔM/Δt) is calculated by:

(1)$$\hskip-5pt \scale 88% { \delta _{\Delta {\rm M}/\Delta {\rm t}} = \sqrt {{\left( {\displaystyle{{\Delta M} \over {\Delta t}}} \right)}^2\left( {{\left( {\displaystyle{{\delta_{\Delta {\rm h}/\Delta {\rm t}}} \over {\Delta h/\Delta t}}} \right)}^2 + {\left( {\displaystyle{{\delta_{\rm S}} \over S}} \right)}^2 + {\left( {\displaystyle{{\delta_\rho} \over \rho}} \right)}^2} \right) + {\left( {\displaystyle{{V_{{\rm pen}}} \over {\Delta t}}{^\ast}\rho} \right)}^2}. } $$

It consists of the error contributions of the elevation change measurements (δ _Δh/Δt), the glacier inventory (δ _S), the ice density (δ _ρ) and an uncertainty due to the differences of the SAR signal penetration (V _pen/Δt). The accuracy of the averaged Δh/Δt measurements (δ _Δh/Δt) is affected by the relative vertical accuracy over the survey glacier area (δ _{Δh/Δt(vert)}) and the uncertainty caused by the hypsometric extrapolation. The relative vertical error δ _{Δh/Δt(vert)} is evaluated by analyzing Δh/Δt on the ice-free areas and masking out water bodies and vegetation (see Section 3.2). Strongly deviation values are removed by applying a 2–98% quantile filter. Since the offsets depend on the slope (Fig. 2; Gardelle and others, 2012; Vijay and Braun, 2016), the area weighted standard deviations (σ _AW) based on the slope distribution (5° bins) are calculated. Additionally, the deviations for each slope bin are filtered by applying 3*NMAD before computing σ _AW. To estimate the accuracy of the glacier-wide averaged elevation change (δ _{Δh/Δt(vert)}), we followed the approach of Rolstad and others (2009) (applied by, e.g. Fischer and others, 2015; Brun and others, 2017) in order to take into account spatial out correlation. For glacier areas (A _gl) larger than the correlation area (A _cor = π × d²_cor; d_cor: correlation length) δ _{Δh/Δt(vert)} is calculated as:

(2)$$\delta _{{\rm \Delta h/\Delta t}\lpar {{\rm vert}} \rpar } = \sqrt {\displaystyle{{A_{{\rm cor}}} \over {5A_{{\rm gl}}}}} \sigma _{{\rm AW}}.$$

Fig. 2. Area on (light blue) and off glacier (red) and off glacier dh/dt distribution (blue dots) based on slope. Error bars indicate NMAD of dh/dt values for each slope bin. Dotted line shows the slope threshold used in the analysis. Note: the glacier area is scaled by a factor of 10.

For A _gl<A _cor as:

(3)$$\delta _{{\rm \Delta h/\Delta t}\lpar {{\rm vert}} \rpar } = \sigma _{{\rm AW}}.$$

Since we are integrating the elevation change over large regions for the regional analysis, the above approach is applied for each contiguous glacierized area and not for individual glaciers. For the glacier-specific mass balances, the individual glacier areas are used. The region-wide δ _{Δh/Δt(vert)} is calculated as the area weighted average of all individual glacierized areas. The spatial auto-correlation of our off-glacier elevation change measurements is analyzed by generating semi-variograms with 1 00 000 random samples, binned in 30 m distance intervals for a maximum distance of 20 km. The correlation range is estimated by applying a spherical semi-variogram function and a mean d _cor of 394 m is found for our data sets (an example is illustrated in Fig. S1).

In order to account for the error due to the hypsometric extrapolation, we followed the approach of Berthier and others (2014) and multiplied the error δ _{Δh/Δt(vert)} by a constant factor (we chose a factor of 2, Braun and others, 2019) for the unmapped areas to finally retrieve δ _Δh/Δt.

SAR signal penetration is evaluated by comparing a DEM derived from a Pléiades acquisition on 28 May 2013 (see Fig. 1) and the TDX DEM in 2013 (TDX data covering the Pléiades DEM is from 31 January 2013). Both DEMs are coregistered (see Section 3.2) and the hypsometric distribution of the elevation differences in glacier areas (inventory from 2013) are analyzed (Fig. S2). A negligible influence on the elevation change measurements by the SAR signal penetration is found (see Section 5.2 for a detailed discussion). Consequently, no correction of the SAR data is applied. However, we considered a potential influence in the total error quantification of the mass balances, since the surface conditions in 2000 and 2016 might have been different as in 2013 and SRTM data were acquired at a lower SAR frequency.

The uncertainty contribution due to an offset in the SAR signal penetration between SRTM and TanDEM-X is estimated according to Malz and others (2018). Deviations in the elevation change rates caused by differences in the SAR signal penetration of C- and X-band data in snow and ice surfaces in areas below the ELA (estimated at an average elevation of 5144 m a.s.l. for Zongo Glacier by Rabatel and others (2012), and assumed to be representative of the regional average in the current study) are neglected, since no penetration is assumed in the typical snow and ice melt affected surfaces. Above the ELA, a linear increase toward 5 m of penetration at the maximum elevation of the study region is applied to estimate V _pen. The respective density of the applied scenario is used to convert V _pen into mass. Moreover, for the second scenario, the potential offset in the mass budgets caused by differences in the SAR signal penetration is ~40% lower due to the lower density considered above the ELA.

4. Results

4.1 Area changes

The glacier outlines obtained for 2000, 2013 and 2016 are plotted in Figure 3 (detailed subsets are plotted in Fig. S3). A reduction of the extent of the glacierized areas is clearly visible. This shrunk from 281 ± 14 km² in 2000 to 246 ± 14 km² in 2013 (−12%) and 200 ± 11 km² in 2016 (−29%) amounting to an area loss of 81 ± 18 km². The area change values of the study area are summarized in Table 1. The average shrinkage rates show a stronger recession toward recent years from −2.7 ± 1.5 to −15 ± 5 km² a⁻¹ in the periods 2000–2013 and 2013–2016, respectively. Moreover, the number of glaciers >0.01 km² reduced from 332 in 2000 to 272 glaciers in 2016 and 34 glaciers (nine glaciers in 2000–2013) completely vanished. Additionally, glacier fragmentation was observed at several glaciers (see Fig. S3). In the Cordillera Real, the glacier area shrunk from 244 km² in 2000 to 175 km² in 2016 (−28%) and in Tres Cruces from 37 km² in 2000 to 24 km² in 2016 (−34%). The links between the area changes of the individual glaciers (in %) and topographic variables, such as area, median elevation and aspect, are illustrated in Figure 4 for the period 2000–2013 (Figs S4 and S5 for the periods 2000–2016 and 2013–2016). Highest retreat rates are generally found for low-lying and small glaciers. Area change measurements for small glaciers are more sensitive to mapping errors due to the typically higher perimeter–area ratio. This partially leads to a wider range and even positive values in relative area changes for smaller glaciers as compared to the larger ones.

Fig. 3. Glacier area changes. Ice divides (black polygons) are from the glacier outlines in 2000. Dashed polygons indicate subsets illustrated in Figure S3 (Supplement) Background: SRTM hillshade © NASA 2000.

Fig. 4. Relative area changes (d_S) (2000–2013) of individual glaciers (dot color) plotted against glacier size (dot size), median elevation (distance from center) and mean aspect (orientation). Red circle: equilibrium line altitude (ELA) from Rabatel and others (2012).

Table 1. Observed glacier changes for different time intervals

Δt: observation period; ΔS: glacier area change; ΔS/Δt: glacier area change rate; Δh _M/Δt: average measured surface lowering rate on S _M; Δh _E/Δt: average extrapolated surface lowering rate on S; S: analyzed glacier area at the beginning of the observation period. Note: slopes >50° are excluded; S _M: fraction of S covered by Δh _M/Δt measurements; ΔM/Δt: mass balance (average and specific); ΔM: total mass change in the observation period.

^a Volume to mass conversion factor of 850 kg m⁻³.

^b Volume to mass conversion factor of 900 and 600 kg m⁻³ for areas below and above ELA, respectively.

4.2 Elevation changes and mass balance

Figure 5 shows the measured (unfiltered) surface elevation change rates for the periods 2000–2013 and 2013–2016 (detailed subsets are presented in Figs S6 and S7). A clear surface lowering at most of the glacier termini is found. Gaps in the Δh/Δt fields are due to incomplete SRTM data coverage (coverage of 67% of the glacierized area) and voids in the TDX DEMs (coverage of 94 and 92% of the glacierized areas in 2013 and 2016, respectively, within the TDX footprints) caused by typical SAR limitations like layover and shadowing (see also Sections 3.2 and 5.2). In the period 2000–2016, the elevation change is measured only on 50% of the glacier area since no complete coverage of the study site by TDX data is available in 2016 (Fig. 1). The hypsometric distribution of the glacier area and the obtained Δh/Δt (filtered) are illustrated in Figure 6 for the period 2000–2013 (Figs S8 and S9 for the periods 2000–2016 and 2013–2016, respectively). Most of the surface lowering is observed at elevations below the ELA. The measured surface elevation gains in the upper glacier areas can be neglected considering the small fraction of coverage (<1%). An increase in the extrapolated mean surface elevation change from −0.51 ± 0.11 m a⁻¹ for 2000–2013 to −0.74 ± 0.34 m a⁻¹ for 2013–2016 is found. All mean Δh/Δt estimates are presented in Table 1. Consequently, a more negative geodetic mass balance of −0.12 ± 0.09 Gt a⁻¹ is revealed for the period 2013–2016 as compared to −0.10 ± 0.03 Gt a⁻¹ in the period 2000–2013. In total, −1.8 ± 0.5 Gt of ice was lost in the observation period.

Fig. 5. Glacier surface elevation changes (unfiltered; left panel: 2000–2013; right panel: 2013–2016). Dashed polygons indicate subsets illustrated in Figures S6 and S7 (Supplement). Background SRTM hillshade © NASA 2000.

Fig. 6. Hypsometric distribution of glacier area (light blue), glacier area with dh/dt measurements (red) and mean, filtered dh/dt values (blue dots) of each hypsometric bin for the observation period 2000–2013. Error bars indicate NMAD of dh/dt for each hypsometric bin.

Glacier-wide mass balances for ~60% of all glaciers could be quantified for the study periods 2000–2013 and 2013–2016. Due to voids in the SRTM DEM and incomplete coverage by the TDX DEM in 2016, individual glacier mass balances of only 40% of all glaciers are obtained in the period 2000–2016. Histograms indicating the distribution of the glacier-specific mass balances, as well as minimum and maximum values, are provided in Figure S10. Figure 7 illustrates links between the mass balances obtained for each glacier and topographic variables for the period 2000–2013 (Figs S11 and S12 for the periods 2000–2016 and 2013–2016, respectively). A trend of more negative mass balances for glaciers with lower median elevation as well as for glaciers with northward orientation is revealed.

Fig. 7. Specific mass balance (spMB) (2000–2013) of individual glaciers (dot color) plotted against glacier size (dot size), median elevation (distance from center) and mean aspect (orientation). Red circle: equilibrium line altitude (ELA) from Rabatel and others (2012). Note: only glaciers with >40% elevation change data coverage, which is spread over >2/3 of the hypsometric distribution are included.

5. Discussion

5.1 Area changes

Our quantified glacier area loss of up to 30% (2000–2016) in the study region is in line with previous observations in this region and highlights the dramatic glacier changes in Bolivia. Several studies (e.g. Rabatel and others, 2006; Soruco and others, 2009b; Réveillet and others, 2015) have reported glacier recession of individual glaciers in the Cordillera Real. Soruco and others (2009a) reported in the central part of the Cordillera Real an area change of −48% (−1.5% a⁻¹) between 1975 and 2006 and −1.7% a⁻¹ between 1997 and 2006, which is comparable to our observed retreat rate in the Cordillera Real of −28% (−1.8% a⁻¹) in the period 2000–2016. In Tres Cruces, Albert and others (2014), Cook and others (2016) and Veettil and others (2018) measured glacier area changes of ~−55% (−1.6% a⁻¹) for 1975–2009, −47.3% (−1.7% a⁻¹) for 1986–2014 and −35.1% (−1.7% a⁻¹) for 1995–2016, respectively. We observed a higher retreat rate of −2.1% a⁻¹ for the period 2000–2016 in Tres Cruces, which agrees with the reported increased glacier retreat in Tres Cruces after 2010 by Cook and others (2016). Cook and others (2016) mapped area changes throughout the Bolivian Cordillera Oriental for the first time and reported recessions of 43.1% (1.5% a⁻¹) for the period 1986–2014, which is similar to the observation of 1.4% a⁻¹ by Veettil and others (2018) for 1995–2016. Cook and others (2016) found high retreat rates of 8.1 km² a⁻¹ (2.5% a⁻¹) in the period 1999–2010 and lower rates of 4.0 km² a⁻¹ (1.2% a⁻¹) in the period 2010–2014, including the Cordillera Apolobamba. Taking into account the different sizes of the studied areas (351 km² in 2014, including the Cordillera Apolobamba), this leads to an average loss rate of 1.4% a⁻¹ for 1999–2014, based on the results of Cook and others (2016). This value is comparable to our observed retreat rates of 0.9% a⁻¹ and 1.8% a⁻¹ in the periods 2000–2013 and 2000–2016, respectively, considering the shift in the observation intervals and a bias by the higher retreat rates in the Cordillera Apolobamba as compared to the Cordillera Real (observed by Cook and others (2016) and Veettil and others (2018)).

In the time interval 2013–2016, we obtained an enhancement of the retreat rates toward 15 ± 5 km² a⁻¹ in our study area. A similar area change trend was observed at Peruvian glaciers by Seehaus and others (2019). The increased rates are comparable to the finding of Cook and others (2016) in the period 1986–1992 (14.5 km² a⁻¹). These high retreat rates are most likely to be associated with the strong El Niño events that occurred during these two observation periods (see Section 5.3 and Seehaus and others (2019) for more details).

Our mapped glacier area for 2000 is similar to that measured by Cook and others (2016). It is ~13.3% smaller than the glacier extents based on the RGI6.0 (2000–2003), which was expected according to the limitation indicated in the Technical Report of the RGI 6.0 (see Section 3.1). Moreover, the minimum glacier altitude of 3706 m a.s.l. according to the RGI 6.0 is caused by a misclassified glacier area. We find a more reasonable minimum glacier altitude of 4418 m a.s.l. Figures 4, S4 and S5 indicate that the highest relative area changes are in general found for small and low-lying glaciers, confirming the findings by Veettil and others (2018). This can be partly attributed to the fact that small glaciers are generally located at lower elevations and thus experience more ablation due to changes in the climatic conditions (see Section 1). Some of the small and low-lying glaciers even lost their accumulation areas due to the shift of the ELA (Vuille and others, 2008). Rabatel and others (2013) reported the disappearance of the monitored Chacaltaya Glacier in the Cordillera Real in 2009. Its extinction was already projected by Ramirez and others (2001) since the estimated mean ELA was above its upper extent. The authors proposed that Chacaltaya Glacier is representative of the Bolivian glaciers, of which ~80% are smaller than 0.5 km² (based on our glacier inventories, 63 and 76% were smaller than 0.5 km² in 2000 and 2016, respectively) with average altitudes close to or below the ELA. All of our identified extinct glaciers (2000–2016) had a median elevation below 5400 m a.s.l. and were smaller than 0.5 km², which supports the predictions of Ramirez and others (2001). At Cerro Charquini (5392 m a.s.l.), Rabatel and others (2006) derived an upward shift of the ELA by 160 m between the LIA maximum (~1660 AD) and 1997 and attributed it to a temperature increase by ~0.6 °C as the main driving factor. For Zongo Glacier, Rabatel and others (2012) estimated an average ELA of 5144 ± 67 m a.s.l. for the period 1991–2006 and Vuille and others (2018) computed an ELA shift toward 5400 and 5700 m a.s.l. by 2100 based on the CIMP5 RCP4.5 and RCP8.5 scenarios, respectively. Consequently, further disappearance of glaciers can be expected considering the large number of glaciers at low altitudes with area losses of more than 60% (2000–2016) and the projected upward shift of the ELA.

The higher shrinkage rates at Tres Cruces, as also observed by Veettil and others (2018), can be related to the observed link between median glacier elevations and shrinkage as well. The median glacier elevation in the study region for the 2000 inventory was 5342 m a.s.l. (5365 m a.s.l. for Cordillera Real only), whereas at Tres Cruces the median glacier elevation was 5285 m a.s.l. Moreover, the uppermost elevation at Tres Cruces is 5728 m a.s.l., clearly below the maximum elevation at Cordillera Real of 6425 m a.s.l. Thus, the glaciers at Tres Cruces are even more strongly affected by the changing climate conditions and the associated uplift of the ELA than the glaciers in Cordillera Real.

Cook and others (2016) have evaluated the formation of proglacial lakes and the glacier lake outburst flood (GLOF) risks in Bolivia from 1986 to 2014. The continuous large-scale glacier shrinkage revealed in this region until 2016 might have led to further extension or formation of new proglacial lakes. Moreover, further rapid glacier shrinkage is expected (e.g. Ramirez and others, 2001; Vuille and others, 2018). Thus, we strongly support the suggestion by Cook and others (2016) to further monitor proglacial lakes and GLOF risks in Bolivia.

5.2 Elevation changes and mass balances

An average glacier surface elevation change of −9.9 m (−0.60 m a⁻¹) and highest surface lowering of up to ~50 m (~3.1 m a⁻¹) at the lower glacier reaches are derived from SAR remote-sensing data for the study period 2000–2016. The estimated accuracy of the elevation change (δ _Δh/Δt) quantified over stable terrain of 0.88–1.04 m (long-term average: 0.08 m a⁻¹ for 2000–2016; short-term average: 0.24 m a⁻¹ for 2013–2016) is comparable to or even better than the reported values of similar studies in other mountain regions (e.g. 0.22–0.28 m a⁻¹ for 2000–2012 by Vijay and Braun (2018), 0.43 m a⁻¹ for 2000–2013 by Vijay and Braun (2016), 0.22 m a⁻¹ for 1999–2008 by Gardelle and others (2012) and 0.21–0.28 m a⁻¹ for 2000–2015/16 by Malz and others (2018)). The evaluation of the elevation differences on glacier areas between the TDX DEM and Pléiades DEM in 2013 revealed an average offset at regions below and above the ELA of −0.7 and −0.1 m, respectively (Fig. S2). The total offset due to SAR signal penetration is small considering the temporal difference of 4 months between both acquisitions (end of wet season and beginning of dry season, accumulation in the upper reaches and high melt rates toward the terminus), the average surface-lowering rates (2000–2016) of −1.0 and −0.7 m a⁻¹ at areas below and above the ELA, respectively, and the uncertainty of the elevation change measurements. The SAR data (SRTM and TDX) were acquired during the transition and wet season (Sections 2 and 3.2), when the melt rates are highest and meltwater is most likely present in the snow pack, strongly limiting the signal penetrations. Moreover, we compared only DEMs based on SAR acquisitions, thus only relative differences in the SAR signal penetration affect the volume change estimation. Therefore, we conclude, that the impact of the SAR signal penetration on elevation change computations is negligible. However, in order to be consistent with previous studies and to account for potential differences in the glacier surface conditions in 2000 and 2016 as well as for the differences in the SAR frequency, a contribution to the error budget of the mass balances is considered. The estimated uncertainty V _pen, due to the difference in the SAR signal penetration between SRTM and TanDEM-X (C-band vs X-band), amounts to 0.288 km³ (area above ELA: 231 km², ~85.5% of the surveyed area) for the periods 2000–2013 and 2000–2016 and 0.273 km³ (area above ELA: 210 km², ~88.6%) for the period 2013–2016. The corresponding mean area weighted penetration depth uncertainty is equal to 0.065 m a⁻¹ for the period 2000–2016, which is comparable to the uncertainty of 0.081 m a⁻¹ reported by Malz and others (2018) for a similar observation period (15.86 years).

A strongly negative mass balance of −1030 ± 830 kg m⁻² a⁻¹ for the RGI6.0 region ‘Low Latitudes’ (Tropics) is reported by Zemp and others (2019) for the period 2006–2016. However, this estimate is based on a very sparse data coverage (<2%), leading to the large uncertainty. The authors expressed the need for more observations in the Tropical Andes, which will be consequently provided by this study. Braun and others (2019) observed a mass balance of −340 ± 80 kg m⁻² a¹ (ρ = 850 kg m⁻³) for the southern wet outer tropics in 2000–2013, which corresponds to ~22% lower mass loss rate than revealed in this study (−437 ± 122 kg m⁻² a⁻¹, first density scenario). As discussed in Section 5.1, the RGI6.0 glacier inventory, which was used by Braun and others (2019), contains many misclassified glacier areas that consequently led to lower specific mass loss. The analysis at the central Cordillera Real by Soruco and others (2009a) revealed a volume loss of 43% in the period 1963–2006 based on the elevation change measurements at 21 glaciers. The authors observed also large differences for the measured glaciers, ranging from −260 to −1380 kg m⁻² a⁻¹ (ρ = 900 kg m⁻³). Moreover, they found strong temporal variations ranging from an average mass balance of 700 kg m⁻² a⁻¹ for 1963–1975 to −500 kg m⁻² a⁻¹ for 1975–1983. For the period 1997–2006, an average mass balance of −450 kg m⁻² a⁻¹ was derived for the surveyed glaciers. This value is similar to our observations of −437 ± 122 and −499 ± 97 kg m⁻² a⁻¹ for the periods 2000–2013 and 2000–2016, respectively (first density scenario). However, for the period 2013–2016, we measured a higher mass loss of −630 ± 450 kg m⁻² a⁻¹, which correlates well with our observed higher glacier retreat rates for this time interval (Fig. 8). An average glaciological mass balances at Zongo Glacier of −279 and −273 kg m⁻² a⁻¹ is reported for the period 2000–2013 and 2000–2016, respectively (WGMS). These values are comparable to our findings of −275 ± 70 and −365 ± 60 kg m⁻² a⁻¹ in the respective intervals at Zongo Glacier. At Charquini Sur Glacier, we quantified a mass balance of −839 ± 160 kg m⁻² a⁻¹ in the period 2000–2016, which fits the average measured glaciological mass balance of −967 kg m⁻² a⁻¹ in the period 2002–2016 (WGMS). Chacaltaya Glacier disappeared within the observation period. Therefore, a comparison with glaciological mass-balance estimates is not meaningful.

Fig. 8. Temporal evolution of glacier area and specific mass balance (spMB) (a and b) as well as Oceanic Nino Index (ONI, c) for the studied period.

At Zongo Glacier, Soruco and others (2009b) found that 80% of the specific mass-balance contribution comes from the strongly ablating low-elevation area (about one-third of the glacier area). This fits well with our observed pattern on elevation change and altitude (Figs 6, S8 and S9) and also explains the high retreat rates of the glacier tongues. Figures 7, S11 and S12 indicate that small, low-lying and northward oriented glaciers showed predominantly more negative mass balances. Similar to Soruco and others (2009b), we correlated the glacier orientation and mass budgets. The applied sinusoidal fits show minimum mass balances for North/North-West orientations (Fig. S14), which is in accordance with observations by Soruco and others (2009b) who reported the highest mass losses on glaciers facing toward North/North-West since the annual mean incoming solar radiation is highest for north-facing slopes at these latitudes. In particular, some of the smallest glaciers show highest specific mass loss variability since they are strongly affected by local settings (Vuille and others, 2008). The air temperature at surrounding rock surfaces can rise above 20 °C and lead to advection of warm air over the glacier (Francou and others, 2003). This so-called ‘edge effect’ becomes important once the glacier shrunk below a critical size and can foster accelerated glacier recession. This pattern also fits with our observed higher relative area changes for small low-lying glaciers (Section 5.1).

Rabatel and others (2013) derived two times higher mass loss rates for glaciers with maximum elevations below 5400 m a.s.l. as compared to glaciers with maximum altitudes above this limit, based on mass-balance series of eight glaciers throughout the Tropical Andes (three glaciers in Bolivia). Such a distinct pattern is not so obvious in the geodetic measurements obtained in this study. However, we revealed an overall trend toward higher mass loss rates for glacier with a lower upper limit, and 12–24% stronger mass losses for glaciers with maximum altitudes below 5400 m a.s.l. (Fig. S13). On the one hand, our analysis covers only Cordillera Real and Tres Cruces, on the other hand, the representativeness of the eight glaciers, analyzed by Rabatel and others (2013), for all glaciers in the Tropical Andes is not guaranteed. Moreover, glaciological mass-balance measurements can be biased by the spatial data sampling (see above). Thus, we attribute the offset in the trends to differences in the spatial coverage and sampling as well as to the different mass-balance method used.

As mentioned in Section 4.2, data voids in InSAR-based DEMs can occur due to topographic limitations like SAR shadow and layover. For the study periods 2000–2013 and 2000–2016, the SRTM DEM is used as a reference. Only 63.3% of the glacierized area in 2000 is covered by SRTM elevation data. In combination with the TanDEM-X coverage in 2013, we retrieve Δh/Δt measurements on 59.2% of the glacier area. This corresponds to coverage of 93.5% of the SRTM data by TanDEM-X data. The largest data voids are present in the Δh/Δt fields for 2000–2016 (52.0%) since no complete TanDEM-X coverage of our study site is available in 2016 (~47 km² ≙ 20% of the glacier area in 2013 is not covered, Fig. 1). However, the amount of voids in the Δh/Δt data does not significantly influence the region-wide mass-balance estimates when applying ‘global hypsometric interpolation’. According to McNabb and others (2019), this method is suitable to compute regional averages from Δh/Δt data sets with up to 60% of voids. For the study period 2013–2016, the mosaicked TanDEM-X DEM from 2013 is used as a reference, leading to the highest coverage of 69.2% (62.3% after filtering) of all analyzed glacierized areas by Δh/Δt data. Moreover, coverage of 86.8% by Δh/Δt data on the glacier areas within the TanDEM-X footprints in 2013 and 2016 (regions where TanDEM-X data from both years is available, see also Fig. 1) is revealed for 2013–2016. Consequently, the achievable data coverage can be increased when comparing TanDEM-X to TanDEM-X DEMs, especially in alpine regions where SRTM data can have large voids. Future studies can benefit from using TanDEM-X DEMs instead of SRTM DEM as a reference since less extrapolation is needed, helping to increase the accuracy (Section 3.3). Moreover, the amount of individual glacier mass-balance estimates is limited by voids in the Δh/Δt fields (Section 3.2). Thus, a higher spatial coverage will also increase the amount and quality of glacier-wide mass budget information. Combining TanDEM-X data from ascending and descending orbits is highly recommended if available for the specific region and time in order to obtain optimum data coverage. Moreover, many published studies used SRTM DEMs to compute elevation change rates. However, it is often not clear how the void-filled areas (no timestamp is available) were handled. Therefore, we emphasize to clearly document for each study which SRTM DEM version is used and especially how the void-filled areas are handled when using SRTM DEMs to assess elevation change rates.

5.3 Climatic effects on glacier changes

Previous surveys at Zongo and Chacaltaya glaciers (e.g. Wagnon and others, 2001; Francou and others, 2003) demonstrated that annual surface mass balance is dominated by ablation and accumulation processes during the wet summer season (70% of the total variance at Zongo Glacier). For the observation period 2000–2013, the glacier mass balance is computed by comparing DEMs from the same seasons (around February, see Table 1) to avoid a potential bias due to intra-annual mass-balance fluctuations and to minimize the impact of differences in the SAR signal penetration due to different surface conditions. Unfortunately, in 2016, TanDEM-X data are only available for the period September–November (see Table 1), the end of the dry season and the beginning of the wet season. Minimal precipitation occurs typically during the dry season and the glacier mass balance is dominated by ablation (Kaser, 2001; Favier and others, 2004). Thus, we would like to point out that the presented mass balances for 2000–2016 and 2013–2016 could potentially be slightly more positive when accounting for the seasonal bias. However, it is difficult to correctly quantify such a bias. Therefore, no seasonal correction is applied.

The inter-annual variability of the surface mass balance is strongly influenced by ENSO events on regional scales (Francou and others, 2003) (see also Section 1). In order to evaluate the ENSO effects on the glacier changes, we used as an indicator the monthly Oceanic Niño Index (ONI) provided by NOAA Climate Prediction Center (http://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php). According to NOAA's definitions, El Niño is present for ONI values above +0.5 and La Niña conditions are indicated for ONI values below −0.5. The temporal evolution of ONI is plotted in Figure 8 for the period 2000–2016. A strong El Niño event is clearly visible for 2015 (ONI values of up to 2.6), which leads to a clear positive average ONI value of 0.42 in the period 2013–2016. In the period 2000–2013, the ONI fluctuations are more balanced with a mean ONI value of −0.17. More negative mass balances, caused by precipitation deficits and higher incoming radiation energy, are typical for El Niño events in the Bolivian Andes (Francou and others, 2003; Maussion and others, 2015). A two-thirds higher runoff at Zongo Glacier was observed by Wagnon and others (2001) during the strong El Niño event in 1997–1998. Additionally, Soruco and others (2009a) revealed that 60% of the mass loss in 1997–2006 occurred during this event. Schoolmeester and others (2018) reported an extreme drought that stressed the water availability in Bolivia in 2015/16. Moreover, we analyzed monthly mean temperature recordings at Zongo and Charquini Sur glaciers (Figs S15 and S16). Higher average temperatures are obtained for the period 2013–2016, as compared to the period 2000–2013, and a clear peak in the temperature trend is obvious during the El Niño event in 2015/16. Thus, it is likely that our observed increase glacier shrinkage in the period 2013–2016, which is also reported by Seehaus and others (2019) for Peruvian glaciers, can be attributed to the strong El Niño event in 2015/16. Since the response time of tropical glaciers to mass-balance variations is very short (e.g. Francou and others, 2004; Rabatel and others, 2013) and the fact that the glaciers in the low latitudes are in average the thinnest worldwide (Farinotti and others, 2019), the increased melt in the period 2013–2016 explains the enhanced glacier recession in this interval. In the period 1986–1992, a mean ONI value of +0.24 (maximum values of 1.71 in January 1992 and 1.70 in August 1987) was obtained. Consequently, the very high glacier retreat rates obtained by Cook and others (2016) of 14.5 km² a⁻¹ in this period, which are similar to ours in the El Niño period 2013–2016, can be most likely attributed to increased El Niño activities. This link is in line with the observations by Morizawa and others (2013) who derived higher recession rates during El Niño and even partially slight area gains during La Niña at Condoriri Glacier, Cordillera Real in the period 1988–2010. Moreover, Silverio and Jaquet (2017) revealed a similar pattern and presented a linear correlation (R ² = 0.8) between ONI and glacier area changes. We correlated glaciological measurements at Zongo and Charquini Sur Glacier (Figs S17 and S18, WGMS) with average ONI values of the respective analysis period (September–August). Tendencies toward smaller accumulation-area-ratios (AAR), higher ELA values and more negative specific mass balances are found (Figs S19 and S20). These revealed tendencies fit the observations by the previously mentioned studies and support our assumption that the strong El Niño event in 2015/16 lead to the increased glacier wastage.

Vuille and others (2018) estimated an average temperature increase by 1°–5 °C until 2100 relative to the 1961–1960 mean value for the Tropical Andes, based on CMIP5 scenarios (RCP 4.5 and RCP 8.5). The authors also evaluated the impact on the ELA and projected for Zongo Glacier a rise toward ~5400 m a.s.l. (RCP4.5) and ~5700 m a.s.l. (RCP8.5). Taking into account our observed median glacier elevations and the relation of area changes to median glacier elevation (Section 5.1), such warming will further accelerate the glacier recession and lead to the disappearance of many small and low-lying glaciers and lead to significant volume losses as projected by Réveillet and others (2015), for example, 40–89% at Zongo Glacier (depending on the model scenario). The impact of future glacier recession on the basin runoff was evaluated by Huss and Hock (2018) on global scales. Their estimates indicated that the Bolivian glaciers have already reached or are close to ‘peak water’ (the point in time with maximum annual glacier runoff) and that a reduction of future glacier runoff in the range of 20–40% until 2090 has to be expected for the major drainage basins of the Bolivian glaciers, the Amazon and Titicaca basins. Comparable reduction in glacier meltwater contribution to the runoff of Zongo watershed by 2100 was modeled by Frans and others (2015). Large-scale and spatially detailed information is needed on the future evolution of glaciers, since the metropolis La Paz/El Alto obtains water from a large catchment area and further extensions are planned (personal communication with local stake holders). Thus, the results obtained in this study provided a consistent database for further regional and detailed analyses of glacier changes in Bolivia.

6. Conclusions

The results presented in this study clearly indicate the rapid glacier shrinkage in Bolivia. We observed area changes of −81 ± 18 km² (−29%, −5.1 km a⁻¹), the disappearance of 34 glaciers (reduction to an area <0.01 km²) and an average mass-balance rate of −403 ± 78 kg m⁻² a⁻¹ (−1.8 ± 0.5 Gt) in the period 2000–2016. An increased area change rate of −15 ± 5 km² a⁻¹ is revealed for 2013–2016 and can be attributed to the strong El Niño event in this period. The more negative mass-balance rate of −488 ± 349 kg m⁻² a⁻¹ in this period (not significant, considering the uncertainties) and reported values for Peruvian glaciers support these findings. The revealed temporal variability of glacier changes highlights the impact of short-term climatic variations on glacier mass balances in the Tropical Andes. The analysis of the area and mass changes of individual glaciers revealed a trend of higher mass losses for northward oriented glaciers and increased glacier wastage for small glaciers at low altitudes confirming the predicted vanishing of small low-lying glaciers in this region. Our observations also confirm that glacier retreat rates are even higher in the Cordillera Tres Cruces than in the Cordillera Real, which we attributed to the lower average glacier elevations.

Additionally, we showed that gaps in the SRTM are the major source of voids in the glacier surface elevation change fields and that by using solely TanDEM-X data, a coverage of 78% could be achieved.

This study provides the first multi-temporal, region-wide and spatially detailed mass-balance estimation for Cordillera Real and Cordillera Tres Cruces combined with a temporally consistent analysis of glacier area changes. These findings form fundamental information for the modeling of future glacier evolution, for developing water resource management plans and for further monitoring of glaciers. The ongoing dramatic glacier wastage in Bolivia will likely cause serious socio-economic challenges in the next decades. On the one hand, it will lead to a reduction in glacier meltwater contribution to discharge and hence affect regional water availability for mining, hydropower production, irrigation or domestic water supply. On the other hand, it can foster the extension and formation of glacial lakes, leading to a potential increase in the risk of GLOFs endangering downstream communities. Thus, further detailed monitoring of the Bolivian glaciers is highly advisable from scientific but also due to socio-economic aspects.