Modelling glacial lake outburst flood impacts in the Bolivian Andes

The Bolivian Andes have experienced sustained and widespread glacier mass loss in recent decades. Glacier recession has been accompanied by the development of proglacial lakes, which pose a glacial lake outburst flood (GLOF) risk to downstream communities and infrastructure. Previous research has identified three potentially dangerous glacial lakes in the Bolivian Andes, but no attempt has yet been made to model GLOF inundation downstream from these lakes. We generated 2-m resolution DEMs from stereo and tri-stereo SPOT 6/7 satellite images to drive a hydrodynamic model of GLOF flow (HEC-RAS 5.0.3). The model was tested against field observations of a 2009 GLOF from Keara, in the Cordillera Apolobamba, and was shown to reproduce realistic flood depths and inundation. The model was then used to model GLOFs from Pelechuco lake (Cordillera Apolobamba) and Laguna Arkhata and Laguna Glaciar (Cordillera Real). In total, six villages could be affected by GLOFs if all three lakes burst. For sensitivity analysis, we ran the model for three scenarios (pessimistic, intermediate, optimistic), which give a range of ~ 1100 to ~ 2200 people affected by flooding; between ~ 800 and ~ 2100 people could be exposed to floods with a flow depth ≥ 2 m, which could be life threatening and cause a significant damage to infrastructure. We suggest that Laguna Arkhata and Pelechuco lake represent the greatest risk due to the higher numbers of people who live in the potential flow paths, and hence, these two glacial lakes should be a priority for risk managers.


3 Introduction
Climate change is driving glacier mass loss in most parts of the world (Zemp et al. 2015). In many cases, glaciers recede into bedrock basins, or behind moraine dams, leading to meltwater ponding in proglacial lakes; meltwater can also pond supraglacially, especially on debris-covered glaciers (Cook and Swift 2012;Cook and Quincey 2015;Carrivick and Tweed 2013). Several studies have documented increases in the number and size of such lakes in recent decades (Komori 2008;Carrivick and Tweed 2013;Hanshaw and Bookhagen 2014;López-Moreno et al. 2014;Wang et al. 2015a). Glacial lakes can burst as a consequence of reduction in the integrity of the dam over time (e.g. from piping, mass movements on the dam flanks, or melting of ice within the dam), impact from ice, snow or rock avalanches entering the lake, glaciers calving into the lake, or seismic activity (Richardson and Reynolds 2000;Westoby et al. 2014b). Consequently, there has been growing concern about the impacts of glacial lake outburst floods (GLOFs) on downstream communities and infrastructure in many deglaciating regions worldwide (e.g. Vilímek et al. 2013;Carrivick and Tweed 2016;Emmer et al. 2016).
GLOF risk in Bolivia has received very little attention despite an event having occurred here in 2009, when an ice-dammed lake drained catastrophically, impacting the village of Keara in the Cordillera Apolobamba (Hoffmann and Weggenmann 2013). Recent work by Cook et al. (2016) found that glacier area in the Bolivian Andes had reduced by ~ 43% from 1986 to 2014, leading to an increase in the number and areal coverage of proglacial lakes. They undertook a rudimentary assessment of GLOF risk in the region, suggesting that 25 lakes could pose a risk to downstream communities. Subsequent work by Kougkoulos et al. (2018) developed a GLOF risk assessment technique based on the use of multicriteria decision analysis (MCDA), which was applied to the lake inventory of Cook et al. (2016). Kougkoulos et al. (2018) found that two lakes represented 'medium' risk and one lake represented 'high' risk. In this study, we aim to model potential GLOF inundation downstream from these three lakes and assess the possible impacts on downstream communities and infrastructure.

Geographical setting
This study focused on the glaciated ranges of the Cordillera Oriental, Bolivia; from north to south these are the Cordillera Apolobamba, the Cordillera Real (including Huayna Potosí, Mururata, Illimani), and the Cordillera Tres Cruces (Fig. 1). The three lakes identified by Kougkoulos et al. (2018) as having the highest risk include Laguna Arkhata, Pelechuco Lake, and Laguna Glaciar (Fig. 1). The only documented GLOF in Bolivia occurred at Keara (Fig. 1), which serves as a test case for our modelling approach. The Keara GLOF took place on 3 November 2009 at about 11:00 a.m. local time and involved the complete drainage of an ice-dammed lake (Hoffmann and Weggenmann 2013). This event flooded cultivated fields, destroyed several kilometres of a local dirt road, washed away pedestrian bridges, and killed a number of farm animals (Hoffmann and Weggenmann 2013). Fortunately, there were no human fatalities.

