Glacier change and glacial lake outburst flood risk in the Bolivian Andes

. Glaciers of the Bolivian Andes represent an important water resource for Andean cities and mountain communities, yet relatively little work has assessed changes in their extent over recent decades. In many mountain regions, glacier recession has been accompanied by the development of proglacial lakes, which can pose a glacial lake outburst ﬂood (GLOF) hazard. However, no studies have assessed the development of such lakes in Bolivia despite recent GLOF incidents here. Our mapping from satellite imagery reveals an overall areal shrinkage of 228 . 1 ± 22 . 8 km 2 (43.1 %) across the Bolivian Cordillera Oriental between 1986 and 2014. Shrinkage was greatest in the Tres Cruces region (47.3 %), followed by the Cordillera Apolobamba (43.1 %) and Cordillera Real (41.9 %). A growing number of proglacial lakes have developed as glaciers have receded, in accordance with trends in most other deglaciating mountain ranges, although the number of ice-contact lakes has decreased. The reasons for this are unclear, but the pattern of lake change has varied signiﬁcantly throughout the study period, suggesting that monitoring of future lake development is required as ice continues to recede. Ultimately, we use our 2014 database of proglacial lakes to assess GLOF risk across the Bolivian Andes. We identify 25 lakes that pose a potential GLOF threat to downstream communities and infrastructure. We suggest that further studies of potential GLOF impacts are urgently required.