3
Information about the lakes discussed in this study (e.g. coordinates, measured areas, estimated volumes) is presented in Table 1. The lakes exist in the same region (Cordillera Oriental) and so are subject to similar environmental conditions (e.g. regional seismic activity, intense precipitation events, high-temperature events). However, these lakes were found by Kougkoulos et al. (2018) to be more susceptible than other glacial lakes in the Fig. 1 Location of glaciers and potentially dangerous glacial lakes in the Bolivian Andes, as well as the 2009 Keara GLOF event. The risk of each glacial lake identified by Cook et al. (2016) was assessed by Kougkoulos et al. (2018) using the MCDA methodology; three lakes were graded as 'medium' (orange) and 'high' risk (red) and form the focus of this study 1 3 region because of local conditions that control potential GLOF triggering mechanisms. We provide some information about these potential triggers in the subsequent subsections.

Pelechuco Lake
Pelechuco lake is located in the Cordillera Apolobamba ( Fig. 1) and, according to Kougkoulos et al. (2018), is considered a medium GLOF risk. This is due to the lake being in contact with steep (> 45°) surrounding slopes that could generate avalanches and/or rockfalls, which in turn could impact the lake and generate a displacement wave; this is the most common GLOF triggering mechanism in three different regions (Cordillera Blanca, North American Cordillera, Himalaya) (Emmer and Cochachin 2013). Further, the parent glacier is also in proximity (< 200 m) to the lake head, raising the possibility of ice calving into the lake, also generating a displacement wave.

Laguna Glaciar
Laguna Glaciar is located in the Cordillera Real ( Fig. 1). According to Kougkoulos et al. (2018), this is considered a medium-risk lake due to a visibly low dam freeboard (< 5 m), as well as the lake being in contact with the retreating parent glacier, which could calve into the lake. The surrounding steep slopes (30-45°) also raise the possibility of avalanches and/or rockfalls impacting the lake.

Laguna Arkhata
Laguna Arkhata is located in the Cordillera Tres Cruces (Fig. 1). According to Kougkoulos et al. (2018), it is considered the highest risk lake. This is mostly due to the steepest slope (> 45°) surrounding the lake, capable of shedding avalanches and/or rockfalls into the lake, as well as the lake being in contact with the retreating parent glacier, which could calve into the lake. Specifically there are numerous steep, hanging areas of glacier ice above the lake.

Topographic data acquisition
Topographic data are central to the modelling of flood inundation. For each lake, we generated DEMs from stereo and tri-stereo SPOT 6/7 satellite imagery (1.5-m horizontal resolution) using the DEM extraction module in ENVI (version 5.3) ( Table 2). For each of the four lakes and their GLOF runout zones, we sought to obtain at least two of the following types of image angles: one forward-looking (F), one nadir-looking (N), and one backward-looking (B) ( Table 2). We acquired images with across-track angles of between − 15 and + 15 degrees. For each of the three lakes and potential GLOF runout zones, we identified approximately thirty tie points over the pairs of images used for each DEM in order to calculate shift parameters between images at each tie point, obtaining an overall root-mean-square error (RMSE) of less than 1 pixel (< 1.5 m). In areas where there is a lack of field GPS data, some studies suggest the use of Google Earth for obtaining ground control points (GCPs) (Watson et al. 2015). Therefore, we collected 15 GCPs for each set of images. The combination of FNB images produced a 1.5-m resolution DEM for each site. The DEMs were resampled to 2-m spatial resolution in order to eliminate artificial roughness due to the stereo processing technique (Kropáček et al. 2015).

3
We also acquired an ASTER GDEM V2 with a 30-m horizontal resolution because, in Agua Blanca, a community situated downstream of Pelechuco lake, our SPOT-derived 2-m DEM contains a gap due to terrain shading. For this small area the ASTER GDEM was used instead. The ASTER GDEM has been used in several GLOF modelling studies (e.g. Gichamo et al. 2012;Wang et al. 2012;Watson et al. 2015) because it is free to download, and covers most of the Earth's landmass. In other cases, the DEM resolution has been resampled into a higher resolution in order to enable a higher level of precision in GLOF modelling, although the underlying DEM accuracy is not improved (e.g. Bajracharya et al. 2007;Allen et al. 2009;Rounce et al. 2016). In this study, we also applied a hydrological correction by filling sinks on both high (SPOT 2 m)-and low (ASTER GDEM 30 m)resolution DEMs.