Introduction
Tropical glaciers are sensitive indicators of climate change and provide important information on climatic trends in locations where meteorological observations are sparse (Kaser, 1999;Vuille et al., 2008;Soruco et al., 2009). In this study, we focus on glacier change in the Bolivian Andes, which contain ∼ 20 % of the world's tropical glaciers (Jordan, 1991(Jordan, , 1998Kaser, 1999). Here, glacial meltwater is an important water resource for major cities, such as La Paz and El Alto, as well as for mountain communities across the region (Vergara et al., 2007;Vuille et al., 2008;Oxfam, 2009;Rabatel et al., 2013;Rangecroft et al., 2013). La Paz and its neighbouring city of El Alto are home to ∼ 2.3 million people, and it is estimated that glaciers supply ∼ 15 % of the potable water supply to these areas . La Paz, like many Andean cities, derives some of its electricity from hydropower generation, which depends to some extent on glacial meltwater generation (Painter, 2007). Some researchers have expressed concern that power generation during the dry season will become unreliable due to low water flows (Painter, 2007;Chevallier et al., 2011;Kaenzig, 2015), although this requires quantitative study.
There is also some evidence that more isolated mountain communities in Bolivia are suffering increasingly from the adverse effects of glacier recession and changing meltwater supply in response to increasing atmospheric temperature (Oxfam, 2009), although it does not yet appear to be a direct driver of rural-to-urban migration (Kaenzig, 2015). Bolivia is the poorest country in South America and hence is very vulnerable to the impacts of climate change (Andersen Published by Copernicus Publications on behalf of the European Geosciences Union. and Verner, 2009;Winters, 2012). Indeed, it is estimated that only 56 % of Bolivia's rural population have access to safe water (Jeschke et al., 2012;Rangecroft et al., 2013), meaning that the sustainability of glaciers is a significant concern in the broader context of poverty and vulnerability to climate change.
Despite the regional importance of Bolivian glaciers, research to monitor their extent and response to climate change has been rather limited. Detailed mass balance and modelling studies have been performed for a few glaciers, such as Zongo Glacier (e.g. Sicart et al., 2011;Réveillet et al., 2015;Soruco et al., 2015). Other studies have documented the demise of Chacaltaya Glacier, which disappeared in 2009 (Ramirez et al., 2001;Soruco et al., 2015). At a broader scale, Jordan et al. (1980) and Jordan (1991Jordan ( , 1998 developed the first inventory of glaciers in Bolivia. Soruco et al. (2009) calculated a volumetric reduction of 43 % for 21 glaciers of the Cordillera Real over the period 1963-2006 (location shown in Fig. 1b), whilst Albert et al. (2014) demonstrated that glaciers of the small range of Tres Cruces had lost approximately half of their surface area between 1975 and 2009 (location shown in Fig. 1b). However, a comprehensive quantification of glacier change across the Bolivian Andes has hitherto not been undertaken, and little is known about broadscale glacier change in the last decade since the 1963-2006 study period of Soruco et al. (2009).
A crucial further issue for consideration is the development of potentially dangerous glacial lakes. Glaciers tend to erode subglacial basins and deposit eroded materials around their margins as lateral-frontal terminal moraines (Cook and Swift, 2012). Recession into these basins and behind impounding moraines causes meltwater to pond as proglacial and supraglacial lakes (Carrivick and Tweed, 2013;Cook and Quincey, 2015). In many mountain ranges, these lakes represent a glacial lake outburst flood (GLOF) hazard, as moraine dam integrity reduces over time, leading to dam failure, and as mass movements of ice, snow, and rock from surrounding valley slopes impact lakes, leading to wave overtopping of moraine and bedrock dams (Clague and Evans, 2000;Richardson and Reynolds, 2000;Hubbard et al., 2005;Carey et al., 2012;Westoby et al., 2014;Haeberli et al., 2016). GLOF hazards have received almost no research attention in the Bolivian Andes, yet Hoffmann and Wegenmann (2013) documented a recent GLOF at Keara in the Apolobamba region in 2009, indicating that the potential impacts of such hazards have been overlooked. In this case, an ice-dammed lake burst, killing several farm animals, destroying cultivated fields, and washing away a road that left the village population cut off for several months. Fortunately, no human casualties occurred. In other locations, several studies have shown that glacier recession has been accompanied by an increase in the number and size of proglacial lakes (e.g. Carrivick and Tweed, 2013;Carrivick and Quincey, 2014;Komori, 2008;Loriaux and Casassa, 2013;Hanshaw and Bookhagen, 2014;López-Moreno et al., 2014;S. Wang et al., 2015), raising concerns that glaciated mountain regions are becoming more hazardous with respect to GLOFs. As yet, however, no study has quantified proglacial lake development in the Bolivian Andes to assess whether such a trend is prevalent here. Furthermore, no studies have assessed the extent to which these lakes represent a GLOF threat to downstream communities and infrastructure within Bolivia.
This study has three primary objectives to address these shortcomings: (1) to quantify glacier change in recent decades (from 1986 to 2014) across the Bolivian Andes; (2) to evaluate the development of proglacial lakes through this period; and (3) to provide an initial assessment of whether any existing proglacial lakes represent a GLOF hazard that may require further monitoring. The time period is chosen based on the availability of satellite imagery, plus a desire to extend the period of observations from previous work on glacier change in the region -the observations of Soruco et al. (2009) extended to 2006 for the Cordillera Real, and those of Albert et al. (2014) extended to 2009 for the Tres Cruces region. The spatial extent of our study is chosen in order to provide a first integrated assessment of glacier change across the whole of the Bolivian Cordillera Oriental.
2 Study region and methodology 2.1 Study region Figure 1 shows a map of Bolivia with the footprint of the satellite imagery used to map glacier change from 1986 to 2014. This study focuses on the glaciated areas of the Cordillera Oriental of Bolivia, which itself can be divided into a series of ranges. The most significant ice cover is in the Cordillera Real, in the middle part of our study region, and the Cordillera Apolobamba in the north, which straddles the Bolivia-Peru border. There are a series of smaller glaciated ranges further to the south of the main glaciated area of the Cordillera Real (including Huayna Potosí, Mururata, Illimani), and we consider all of these together as the Cordillera Real (after Jordan, 1991Jordan, , 1998. The southernmost Tres Cruces region holds another small, glaciated range ( Fig. 1). We follow previous studies (Jordan, 1991(Jordan, , 1998Soruco et al., 2009;Albert et al., 2014) in mapping across the three regions covered by the three satellite footprints: Cordillera Apolobamba, Cordillera Real, and Tres Cruces.

Mapping glacier and lake change
Landsat satellite imagery, with a spatial resolution of 30 m, was obtained from the United States Geological Survey (USGS) using the Earth Explorer interface (http:// earthexplorer.usgs.gov/). Data were obtained for the three footprints outlined in Fig. 1 for the years 1986, 1992, 1999, 2010, and 2014 (Table 1). Wherever possible, images were selected with minimal cloud cover, and all images were selected from the dry season when snow cover would have been at a minimum (which could be confused for glacier ice). Ice cover was identified automatically in the first instance by us-ing the TM3 / TM5 band ratio (i.e. the ratio of red to shortwave infrared) (e.g. Bolch et al., 2010) on atmospherically corrected imagery (Burns and Nolin, 2014). Some manual editing of the TM3 / TM5 output was also required. Firstly, lakes that had been misidentified as ice cover were removed. Secondly, all glacier polygons smaller than 0.05 km 2 were removed as they probably represented snow patches rather than glacier ice (Bolch et al., 2010). Thirdly, the areas identified as ice cover were checked against imagery in Google Earth, which has historic imagery stretching back to 2002 for much of the area considered in this study. Any other misidentified features (e.g. snow patches, areas of moraine or debriscovered ice) were then edited manually. We estimated uncertainty in our mapping following the method outlined by Hanshaw and Bookhagen (2014). They www.the-cryosphere.net/10/2399/2016/ The Cryosphere, 10, 2399-2413, 2016 assume that mapping errors are Gaussian, i.e. that ∼ 68 % (1σ ) of pixels will be subject to error. Uncertainty is calculated as where P is the measured glacier perimeter and G is the grid cell size. Uncertainties calculated using Eq. (1) are ∼ 10 %. Paul et al. (2013) found that error for measuring clean-ice glaciers, such as in this study, are of the order of ∼ 5 %. Lake extents were digitised manually with reference to both Landsat false-colour composites and NDWI (normalised difference water index) rasters (Huggel et al., 2002;Bolch et al., 2010). All visible lakes were mapped where they occurred within approximately 2 km of the glacier margins. Ultimately, however, we focussed on those lakes located within 500 m of the glacier margins because these lakes are likely to be of most relevance in terms of GLOF risk (see below). As with glacier mapping, we estimate that our lake mapping uncertainty is ∼ 5-10 %.

Assessing glacial lake outburst flood risk
Several studies have proposed schemes or criteria by which to identify potentially dangerous glacial lakes or to assess the consequences of GLOF events (e.g. Allen et al., 2009;Huggel et al., 2004;McKillop and Clague, 2007;Bolch et al., 2008). In this study, we first identified lakes from the 2014 imagery that have the potential to burst by using some of the criteria outlined by Bolch et al. (2008); we then undertook an initial assessment of the potential severity of the resultant flood event (i.e. by estimating flood peak discharge) following procedures outlined by Huggel et al. (2004) and McKillop and Clague (2007). Some data identified in these studies (such as moraine height and glacier velocity) were not available to us and so could not be integrated into the risk assessment at this stage. Figure 2 summarises the process of GLOF risk assessment.
We considered potential GLOF sources to be (1) lakes in direct contact with glaciers and (2) lakes within 500 m of the glacier margin that were also within 500 m of a steep (45 • or steeper) slope. Ice-contact lakes could be directly affected by glacier calving events that generate waves with the potential to overtop impounding moraine and rock basin slopes. Previous studies have considered lakes within 500 m of a glacier to be potential GLOF sources (e.g. W. Wang et al., 2011;S. Wang et al., 2015). Both ice-contact lakes and lakes within 500 m of a glacier could be impacted by ice and snow avalanches, which could also generate overtopping waves. These immediate proglacial areas are also affected by slope debuttressing associated with recent deglaciation and hence are more likely to experience paraglacial slope failures that could impact the lakes (e.g. Hubbard et al., 2005;Cook et al., 2013;Haeberli et al., 2016). The potential for mass movements is likely to increase as permafrost thaws (e.g. Haeberli et al., 2016); recent predictions by Rangecroft  2016) suggest that permafrost in the Bolivian Andes will disappear almost entirely by the 2080s. Bolch et al. (2008) considered slopes steeper than 45 • to be particularly hazardous in this regard. Hence, we generated a slope map from a 30 m resolution Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model Version 2 (GDEM2) (downloaded from http://reverb.echo.nasa.gov/reverb/) to identify slopes of 45 • or steeper that could shed material into proglacial lakes. Runout distances of ice and snow avalanches, rockfalls, rock avalanches, and debris flows vary widely but appear to cluster of the order of 10 2 to 10 3 m (Alean, 1985;Rickenmann, 1999Rickenmann, , 2005Copons et al., 2009). In the absence of detailed modelling of mass movement runout distances, we considered any proglacial lake within 500 m of a 45 • or steeper slope to be potentially dangerous, although we emphasise that the selection of these values is somewhat subjective.
Next, we removed lakes that were unlikely to yield large flood events if they were to burst. This assessment was based principally on lake area, assuming that the smaller the area, the less significant the potential flood volume, although clearly other factors, such as the triggering mechanism and downstream conditions, play crucial roles in determining flood magnitude. Unfortunately, there are no size criteria to determine the smallest lake size to include in any inventory of GLOF risk. We know that the damaging flood at Keara (Hoffmann and Wegenmann, 2013) resulted from the almost complete drainage of an ice-dammed lake that had a surface area of ∼ 34 000 m 2 . Hence, we set a lower lake size threshold of 30 000 m 2 to capture lakes of similar size to the one at Keara, although smaller lakes could generate damaging floods. This should be sufficient to capture most potentially dangerous lakes because all of the mapped lakes either are moraine-dammed or sit within glacially overdeep- The Cryosphere, 10, 2399-2413, 2016 www.the-cryosphere.net/10/2399/2016/ ened rock basins and hence are less likely to drain completely as was the case for the ice-dammed lake at Keara. Having identified all of the lakes with a potential outburst risk, we then identified all such lakes within that population that posed a risk of damage to human interests (e.g. homes, roads or infrastructure, cultivated fields). For the most part, this was achieved using the GIS database of roads and settlements available freely from GeoBolivia (http://geo.gob.bo/). We also cross-referenced this with a visual assessment of human features in Google Earth and Bing Maps, plus our own observations of a few of the sites. As a first assessment, we judged lakes to be dangerous where they had a direct hydrological connection to downstream infrastructure and communities (e.g. where a road or village was in the direct path of a GLOF event), and in general we searched for such infrastructure within ∼ 20 km downstream of the proglacial lakes (the complete drainage of the ice-dammed lake at Keara affected the channel for ∼ 10 km downstream). However, floods could propagate further than 20 km. Mapping of hydrological connectivity was achieved using the hydrological (i.e. flow-routing) tools within ArcGIS 10.2.2, with the ASTER GDEM2 as input data.
We used Google Earth imagery to assess whether the lakes identified from previous steps were moraine-or rockdammed. Moraine-dammed lakes are considered more dangerous because an initial trigger event, such as an avalancheinduced wave, could lead to breach incision in the moraine and hence enhanced drainage of the lake (e.g. Westoby et al., 2015). Lakes confined within rock basins are less likely to experience breach incision, but there are notable examples of bedrock-dammed lake outbursts (e.g. Dortch et al., 2011;Vilímek et al., 2015), and these lakes remain vulnerable to wave overtopping (e.g. Haeberli et al., 2016). Thus, both moraine-and rock-dammed lakes represent potential sources of GLOFs.
To assess the severity of GLOF events from our inventory of potentially dangerous lakes, we first estimated lake volumes. Cook and Quincey (2015) reviewed the empirical approaches that have been adopted in previous studies to model lake volume. The most popular approaches (e.g. Huggel et al., 2002) estimate lake volume based on a measurement of lake area. Cook and Quincey (2015) highlighted some issues in this approach. Specifically, volume-area relationships suffer from autocorrelation between volume and area, which gives an unrealistic impression of the strength of correlation between the two. Lake volume is typically calculated by multiplication between measured lake area and averaged or integrated lake depth. Consequently, it is more appropriate to use depth and area data, which can be measured independently, to estimate lake volume. Cook and Quincey (2015) also suggested that geomorphological context be considered when applying empirical techniques to estimate volume because different lake types (i.e. ice-dammed, moraine-dammed, thermokarst) can have different morphological characteristics. All of the lakes identified in our in-ventory either are moraine-dammed or sit within rock basins. Given these suggestions, we used the data set from Cook and Quincey (2015) to derive an empirical relationship between lake depth and area that is specific to moraine-dammed lakes of a similar area range to those found in this study: where D is mean lake depth in metres (m) and A is lake surface area in square metres (m 2 ). Equation (2) can be used to predict mean lake depth, which can be multiplied by measured lake area to calculate lake volume. We are not aware of any empirical formula that can predict the volume of lakes situated within rock basins. In the absence of any such formula, we use Eq.
(2) to provide a first-order estimation of their volumes as well as moraine-dammed lakes. Lake volume can be used to estimate peak discharge (Q max ). Huggel et al. (2002) collated several empirical models for estimating GLOF peak discharge from lake volume but ultimately adapted the relationship of Haeberli (1983) to give where t, time, is equal to 1000 s. We used Eq.
(3) to estimate peak discharge for the lakes identified as being potentially dangerous. It should be noted that peak discharges calculated using Eq.
(3) represent worst-case estimates. Lakes are ranked from lowest to highest peak discharges in Table 2.

Glacier change 1986-2014
Our results reveal that total glacier areal cover across the Bolivian Cordillera Oriental in 1986 was 529.3 ± 52.9 km 2 and that by 2014 this area had reduced to 301.2 ± 30.1 km 2 (Fig. 3a). This represents a total areal reduction of 43.1 % over the 28-year study period. If the Peruvian Cordillera Apolobamba glaciers are included in the data set, then the total glacier cover in 1986 was 626.5 ± 62.7 km 2 and in 2014 it was 351.7 ± 35.2 km 2 , representing a 43.9 % reduction. Figure 3a illustrates the reduction in overall glacier cover across this period. Rates of ice loss appear to vary across the study period, with an initially rapid shrinkage between 1986 and 1992 (14.5 km 2 a −1 ), relatively modest losses between 1992 and 1999 (5.1 km 2 a −1 ), strong ice shrinkage between 1999 and 2010 (8.1 km 2 a −1 ), and modest losses between 2010 and 2014 (4.0 km 2 a −1 ) (except for the Tres Cruces region).
For consistency with earlier studies, we present results in Fig. 3 of glacier areal change for separate glaciated mountain ranges, and Figs. 4-6 illustrate glacier change as a series of maps for each region. All mountain ranges show decreases in overall glacier area across the study period with a total loss of 43.1 % glacier cover in the Bolivian Cordillera Apolobamba www.the-cryosphere.net/10/2399/2016/ The Cryosphere, 10, 2399-2413, 2016

Proglacial lake development 1986-2014
Figure 7 illustrates how the number and areal cover of icecontact and proglacial lakes has developed between 1986 and 2014 across the three regions of the Bolivian Andes. Data are presented for ice-contact lakes (Fig. 7a, b) and for lakes within 500 m of the 1986 ice margin (Fig. 7c, d), which illustrates the cumulative change in proglacial lake number and area as the ice receded to its 2014 position. Figure 8 shows the total change in lake number and area across the Bolivian Andes. Figure 7a indicates that that there is no clear pattern in the number of ice-contact lakes that developed from 1986 to 2014, although there has been an overall decline in the number of ice-contact lakes for all three regions. The Tres Cruces region shows the clearest trend of decline, although there were very few ice-contact lakes here throughout the study period (4 in 1986 and 0 in 2014). For the Cordillera Apolobamba, the number of lakes increased from 4 to 10 lakes between 1986 and 1999, before declining to 1 lake in 2014. The Cordillera Real region experienced the reverse situation to the Cordillera Apolobamba between 1986 and 2010, with a decline from 15 to 8 ice-contact lakes between 1986 and 1999, and an increase to 12 lakes in 2010, before falling to 6 lakes in 2014. Figure 8a illustrates the overall decline in the number of ice-contact lakes across all regions, from 23 in 1986 to 7 in 2014. Figure 7b shows the change in area of ice-contact lakes for all three regions. The trend in ice-contact lake area for the Cordillera Apolobamba and Tres Cruces region follows the trend in the number of ice-contact lakes shown in Fig. 7a. The Cordillera Real experienced an overall 22 % increase in ice-contact lake area, from 0.9 ± 0.09 km 2 in 1986 to 1.1±0.11 km 2 in 2014, but peaked at 1.2±0.12 km 2 in 2010. This increase in lake area, even though lake number has fallen across the same period, is driven by the growth of a few large ice-contact lakes (e.g. Laguna Glaciar at the northern tip of the Cordillera Real and the large lake, Laguna Arkhata, beneath the summit of Mururata - Fig. 5). Figure 8b shows the overall trend in ice-contact lake area across the study period, indicating a very slight (10 %) overall increase from 1.0 ± 0.1 km 2 in 1986 to 1.1 ± 0.11 km 2 in 2014 (represented mostly by lakes in the Cordillera Real), with a peak at 1.6 ± 0.16 km 2 in 2010. The number of ice-contact lakes stayed relatively stable from 1992 to 2010 (Fig. 8a), yet the total area covered by these lakes increased (Fig. 8b). This is explained by growth of ice-contact lakes until 2010, followed by detachment, leading to an overall decrease in both lake number and area by 2014. Figure 7c shows that the number of proglacial lakes has increased from 1986 to 2014, for all regions, as the ice has drawn back from its 1986 position. These lakes fall within 500 m of the 1986 margin and hence represent the cumulative total of lakes for the three regions as glaciers have receded. The greatest number of lakes exist in the Cordillera Real, which saw a 47 % increase from 92 lakes in 1986 to 135 lakes in 2014. The Tres Cruces region saw a 67 % increase from 24 lakes in 1986 to 40 lakes in 2014. The Cordillera Apolobamba has seen an overall increase of 72 % from 29 to 50 lakes across the study period, although there was a peak in 2010 at 53 lakes. Figure 8a reveals that there has been an overall total increase of 55 % in proglacial lakes from 145 to 225 lakes. Figure 7d shows that the area covered by lakes within 500 m of the 1986 ice margin has increased, and it broadly reflects the pattern of lake number change illustrated in Fig. 7c. In the Cordillera Real, there has been a 54 % increase in proglacial lake area from 2.7 ± 0.27 km 2 to 4.1 ± 0.41 km 2 . The Tres Cruces region has seen a rather more modest increase of 15 % from 2.6 ± 0.26 to 2.9 ± 0.29 km 2 , and the trend has levelled off since 2010. Proglacial lakes in the Apolobamba region have seen a 51 % increase in area from 1.1 ± 0.11 to 1.7 ± 0.17 km 2 , although the peak occurred in 2010 and lake area has since decreased. Figure 8b shows that total lake area has increased by 38 % from 6.33 ± 0.63 to 8.73 ± 0.87 km 2 .

Identification of potentially dangerous lakes
In 2014, there were 137 lakes within 500 m of the 2014 ice margin, with a total area of 5.7 ± 0.57 km 2 . From this database, we identified 25 lakes that have the potential to be the source of damaging GLOFs (Table 2). Table 2 shows that these lakes either are moraine-dammed or sit within rock basins and have surface areas that range across 3 orders of magnitude from 32 800 to 1 355 700 m 2 . Complete drainage of the smallest lake would yield a flood with an estimated peak discharge of ∼ 600 m 3 s −1 . Complete drainage of the largest lake would yield an estimated peak discharge of 126 710 m 3 s −1 , although the nine largest lakes are contained within rock basins and would be unlikely to drain completely. We emphasise that Eq. (3), which is used to calculate peak discharges, yields worst-case estimates and that such severe floods are probably unlikely. Many villages, farms, and roads could be impacted by floods from the identified lakes. The lakes are also shown in Figs. 4-6, although, because of their size relative to the scale of the map, these are best viewed in the Supplement.

Bolivian glacier change
We make some comparisons with the limited previous research on Bolivian glacier change, although this is complicated to some extent because of inconsistent methodologies, different study periods, and inclusion or exclusion of glaciated areas in different inventories. Glacier change in the Cordillera Apolobamba has not been investigated previously, which represents a significant gap in our understanding of Bolivian (and some Peruvian) glaciers. We have provided the first assessment of glacier change in this region. Jordan (1991Jordan ( , 1998 reported glacier areal coverage in this region for 1984 as 219.8 km 2 . Our results from 1986 indicate ice coverage of 172.3 ± 17.2 km 2 for the Bolivian Cordillera Apolobamba, but this figure rises to 269.5 ± 27.0 km 2 when Peruvian glaciers are included (Figs. 3b, 4). The discrepancy between our results and those of Jordan (1991Jordan ( , 1998 could be explained to some extent by further ice recession between 1984 and 1986. However, on closer inspection, it appears that Jordan (1991Jordan ( , 1998 included all ice across the Chaupi Orko range (see Fig. 4 for location), both on the Bolivian and Peruvian sides of the border. When we include the same areas, the total for the Apolobamba Range is 221.3 ± 22.1 km 2 , which is more consistent with Jordan's (1991Jordan's ( , 1998 work. The slightly higher value is explained by our inclusion of some relatively small glaciers separate from the main glaciated ranges of the Cordillera Apolobamba, which were not included in Jordan's (1991Jordan's ( , 1998 mapping. Since 1986, the trend of glacier loss has been sustained throughout the study period, similar to other glaciated mountain ranges in Bolivia (Figs. 3, 4, 5, 6), with an overall glacier ice shrinkage of 43.1 % for the Bolivian Apolobamba and 45.6 % for the combined Bolivian-Peruvian Cordillera Apolobamba.
Whilst individual glaciers of the Cordillera Real have been the subject of intensive study over many years (e.g. Ramirez et al., 2001;Sicart et al., 2011;Réveillet et al., 2015 Jordan (1991Jordan ( , 1998 gives glacier ice cover of 323.6 km 2 for the Cordillera Real, which is consistent with our 1986 value of 315.2 ± 32.4 km 2 (Figs. 3c, 5). Glacier change in the Tres Cruces region has been investigated by Albert et al. (2014Albert et al. ( ) from 1975Albert et al. ( to 2009. Their results indicated ∼ 55 % areal loss over this study period, with a marked reduction in ice cover between 1975 and 1986, before the start of our monitoring period. Between 1986 and 2009, their results showed that ice cover had shrunk from ∼ 36 to ∼ 25 km 2 , similar to our results (Fig. 3d). Since then, our results have shown a further reduction of 12.4 % to ∼ 22 km 2 for 2014.
Overall, the glacier retreat rate across the Cordillera Oriental is 1.54 % a −1 , excluding Peruvian glaciers of the Cordillera Apolobamba, or 1.57 % a −1 if Peruvian glaciers are included. Regionally, the Tres Cruces experienced the highest retreat rates (1.69 % a −1 ), followed by the Cordillera Apolobamba (1.54 % a −1 excluding Peruvian glaciers and 1.63 % a −1 including Peruvian glaciers), with the Cordillera Real experiencing the lowest retreat rates (1.50 % a −1 ). These values are comparable to retreat rates measured elsewhere in the Andes. For example, Rabatel et al. (2013) report retreat rates in Ecuador between 1962 and 1997 of 1.6 % a −1 , and 2 % a −1 for Columbian glaciers from the late 1970s to early 2000s; retreat rates of between 0.34 and 2.05 % a −1 are reported by Vaughan et al. (2013) for glaciers in Peru across similar time periods. Retreat rates further south in Patagonia are reported to be much lower at 0.14-0.66 % a −1 (Vaughan et al., 2013); in extra-tropical mountain ranges, such as the Alps, comparable rates of retreat of between 0.59 and 2.07 % a −1 are reported, although retreat is generally reported to be lower than in the present study (Vaughan et al., 2013).
The trend in glacier shrinkage across the Cordillera Oriental is of some concern in terms of water resources across the region, and particularly for the major cities of La Paz and El Alto. According to Soruco et al. (2015), the 50 % loss of glacier cover between 1963 and 2006 had not led to a decrease in runoff at La Paz, indicating that the reduction in glacier area had, so far, been offset by increased melt rates. However, complete glacier disappearance would reduce runoff by ∼ 12 % annually and by up to 24 % during the dry season . Although glacier shrinkage is not yet known to have driven rural-to-urban migration (Kaenzig, 2015), this could be a further pressure in the future, and it is likely that the effects of glacier shrinkage will be felt both in large cities, which rely to some extent on glacial meltwater (Vuille et al., 2008;Rangecroft et al., 2013;, and in rural communities (Andersen and Verner, 2009;Oxfam, 2009;Winters, 2012). Even in the Cordillera Apolobamba, which is rather sparsely populated, there are still ∼ 5500 people who live within 10 km of the glaciers (according to GeoBolivia GIS data, http://geo.gob.bo/). Approximately 13 700 people live within 10 km of the Tres Cruces glaciers, and ∼ 30 000 people for the Cordillera Real. Future changes to glacial water supply are likely to be felt keenly within these immediate rural areas, where communities may depend to some extent on meltwater during the dry season for drinking water, crop irrigation, and sustaining livestock. However, glacial meltwater also supplies populations in villages and cities beyond the immediate vicinity of the glaciated mountains. Another adverse impact would be toward the bofedales ecosystems (high-altitude peat bogs and wetlands), which are also fed by glaciers and which represent important water stores (Garcia et al., 2007;Squeo et al., 2006). Long-term glacier monitoring of the ice masses that supply water to La Paz and El Alto have been crucial in terms of understanding the sustainability of glacier meltwater in the region (e.g. Réveillet et al., 2015), yet similar studies have not been undertaken in other mountain ranges across Bolivia (e.g. in the Cordillera Apolobamba and Tres Cruces) where local populations could be very vulnerable to The Cryosphere, 10, [2399][2400][2401][2402][2403][2404][2405][2406][2407][2408][2409][2410][2411][2412][2413]2016 www.the-cryosphere.net/10/2399/2016/ future glacier shrinkage (Andersen and Verner, 2009;Oxfam, 2009;Winters, 2012).
To provide a first-order estimate of future glacier evolution across Bolivia, we used the data presented in Figure  3a to derive an exponential function (because it provided the best fit with the data compared to linear and other bestfit lines) that could be used to model future glacier decay. This method indicates that glaciers across Bolivia will have shrunk to around 10 % of their 1986 area by ∼ 2100. Extrapolation of glacier areal decline trends can only represent a first-order approximation and masks the complex array of factors that determine glacier mass balance and volume, but our estimate suggests that further work is urgently required to accurately model glacier change and to assess the consequences of that change on people and mountain ecosystems. There are few studies that model glacier demise in Bolivia, but Réveillet et al. (2015) modelled the future evolution of Zongo Glacier, forcing the model with temperature changes predicted by the Coupled Model Intercomparison Project phase 5 (CMIP5). Their results indicated that Zongo Glacier would lose 69 ± 7 % of its volume by 2100 with the intermediate CMIP5 scenario, and 40 ± 7 and 89 ± 4 % with the extreme scenarios. Although we are comparing an individual glacier with glacier demise across the whole Bolivian Andes, our results are consistent with the upper end of these predictions.

Proglacial lake development
We have made the first evaluation of proglacial lake development across Bolivia. In general, the number of proglacial lakes and their areas (i.e. those that formed within 500 m of the 1986 ice margin) increased as glaciers have receded from their 1986 positions (Figs. 7c,7d,8a,8b). This shows that ice recession has revealed further basins that have filled with meltwater. Several studies have described similar trends of increasing number and size of proglacial lakes from a range of locations, both from ice-sheet and valley glacier contexts (e.g. Carrivick and Tweed, 2013;Carrivick and Quincey, 2014;Komori, 2008;Loriaux and Casassa, 2013;Hanshaw and Bookhagen, 2014;López-Moreno et al., 2014;Schomacker, 2010;W. Wang et al., 2015).
Conversely, there has been an overall decrease in the number of ice-contact lakes across the study period for all three regions (Figs. 7a, 8a). However, this change has been highly variable both spatially (across the three regions) and temporally. The changes in lake area have followed the changes in lake number for the Cordillera Apolobamba and Tres Cruces regions, but the Cordillera Real has experienced overall areal growth despite a reduction in lake number (Fig. 7b). This is because there are a few large ice-contact lakes in the Cordillera Real (e.g. Laguna Glaciar, Laguna Arkhata - Fig. 5) that have been growing rapidly as ice has receded, in accordance with most findings in other deglaciating regions (e.g. Carrivick and Tweed, 2013;Komori, 2008;Lo-riaux and Casassa, 2013;Hanshaw and Bookhagen, 2014;López-Moreno et al., 2014;Schomacker, 2010;W. Wang et al., 2015). Reductions in ice-contact lake number are explained by lakes becoming disconnected from the glacier. It is nonetheless intriguing that lake number has generally decreased across all regions (Figs. 7a, 8a), and lake area has decreased in the Cordillera Apolobamba and Tres Cruces regions (Figs. 7b, 8b).
A few studies have also found that ice-contact lakes have reduced in number and/or size over time. For example, Gardelle et al. (2011) found that proglacial lakes in the Karakoram had shrunk because glaciers had surged or experienced reduced mass loss. This cannot explain the trends observed in Bolivia, where there are no surge-type glaciers and all glaciers are shrinking. Emmer et al. (2015) found that some lakes in western Austria had shrunk as a consequence of sedimentation and changes in water supply from the glacier. This explains the loss of some of the lakes within 500 m of the 1986 ice margin but is less important for the evolution of ice-contact lakes. Apart from the documented GLOF case for the Apolobamba (Hoffmann and Wegenmann, 2013), we did not observe any further evidence for significant lake drainage events. Perhaps the simplest explanation is that, as the ice has receded, there is now a lower contact area between the ice and ice-marginal zone, and hence much less space within which ice-contact lakes can exist (the ice-contact perimeter reduced by 32.6 % between 1986 and 2014). Another possibility is that glaciers have now receded far behind their Holocene erosional maxima, where basins were carved out under thicker, faster-flowing ice (cf. Cook and Swift, 2012). Hence, as glaciers continue to recede, there are fewer deep basins being revealed that could provide accommodation space for ice-contact lakes to develop. Likewise, unlike the debris-covered glaciers of the Himalaya, which develop large terminal moraine complexes that enclose lakes as ice recedes (e.g. Hambrey et al., 2009), Bolivian glaciers are mostly clean-ice glaciers that do not generally develop large terminal moraines. Indeed, most of the potentially dangerous lakes are found in rock basins (Table 2).
Despite the recent trend of reducing ice-contact lake development, it should be emphasised that the change in icecontact lake number and area has been highly variable throughout the study period, indicating that future lake development could be unpredictable without further efforts to investigate proglacial lake appearance and evolution. New and large lakes could develop in the future. One promising avenue would be to measure or model glacier-bed topography in an effort to identify future lake locations (e.g. Frey et al., 2010;Linsbauer et al., 2016), some of which could represent a GLOF risk.

Glacial lake outburst flood risk
An emerging issue in the Bolivian Andes is the threat of possible GLOF events (Hoffmann and Wegenmann, 2013). From our 2014 data set of proglacial lakes, we undertook a first-pass assessment of lakes that could represent the greatest outburst flood hazard (Table 2; Supplement) and hence should be the subject of future monitoring and GLOF modelling studies. A total of 25 lakes were identified that were large enough, and sufficiently close to potential sources of ice or rock avalanches, to be considered a potential GLOF risk to downstream communities or infrastructure. Nine of the potentially dangerous lakes are moraine-dammed. Our estimations of peak discharge indicate potentially very damaging floods from all of the lakes identified in Table 2, although these values represent worst-case scenarios. Indeed, even relatively small lakes can generate damaging floods. Specifically, the ∼ 34 000 m 2 ice-dammed lake that drained at Keara in 2009 (Hoffmann and Wegenmann, 2013) produced a peak discharge of approximately 400 m 3 s −1 . This discharge value is calculated using the supplementary icedammed lake data in Cook and Quincey (2015) to derive an equation similar to Eq. (2), and then using Eq.
(3) to estimate peak discharge. The geomorphological evidence for the Keara flood can be seen in Google Earth for ∼ 10 km downstream. Hence, we recommend that future studies focus on more detailed modelling of potential floods from these lakes; more detailed hazard analysis taking into account the potential runout and inundation of the floods, as well as the vulnerability of the affected communities (e.g. Carey et al., 2012); and the monitoring of lake evolution and the development of new lakes. Bathymetric studies of the lakes would also be welcome in order to improve volume estimations (Cook and Quincey, 2015).

Conclusions
Glaciers of the Bolivian Andes represent an important regional water source, and hence there is significant concern with respect to the sustainability of that water supply in a changing climate. We have performed the first integrated study of glacier change across the Bolivian Andes. Our mapping from 1986 to 2014 revealed that there has been a reduction in glacier area from 529.3 ± 52.9 to 301.2 ± 30.1 km 2 across the study period, equivalent to a 43.1 % shrinkage. Proportionally, ice loss was greatest in the southernmost part of our study area (Tres Cruces), where glaciers lost 47.3 % of their area between 1986 and 2014. The Cordillera Real (middle part of our study area) represents the largest area of glaciation in Bolivia and lost 41.9 % of its ice cover, while the Cordillera Apolobamba in the north lost 43.1 % of its ice cover (or 45.6 % if the glaciers of the Peruvian Cordillera Apolobamba are also included). The trend in glacier reces-sion has generally been rapid and continuous throughout the study period. We undertook the first assessment of proglacial lake development in Bolivia. There has been a general increase in the number and size of proglacial lakes across the study period, in accordance with several studies of proglacial lake development elsewhere. However, the change in ice-contact lake number and area has been more complex throughout the study period. All regions show a net decrease in the number of ice-contact lakes through the study period, although this change has been highly variable. Whilst ice-contact lake area has experienced a net increase in the Cordillera Real, consistent with previous studies performed in other deglaciating mountain ranges, the area coverage in the Cordillera Apolobamba and Tres Cruces regions has decreased overall. It is unclear why these patterns have emerged, but the variability in lake number and area indicates that ongoing monitoring of proglacial lake development is required, especially in light of the potential for these lakes to burst and initiate glacial lake outburst flood events.
GLOFs represent a potentially serious threat in Bolivia, with recent reports of GLOF events from the Apolobamba region (Hoffmann and Wegenmann, 2013). We used our 2014 (most recent) inventory of proglacial lakes to provide a first assessment of the potential for GLOFs in Bolivia. Overall, we identified 25 lakes that pose a potential GLOF risk to downstream communities or infrastructure. Estimated peak discharges from these lakes range from ∼ 600 to ∼ 127 000 m 3 s −1 , although the values derived from modelling of peak discharges represent worst-case scenarios. Nine of the lakes are moraine-dammed, and these could be susceptible to complete drainage. Sixteen lakes are rockdammed, including the nine largest potentially dangerous lakes identified in this study. We recommend further monitoring of potentially dangerous lakes and modelling of GLOF hazards across Bolivia.

Data availability
Glacier shapefiles will be available soon from the World Glacier Monitoring Service and Randolph Glacier Inventory. In the meantime, all data are available from the lead author, Simon Cook -contact details are given on the title page.
The Supplement related to this article is available online at doi:10.5194/tc-10-2399-2016-supplement.
thank Wilfried Haeberli and an anonymous reviewer for their constructive reviews of our manuscript, as well as the thoughtful interactive comment by Mauri Pelto and editorial comments from Etienne Berthier.
Edited by: E. Berthier Reviewed by: W. Haeberli and one anonymous referee