Estimating lake volume
Lake volume defines the maximum amount of water that could be involved in a GLOF. Since bathymetric data were unavailable for the lakes examined in this study, lake volume had to be estimated using lake area measured from satellite imagery and a mean depth based on empirical equations. Specifically, Cook and Quincey (2015) found that the relationship between lake area, depth and volume varied depending on the geomorphological context of the lake. Using the empirical data set compiled by Cook and Quincey (2015), the mean depth (Dm) of ice-dammed lakes (i.e. Keara in this study) is predicted by the equation: where A is lake area. This relationship yields an R 2 value of 0.90 and a P value < 0.0001 (Appendix A1).
The mean depth of moraine-dammed lakes (i.e. Pelechuco and Laguna Glaciar in this study) is predicted by the equation: This relationship yields an R 2 value of 0.83 and a P value < 0.0001 (Appendix A1). Laguna Arkhata is bedrock-dammed, but as yet there are no depth-area relationships that exist for such lakes. Therefore, lake depth for Laguna Arkhata was calculated using Eq. 2.
We calculated the errors on the constants within each of the two linear regression models in Eqs. 1 and 2 (at a confidence level of 95%) and applied those errors to determine a range of lake depths for the three lakes under investigation (Appendix A1, Table 1). Lake volume (V) was estimated by multiplying the measured lake area (in Google Earth Pro) by the derived lake depth from Eq. 1 or 2, and the errors associated with lake volume were estimated by propagating the errors associated with lake depth.

Estimating peak discharge
Several empirical equations have been applied in the past to estimate peak discharge from natural-dam (e.g. moraine, bedrock) failures (see Table 3 in Westoby et al. 2014a). These relationships require inputs such as lake area, lake volume, water depth at dam, dam height (e.g. Williams 1978;Froehlich 1995;Pierce et al. 2010). Some of these parameters can be (2) Dm = 3 × 10 −5 A + 12.64 Table 3 Potential GLOF impact scenarios for all six villages under threat 1 For Agua Blanca, due to gaps on the 2-m DEM extraction we made use of ASTER 30-m GDEM to estimate potential impacts 2 2015 data set from GeoBolivia, which refers to the 2012 population census from the Bolivian National Statistical Institute 3 In parenthesis, > 2 m depth flood impacts on buildings and population Lake  (432) 979 (944) 1140 (843) 1260 (956) 2202 (2119) 1 3 assessed remotely (e.g. from satellite imagery or DEMs), whereas others require field investigation. For simplicity, we chose to use two commonly employed relationships to model peak discharge from ice-dammed and moraine-dammed lakes. The equation of Walder and Costa (1996) for non-tunnel floods was used to model peak discharge (Q max ) of the ice-dammed lake failure at Keara, which we used as a test case for our hydrodynamic modelling (see below).
The relationship of Evans (1986) was used for the potentially dangerous morainedammed lakes of Pelechuco and Laguna Glaciar, as well as the bedrock-dammed Laguna Arkhata: All pro-glacial lake size estimations (e.g. area, depth, volume) can be found in Table 1.

Dam-breach hydrograph and modelling parameters
Numerous parameters (e.g. peak discharge, sediment load, channel roughness) can modify the flood inundation and flow behaviour of a simulated flood, leading to uncertainty in flood modelling and risk assessment (Anacona et al. 2015;Wang et al. 2015b;Kropáček et al. 2015). Table 1 presents a summary of the modelling scenarios undertaken, together with the accompanying values for each model parameter (i.e. lake volume, percentage of lake drainage, outburst duration, Manning's value, peak discharge). We discuss each of these parameters in turn.
We followed Fujita et al. (2013) in calculating potential flood volume (PFV) as the product of lake area and either the mean depth (Dm) or the potential lowering height (Hp); Fujita et al. (2013) recommend using the lower of these two values to calculate PFV. We compared our values for Dm (Sect. 2.3; Table 1) against those calculated for Hp. Hp is the amount of lake lowering expected during a GLOF assuming that, after GLOF initiation, incision into the moraine dam will proceed until the angle between the lake and the downstream area is lowered to 10° (termed the 'depression angle'). We created depression angle maps (Appendix A2) and mapped the steep lakefront area (SLA) 1-km downstream of each lake, following Fujita et al. (2013), in order to calculate Hp (see Appendix A2). In all cases, the value of Dm is lower than that calculated for Hp, and so Dm is used in the calculation of lake volume.
The volume of water drained from the lake is highly uncertain principally because this will be determined by the triggering mechanism, which could vary in nature and severity, and the resilience and integrity of the impounding dam. Hence, three different drainage volumes were modelled, which represent optimistic, intermediate, and pessimistic scenarios. Different values for the extent of lake drainage were defined, ranging from 20% of the lake volume for the optimistic scenario, 50% for intermediate, and 100% for the pessimistic scenario. These scenarios are guided by previous work on GLOF drainage volumes. For example, Anacona et al. (2014) studied 14 Patagonian moraine-dammed lakes before and after GLOF events and found that two lakes emptied entirely (100% drainage), four emptied by more than 50%, two by around 50%, and six lakes generated smaller outbursts with less than 20% drainage. In another study, Kropáček et al. (2015) modelled 125% lake drainage for their worst-case (pessimistic) 1 3 scenario, with the additional 25% drainage accounting for possible future lake expansion in response to continued glacier recession. However, all lakes in our study are already at or very close to their maximum filling capacity, and so it is not necessary to model larger lakes. The lake drainage value for Keara was set to 100% since the observational evidence indicates that the lake drained completely (Hoffmann and Weggenmann 2013). For the other three lakes, we suggest that it is unlikely that they would drain completely, but we use the 100% drainage scenario as the worst possible case.
Previous studies have recommended a flood duration between 1000 and 2000s based on empirical data from the Swiss Alps (Haeberli 1983;Huggel et al. 2002). Therefore, we used 1000 s and 2000 s as outburst duration for the pessimistic and optimistic scenarios, respectively, and 1500 s for the intermediate scenario. As Worni et al. (2012) and Westoby et al. (2014a) suggest, if there are no data on lake outflow duration and lake discharge, the outflow hydrograph can only be validated indirectly (e.g. with the Keara event as a validation point in this case). Although simulations of drainage duration can be tuned to fit an observed hydrograph, hydrograph forecasting is difficult because nonlinear flood physics make drainage duration sensitive to the initial conditions (Ng and Björnsson 2003), which are usually uncertain. All three of our potentially dangerous glacial lakes are either bedrock-or moraine-dammed, so we made the assumption that flood discharge would increase linearly to a peak, after which it decreases linearly to 0 m 3 /s over a time span equal to that of the rising limb; in other words, hydrographs were assumed to be triangular in shape as has been applied in previous GLOF modelling studies (Anacona et al. 2015;Kropáček et al. 2015;Somos-Valenzuela et al. 2015;Wang et al. 2015b) (Fig. 2). Some previous studies have considered that the higher the peak discharge, the longer the flood duration (Anacona et al. 2015;Wang et al. 2015b), whereas others assume a shorter flood duration for a higher peak discharge (Somos-Valenzuela et al. 2015). We chose to use the latter option such that our worst-case pessimistic scenario has the highest peak discharge and shortest flood duration, and the optimistic scenario has the lowest peak discharge and longest flood duration. Ice-dammed lake failures usually generate flood hydrographs with a relatively slow, exponentially rising limb, and a rapidly falling limb (Kingslake 2013). Nevertheless, we chose to use a triangular-shaped hydrograph for the Keara ice-dam breach because no detailed data were available concerning the event, and we can compare more readily with results from the three lakes that are yet to generate GLOFs. In-channel and floodplain Manning's values for each scenario were set to represent mountain streams without vegetation (Table 1), although we could only verify the type of terrain by performing an inspection on Google Earth and from our own fieldwork observations at Pelechuco, in the Cordillera Apolobamba.

Hydrodynamic modelling of GLOFs
Modelling of GLOF hazard and risk has been undertaken in several regions around the world, using a variety of different approaches. Some studies have used simple geometric models such as the modified single flow (MSF) (Huggel et al. 2002;Allen et al. 2009;Prakash and Nagarajan 2017), the random walk process (Mergili and Schneider 2011), or the Monte Carlo Least Cost Path (MC-LCP) (Watson et al. 2015;Rounce et al. 2016Rounce et al. , 2017 to make a rapid assessment of potential flood inundation. Such models require very little input data and can be applied to many lakes, but serve as a first-order assessment that cannot go as far as to produce realistic flood maps. Other studies have focused on more detailed modelling approaches using HEC-RAS 1D (Bajracharya et al. 2007;Dortch et al. 2011;Klimeš et al. 2014;Watson et al. 2015), HEC-RAS 2D (Anacona et al. 2015;Wang et al. 2015a, b), FLO 2D (Petrakov et al. 2012;Somos-Valenzuela et al. 2015), and BASE-MENT (Worni et al. 2012(Worni et al. , 2013Somos-Valenzuela et al. 2016), all of which require further data such as channel bed roughness, and the volume of water that could drain from the lake. These techniques are capable of generating inundation maps (area affected, runout distance, and depth of flow) that can be used to inform risk management and mitigation strategies.
Here, GLOFs were modelled (inundation area, arrival time, depth, and velocity) using the 2D US Army Corp of Engineers model, HEC-RAS 5.0.3 (http://www.hec.usace .army.mil/). HEC-RAS was used because it is downloadable free of charge, and has been employed successfully to model GLOF inundation in a number of previous studies (e.g. Bajracharya et al. 2007;Dortch et al. 2011;Klimeš et al. 2014;Anacona et al. 2015;Wang et al. 2015a, b;Watson et al. 2015). The 2D version of this model can simulate multidirectional and multi-channel flows, which are characteristic of GLOFs (Westoby et al. 2014a, b;Wang et al. 2015b;Watson et al. 2015). The model set-up included the definition of upstream and downstream boundary conditions, the creation of a grid with elevation data, the selection of Manning's roughness values, slope parameters, and the model spatial domain. The unsteady flow simulation was performed for all four Bolivian GLOF case studies in order to observe: (1) peak flow propagation, (2) flood inundation extent, and (3) the flood water depth.
We used our simulated dam-breach hydrographs for the upstream boundary condition of each model (see Sect. 2.5 for more details). The normal depth option was used for the downstream boundary condition. This latter option uses Manning's equation to estimate a stage for each computed flow. To use this method, the user is required to enter a friction slope for the reach close to the boundary condition. If no detailed data exist, the slope of water surface can be used as a good estimate for the friction slope (Brunner 2010). This type of boundary condition should be placed far enough downstream of the study reach (i.e. potentially impacted communities) such that any errors it produces will not affect the results of the GLOF runout area (Brunner 2010;Watson et al. 2015). Hence, we placed the boundary 2-km downstream of the last community potentially affected by a GLOF from each lake. In addition, even though HEC-RAS has the ability to simulate debris flows, GLOFs were simulated as clear-water flows due to a lack of information about the nature of stream beds. We acknowledge that GLOFs are very likely to erode and entrain debris as they propagate downstream, and possibly evolve into debris flows, although many other studies also model GLOFs as clear-water flows due to data constraints (e.g. Anacona et al. 2015;Kropáček et al. 2015;Wang et al. 2015b;Watson et al. 2015).
Our overall approach was to test the model against field observations of the 2009 Keara event to ensure that realistic flood depths and inundation extent were achieved. The same modelling methodology was then applied to the three potentially dangerous lakes identified by Kougkoulos et al. (2018).

Population and infrastructure data
Population and infrastructure data are required to assess potential GLOF impacts. These data were acquired from the GeoBolivia portal (http://geo.gob.bo/porta l/), which offers open access to the 2012 population census and infrastructure data of the Bolivian National Statistical Institute. To quantify the downstream impacts, we manually counted the number of buildings affected by the flood. For each community, we also divided the total population of the community by the number of buildings in the community to estimate the number of people per building. This enabled us to estimate the number of people impacted by the flood. However, we acknowledge that there are likely to be spatial differences in population within each community, and, because of seasonal migration within Bolivia (Oxfam 2009), there will be temporal population variations (which we have not considered further in this study).
Field observations in the Cordillera Apolobamba region in July 2015 demonstrate that building structures are usually single-storey residential dwellings made of unreinforced brick walls (Fig. 3). Roofs are mostly constructed from corrugated steel sheets. According to Reese et al. (2007), this type of structure is vulnerable to flood depths of greater than 2 m, and from observations of other types of extreme flood events, such as tsunamis, lahars, and debris flows, only the concrete floor is likely to remain intact after the passage of a ≥ 2 m flood (Reese et al. 2007). Hence, for each one of the three scenarios simulated for each lake, we estimate the impact on the downstream communities taking into account two inundation depths (1) the extent of the flood from the > 0 m flood depth and (2) the extent of the flood from the > 2 m flood depth.

Modelling the 2009 Keara GLOF
The pre-GLOF lake area for Keara was 0.034 km 2 , measured from an August 2005 satellite image in Google Earth. When multiplied with the lake depth (Eq. 1), the resulting lake volume is between 0.14 × 10 6 m 3 and 0.39 × 10 6 m 3 , where the range in values accounts for potential errors in the relationship between depth, area, and volume (see Sect. 2.2) ( Table 1). Figure 4 illustrates a simulation of the Keara event with the intermediate scenario (Table 1), as well as images taken shortly after the GLOF. Comparison of the incision and inundation displayed in the field photographs with the flow depths and inundation extents of the numerical model 1 3 shows good agreement. In location A, the flood depth was 4-6 mat the same location as the accompanying photograph, which shows a ~ 5-m-deep incision through proglacial till (Fig. 4). In location B, the model replicates flood inundation around farm buildings and walls with favourable comparison with the accompanying field photograph. In zone C (Figs. 4 and 5), we observe that infrastructure is only affected by flows that do not exceed 2 m, showing favourable comparison with the real event since no dwellings were destroyed during the GLOF.

Pelechuco lake
Two communities (Agua Blanca and Pelechuco) are situated downstream of Pelechuco Lake (Fig. 6). Agua Blanca is ~ 8.5 km from the lake and has a total population of 379  (Table 3). Agua Blanca is the only threatened community for which the ASTER GDEM is used in the hydraulic modelling because the SPOT-derived 2-m DEM contains a gap at this location due to terrain shading. Agua Blanca is also the smallest community (~ 30,600 m 2 ; Table 3), meaning that the 30-m ASTER GDEM offers a rather crude estimation of flood impact. Model results indicate that around 111 to 372 people could be affected by the flood, and between 69 and 335 of those people could be affected by potentially damaging flood waters of ≥ 2 m depth (Table 3). Damaging floods are considered here to be life-threatening events or events causing a significant damage to infrastructure. Figure 6 provides a visualisation of the inundated areas downstream of Pelechuco lake.
The village of Pelechuco is about three times larger than Agua Blanca in terms of population number (981 individuals) and sits ~ 12.5 km from the lake. The three scenarios show that between 410 and 959 individuals could be affected by a GLOF, and that between 259 and 934 people could be affected by damaging (≥ 2 m depth) floods (Table 3).

Laguna Glaciar
Laguna Glaciar is situated above the community of Sorata, which is home to ~ 2788 individuals and is also a popular hiking and climbing destination for tourists. The largest part of the community is situated between ~ 100 and 250 m above the meltwater channel and floodplain and so is largely safe from any immediate GLOF impacts. However, a small part of the town sits close to the channel. Our modelling reveals that between 67 and 90 people could  (Table 3; Fig. 7).

Laguna Arkhata
Of all the glacial lakes in Bolivia, Kougkoulos et al. (2018) found that Laguna Arkhata represented the greatest risk to downstream communities. Totoral Pampa is the closest community to the lake (~ 6-km downstream) and is situated mostly within the floodplain of the meltwater stream that drains from the lake. Our modelling revealed that a GLOF could affect between 256 and 272 people (Table 3; Fig. 8). Therefore, even the optimistic scenario could affect most of the community.
Tres Rios, situated ~ 2-km downstream from Totoral Pampa, is located entirely on the floodplain. Thus, our modelling revealed that most of the community could be affected by a flood and 250 to 300 individuals could be exposed to potentially damaging floods.
Khanuma, which is situated another ~ 2-km downstream of Tres Rios, could be impacted to a lesser extent than the other two villages close to Laguna Arkhata. Modelling revealed that between 33 and 206 individuals could experience some flooding, and 25 and 201 of those people could be affected by damaging floods (Table 3; Fig. 8).

Modelling approach
For the first time, we have undertaken a hydrodynamic modelling study of GLOFs in the Bolivian Andes. For the most part, this assessment was achieved using a 2-m resolution DEM derived from SPOT 6/7 imagery as data input and HEC-RAS 5.0.3 to model clearwater flood flows. Our approach of producing and using a 2-m resolution DEM in favour of using coarser resolution, but freely available DEMs, such as ASTER GDEM or SRTM, offered us the best spatial resolution data set available to drive our model without having to resort to field-based topographic data acquisition through differential-GPS (d-GPS), LiDAR, or Structure-from-Motion (SfM) surveys. Coarser resolution DEMs have been used routinely in GLOF modelling (e.g. Anacona et al. 2015;Wang et al. 2015a), but our approach offers arguably more robust results for relatively narrow streams and valleys that may, in some parts of the reach, be on a scale at or below the spatial resolution of coarser DEMs. Hence, our 2-m DEMs offer a good balance between spatial resolution and convenience of data access, while maximising the veracity of numerical modelling output and facilitating a more robust risk assessment of downstream GLOF impacts.
We employed HEC-RAS 2D to model GLOFs because it is available free of charge, works directly in ArcMap, which is convenient, and has been used successfully in previous studies to model GLOF inundation (e.g. Bajracharya et al. 2007;Dortch et al. 2011;Klimeš et al. 2014;Anacona et al. 2015;Wang et al. 2015a, b;Watson et al. 2015). The 2D version of the model has a number of advantages over 1D models. Firstly, it offers the ability to simulate multi-directional and multi-channel flows (in contrast, 1D models are only capable of routing flow in one direction, i.e. downstream). Second, it models superelevation of flow around channel bends, hydraulic jumps (i.e. in-channel transitions between supercritical and subcritical flow regimes), and turbulent eddying (Westoby et al. 2014a). These conditions are characteristic of GLOFs, so it is advantageous to be able to model them effectively. Furthermore, with a 2D model, assessing people and properties at risk becomes a simpler task as there is no need to interpolate 1D results across the floodplain.
An issue for consideration was that we used HEC-RAS to model clear-water flows when it is likely that sediment will be entrained into the flow, changing its rheology, and thereby affecting its travel distance and downstream impacts. Westoby et al. (2014a) indicate that debris flows (with sediment comprising at least ~ 20% of flow volume) are more common than clear-water flows in GLOFs from moraine-dammed lakes. Furthermore, debris-laden flows will travel shorter distances than clear-water flows (Huggel et al. 2004). Therefore, in some cases they may not reach as far as our modelled clear-water flows, but will induce more damage to population and infrastructure situated in close proximity to the pro-glacial lakes (Reese et al. 2007;Jakob et al. 2012). Nonetheless, our study is not unusual in modelling GLOFs as clear-water flows (e.g. Anacona et al. 2015;Kropáček et al. 2015;Wang et al. 2015b;Watson et al. 2015). HEC-RAS 2D offers the possibility to include sediment entrainment in the unsteady flow simulations, which may improve the representativeness of the modelled floods. This would require field data on the sedimentary nature of potential flood channels and dams, which would be worthwhile obtaining in future studies. Given these uncertainties, we tested the robustness of our modelling approach against field data following the ice-dammed lake outburst at Keara, in the Apolobamba region (Figs. 1, 4). Comparison of our modelled flood inundation and flow depth results with documentary evidence of GLOF impacts illustrates that the model produced realistic results (Figs. 4,5).
A limitation of the HEC-RAS model is that, although it can model sediment entrainment, it cannot model debris flows that may result from a GLOF (Wang et al. 2015a, b). Unfortunately, cost constraints and a lack of field data on channel and moraines sediments did not permit the use of models such as RAMMS, which have been used in previous studies to model GLOFs and debris flows originating from such events (e.g. Frey et al. 2016). Other options that have been used to model GLOFs include FLO-2D (e.g. Somos-Valenzuela et al. 2016), the PRO version of which comes at a cost, but which can be used to model debris flows, and BASE-MENT (e.g. Worni et al. 2012Worni et al. , 2013Lala et al. 2017), which is downloadable free of charge, but cannot model debris flows. A direct comparison of the outputs of each of these modelling tools is beyond the scope of this study, but it would be worthwhile undertaking such work in the future for the three lakes investigated here. Further, recent GLOF studies have attempted to model complex process chains in a multi-hazard proglacial lake context (e.g. Lala et al. 2017;Mergili et al. 2018), and it is suggested here that investigation of process chains for GLOFs from proglacial lakes in Bolivia would be a valuable undertaking in the future.
Some uncertainties in our modelling include the volume of the lakes, the appropriate Manning's value to use for the channel, and the population living within each community. Hence, a field study of lake bathymetries, the nature of stream sediments, and a survey of the number of community occupants, which may change seasonally, would be welcome. Kougkoulos et al. (2018) identified three glacial lakes across Bolivia that represent medium and high GLOF risk. Table 3 provides a summary of the potential impacts of GLOFs from those three lakes, and Figs. 6, 7, and 8 give a visualisation of flood inundation for each affected community. In total, there are six communities at risk, which are home to ~ 4975 people and 2209 buildings. Our modelling results indicate that between 1140 and 2202 people could be affected by GLOF events if all three lakes burst and between 843 and 2119 of those people could be impacted by potentially damaging floods where the flow depth is modelled to be ≥ 2 m. Figure 9 provides a summary and comparison of GLOF impacts across the study sites.

Pelechuco Lake
Our results indicate that Pelechuco lake should be a priority for risk mitigation efforts. Kougkoulos et al. (2018) ranked this lake as 'medium' risk, and the modelling results presented in the current study show that part of both downstream communities (Agua Blanca and Pelechuco) could be severely damaged by potential flooding. Pelechuco in particular stands out on Fig. 9 as having the greatest number of potentially impacted people across all three lakes investigated. This is despite the fact that Pelechuco lake is the smallest of the three lakes examined here, and hence generates the smallest floods (Table 3). The high number of potentially impacted people results from the high proportion of the village area that is situated along the riverbanks (Fig. 6). The impact of a flood from Pelechuco lake on the community of Agua Blanca may require further detailed field investigation because only the 30-m resolution ASTER GDEM was available to drive the HEC-RAS model at this section of the channel reach (Fig. 6).

Laguna Glaciar
Of the three lakes, Laguna Glaciar arguably represents the lowest priority for risk managers in Bolivia. Kougkoulos et al. (2018) ranked the lake as medium risk, and our results indicate that, because a large part of the community of Sorata is situated far from the floodplain, the impacts of a GLOF would be relatively low (Figs. 7,9). Cost-benefit analysis between relocation of the part of the community potentially hit by a flood, and a combination of sustaining GLOF awareness programmes and early warning systems could be worthwhile undertaking for local risk managers and stakeholders. Ongoing monitoring of the lake and its surroundings as the glacier continues to recede would certainly be valuable.

Laguna Arkhata
Laguna Arkhata is a bedrock-dammed lake, and hence the probability of the dam failing and the lake draining completely is low; the most likely scenario would be an ice, snow, or rock avalanche impacting the lake and generating a wave that overtops the bedrock dam. However, irrespective of the drainage scenario (optimistic to pessimistic), the numbers of people affected by GLOFs vary relatively little for Totoral Pampa and Tres Rios; the differences are more marked for Khanuma. All three communities are affected by floods because they are located on floodplains. The communities of Totoral Pampa and Tres Rios look particularly vulnerable to GLOFs Totoral Pampa would be the first to be hit by a GLOF, and no matter what scenario is used in the modelling (from pessimistic to optimistic), Tres Rios is almost completely inundated from all flood scenarios because of its position on a floodplain (Figs. 8,9). This provokes the question of whether any structural or non-structural measures could be applied in these communities (e.g. GLOF awareness programmes, flood defence construction, community relocation). Future studies could focus on analysing the benefit of a managed retreat, by basically relocating most of the community outside the risky areas which is an increasingly discussed technique in the context of climate change for numerous regions facing various types of natural hazards around the world (Hino et al. 2017).

Conclusion
For the first time, we have modelled potential glacial lake outburst floods (GLOFs) in the Bolivian Andes. This was achieved using high-resolution (2 m per pixel) DEMs derived from (tri-) stereo SPOT satellite imagery to drive a 2D hydrodynamic model (HEC-RAS 5.0.3). The modelling approach was first tested against field evidence from a documented GLOF event at Keara in the Apolobamba region in 2009 and was shown to reproduce realistic flood depths and inundations. The model was then applied to three lakes that have been identified previously as representing a GLOF risk to downstream communities (Kougkoulos et al. 2018). These are Pelechuco lake, in the Apolobamba region, and Laguna Glaciar and Laguna Arkhata, in the Cordillera Real. GLOFs from all three lakes were shown to represent potentially life-threatening and/or damaging events. In total, GLOFs from the three lakes could impact six downstream communities. As a sensitivity analysis, we ran our model for three scenarios, representing pessimistic, intermediate, and optimistic scenarios, by varying model parameterisation. These scenarios give a range in the number of impacted people. In total, between 1140 and 2202 people could be affected by flooding if all lakes were to burst, and between 843 and 2119 of those people could be exposed to damaging floods (flow depth ≥ 2 m). Laguna Arkhata and Pelechuco lake represent the greatest GLOF risk due to the large numbers of people who live in the potential flow paths, and hence should be a priority for risk managers.