2.1 Quilcay catchment and Cojup Valley

Lake Palcacocha is part of a proglacial system in the headwaters of the Cojup Valley in the Cordillera Blanca, Peru (Fig. 1). This system was – and is still – shaped by the glaciers originating from the southwestern slopes of Nevado Palcaraju (6264 m a.s.l.) and Nevado Pucaranra (6156 m a.s.l.). A prominent horseshoe-shaped ridge of lateral and terminal moraines marks the extent of the glacier during the first peak of the Little Ice Age, dated using lichenometry to the 17th century (Emmer, 2017). With glacier retreat, the depression behind the moraine ridge was filled with a lake, named Lake Palcacocha. A photograph taken by Hans Kinzl in 1939 (Kinzl and Schneider, 1950) indicates a lake level of 4610 m a.s.l., allowing surficial outflow (Fig. 2a). Using this photograph, Vilímek et al. (2005) estimated a lake volume between 9 and 11 million m³ at that time, whereas an unpublished estimate of the Autoridad Nacional del Agua (ANA) arrived at approx. 13.1 million m³. It is assumed that the situation was essentially the same at the time of the 1941 GLOF (Sect. 2.2).

Figure 1Location and main geographic features of the Quilcay catchment with Lake Palcacocha and the former Lake Jircacocha.

Figure 2The Quilcay catchment from Lake Palcacocha down to Huaráz. (a) Lake Palcacocha in 1939, two years prior to the 1941 event. (b) The site of former Lake Jircacocha with the breached landslide dam and the former lake level. (c) Breached moraine dam and 1941 GLOF deposits, seen from downstream. (d) Left lateral moraine of Lake Palcacocha with landslide area of 2003. (e) Panoramic view of Lake Palcacocha, with the breach in the moraine dam and the modern lake impounded by a smaller terminal moraine and two artificial dams. (f) Panoramic view of Huaráz with city centre and approximate impact area of the 1941 event. Note that a small part of the lowermost portion of the impact area is hidden behind a hillslope. Photos: (a) Hans Kinzl, 1939 (Kinzl and Schneider, 1950); (b) Martin Mergili, July 2017; (c) Gisela Eberhard, July 2018; (d) and (f) Martin Mergili, July 2017.

The Cojup Valley is part of the Quilcay catchment, draining towards southwest to the city of Huaráz, capital of the department of Ancash located at 3090 m a.s.l. at the outlet to the Río Santa Valley (Callejon de Huaylas). Huaráz had a population of <20 000 at that time (Frey et al., 2018). The distance between Lake Palcacocha and Huaráz is approx. 23 km, whereas the vertical drop is approx. 1500 m. The Cojup Valley forms a glacially shaped high-mountain valley in its upper part whilst cutting through the promontory of the Cordillera Blanca in its lower part. Downstream from Lake Palcacocha by about 8 km (15 km upstream of Huaráz), the landslide-dammed Lake Jircacocha (4.8 million m³; Vilímek et al., 2005) existed until 1941 (Andres et al., 2018). The remnants of this lake are still clearly visible in the landscape in 2017, mainly through the change in vegetation and the presence of fine lake sediments (Fig. 2b). Table 1 summarizes the major characteristics of Lake Palcacocha and Lake Jircacocha before the 1941 GLOF.

Table 1Characteristics of Lake Palcacocha (1941 and 2016) and Lake Jircacocha (1941), and changes due to the 1941 GLOF. Topographic reconstruction according to field observations and historic photographs. Additional information from Vilímek et al. (2005) and ANA (2016).

¹ Reference values differ among sources: according to Vilímek et al. (2005), the volume of Lake Palcacocha in 1941 was 9–11 million m³, whereas a reconstruction of ANA resulted in 13.1 million m³. In contrast, Vilímek et al. (2005) estimate a pre-failure volume of 4.8 million m³ for Lake Jircacocha, whereas, according to ANA, the volume was only 3.0 million m³. ² Computed from the difference between the pre-1941 lake level and the modern lake level (before mitigation) of 4563 m. A reconstruction of ANA in 1948 resulted in a residual lake volume of approx. 100 000 m³ and a residual depth of 17 m, both much smaller than values derived through the reconstruction in the present work. One of the reasons for this discrepancy might be the change of the glacier in the period 1941–1948. ³ This value is highly uncertain and might represent an overestimation: the maximum depth of the lake strongly depends on the exact position of the glacier terminus, which was most likely located in an area of increasing lake depth in 1941.

Download Print Version | Download XLSX

Figure 3Situation in 1948, seven years after the 1941 event. (a) Residual Lake Palcacocha, and traces of the 1941 event. (b) Huaráz with the impact area of the 1941 event. Imagery source: Servicio Aerofotogramétrico Nacional, Peru.

2.2 1941 multi-lake outburst flood from Lake Palcacocha

On 13 December 1941 part of the city of Huaráz was destroyed by a catastrophic GLOF-induced debris and mud flow, with thousands of fatalities. Portocarrero (1984) gives an estimate of 4000 deaths, Wegner (2014) an estimate of 1800; but this type of information has to be interpreted with care (Evans et al., 2009). The disaster was the result of a multi-lake outburst flood in the upper part of the Cojup Valley. A sudden breach of the dam and the drainage of Lake Palcacocha (Fig. 2c and e) led to a mass flow proceeding down the valley. Part of the eroded dam material, mostly coarse material, blocks, and boulders, was deposited directly downstream from the moraine dam, forming an outwash fan typical for moraine dam failures (Fig. 2c), whereas additional solid material forming the catastrophic mass flow was most likely eroded further along the flow path (both lateral and basal erosion were observed; Wegner, 2014). The impact of the flow on Lake Jircacocha led to overtopping and erosion of the landslide dam down to its base, leading to the complete and permanent disappearance of this lake. The associated uptake of the additional water and debris increased the energy of the flow, and massive erosion occurred in the steeper downstream part of the valley, near the city of Huaráz. Reports by the local communities indicate that the valley was deepened substantially, so that the traffic between villages was interrupted. According to Somos-Valenzuela et al. (2016), the valley bottom was lowered by as much as 50 m in some parts.

The impact area of the 1941 multi-GLOF and the condition of Lake Palcacocha after the event are well documented through aerial imagery acquired in 1948 (Fig. 3). The image of Hans Kinzl acquired in 1939 (Fig. 2a) is the only record of the status before the event. Additional information is available through eyewitness reports (Wegner, 2014). However, as Lake Palcacocha is located in a remote, uninhabited area, no direct estimates of travel times or associated flow velocities are available. Also the trigger of the sudden drainage of Lake Palcacocha remains unclear. Two mechanisms appear most likely: (i) retrogressive erosion, possibly triggered by an impact wave related to calving or an ice avalanche, resulting in overtopping of the dam (however, Vilímek et al., 2005 state that there are no indicators for such an impact), or (ii) internal erosion of the dam through piping, leading to the failure.

2.3 Lake evolution since 1941

As shown on the aerial images from 1948, Lake Palcacocha was drastically reduced to a small remnant proglacial pond, impounded by a basal moraine ridge within the former lake area, at a water level of 4563 m a.s.l., 47 m lower than before the 1941 event (Fig. 3a). However, glacial retreat during the following decades led to an increase in the lake area and volume (Vilímek et al., 2005). After reinforcement of the dam and the construction of artificial drainage in the early 1970s, a lake volume of 514 800 m³ was derived from bathymetric measurements (Ojeda, 1974). In 1974, two artificial dams and a permanent drainage channel were installed, stabilizing the lake level with a freeboard of 7 m to the dam crest (Portocarrero, 2014). By 2003, the volume had increased to 3.69 million m³ (Zapata et al., 2003). In the same year, a landslide from the left lateral moraine caused a minor flood wave in the Cojup Valley (Fig. 2d). In 2016, the lake volume had increased to 17.40 million m³ due to continued deglaciation (ANA, 2016). The potential of further growth is limited since, as of 2017, Lake Palcacocha is only connected to a small regenerating glacier. Further, the lake level is lowered artificially, using a set of siphons (it decreased by 3 m between December 2016 and July 2017). Table 1 summarizes the major characteristics of Lake Palcacocha in 2016. The overall situation in July 2017 is illustrated in Fig. 2c.

2.4 Previous simulations of possible future GLOF process chains

Due to its history, recent growth, and catchment characteristics, Lake Palcacocha is considered hazardous for the downstream communities, including the city of Huaráz (Fig. 2e). Whilst Vilímek et al. (2005) point out that the lake volume would not allow an event comparable to 1941, by 2016 the lake volume had become much larger than the volume before 1941 (ANA, 2016). Even though the lower potential of dam erosion (Somos-Valenzuela et al., 2016) and the non-existence of Lake Jircacocha make a 1941-magnitude event appear unlikely, the steep glacierized mountain walls in the back of the lake may produce ice or rock-ice avalanches leading to impact waves, dam overtopping, erosion, and subsequent mass flows. Investigations by Klimeš et al. (2016) of the steep lateral moraines surrounding the lake indicate that failures and slides from moraines are possible at several sites but do not have the potential to create a major overtopping wave, partly due to the elongated shape of the lake. Rivas et al. (2015) elaborated on the possible effects of moraine-failure-induced impact waves. Recently, Somos-Valenzuela et al. (2016) used a combination of simulation approaches to assess the possible impact of process chains triggered by ice avalanching into Lake Palcacocha on Huaráz. They considered three scenarios of ice avalanches detaching from the slope of Palcaraju (0.5, 1.0, and 3.0 million m³) in order to create flood intensity maps and to indicate travel times of the mass flow to various points of interest. For the large scenario, the mass flow would reach the uppermost part of the city of Huaráz after approx. 1 h 20 min, for the other scenarios this time would increase to 2 h 50 min (medium scenario) and 8 h 40 min (small scenario). For the large scenario in particular, a high level of hazard is identified for a considerable zone near the Quilcay River, whereas zones of medium or low hazard become more abundant with the medium and small scenarios, or with the assumption of a lowered lake level (Somos-Valenzuela et al., 2016). In addition, Chisolm and McKinney (2018) analysed the dynamics impulse waves generated by avalanches using FLOW-3D. A similar modelling approach was applied by Frey et al. (2018) to derive a map of GLOF hazards for the Quilcay catchment. For Lake Palcacocha the same ice avalanche scenarios as applied by Somos-Valenzuela et al. (2016) were employed, with correspondingly comparable results in the Cojup Valley and for the city of Huaráz.

5.1 Scenario A – event induced by overtopping; fluid without yield strength

Outflow from Lake Palcacocha starts immediately, leading to (1) lowering of the lake level and (2) retrogressive erosion of the moraine dam. The bell-shaped fluid discharge curve at the hydrograph profile O1 (Fig. 4) reaches its peak of 18 700 m³ s⁻¹ after approx. 780 s and then decreases to a small residual (Fig. 5a). Channel incision happens quickly – 53 m of lowering of the terrain at the reference point R1 occurs within the first 1200 s, whereas the lowering at the end of the simulation is 60 m (Fig. 6a). This number represents an underestimation, compared to the reference value of 76 m (Table 2). The lake level decreases by 42 m, whereby 36.5 m of the decrease occurs within the first 1200 s. The slight underestimation, compared to the reference value of 47 m of lake level decrease, is most likely a consequence of uncertainties in the topographic reconstruction. A total amount of 1.5 million m³ is eroded from the moraine dam of Lake Palcacocha, corresponding to an underestimation of 22 %, compared to the reconstructed breach volume. Underestimations mainly occur at both sides of the lateral parts of the eroded channel near the moraine crest – an area where additional post-event erosion can be expected, so that the patterns and degree of underestimation appear plausible (Fig. 7a). In contrast, some overestimation of erosion occurs in the inner part of the dam. For numerical reasons, some minor erosion is also simulated away from the eroded channel. The iterative optimization procedure results in an entrainment coefficient $C_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{6.75}}$.

Figure 6Evolution of flow height and basal topography at the outlets of Lake Palcacocha (reference point R1 in Fig. 4b), and Lake Jircacocha (reference point R2 in Fig. 4c). The reference points are placed in a way to best represent the evolution of the breach in the dam for Lake Palcacocha, and the evolution of the impact wave for Lake Jircacocha. Additionally, the evolution of the lake level is shown for Lake Palcacocha. Note that the result for Scenario B is only displayed for Lake Palcacocha (e), whereas the result for Scenario C is only illustrated for Lake Jircacocha (f). The vertical distance displayed on the y axis refers to the terrain height or the lake level at the start of the simulation, respectively, whereby the flow height is imposed onto the topography. In Scenario B, the initial impact wave at the dam of Lake Palcacocha is only poorly represented due to the low temporal resolution of the simulation, and due to blurring by numerical effects (e).

Download

Figure 7Simulated versus reconstructed entrainment patterns for scenarios A and AX. The total entrained height and the difference between simulated and reconstructed entrainment (error) are shown. (a) Lake Palcacocha, Scenario A. (b) Lake Jircacocha, Scenario A. (c) Lake Palcacocha, Scenario AX. (d) Lake Jircacocha, Scenario AX.

The deposit of much of the solid material eroded from the moraine dam directly downstream from Lake Palcacocha, as observed in the field (Fig. 2c), is reasonably well reproduced by this simulation, so that the flow proceeding down-valley is dominated by the fluid phase (Fig. 8). It reaches Lake Jircacocha after t=840 s (Fig. 5b). As the inflow occurs smoothly, there is no impact wave in the strict sense, but it is rather the steadily rising water level (see Fig. 6b for the evolution of the water level at the reference point R2) inducing overtopping and erosion of the dam. This only starts gradually after some lag time, at approx. t=1200 s. The discharge curve at the profile O2 (Fig. 4) reaches its pronounced peak of 750 m³ s⁻¹ solid and 14 700 m³ s⁻¹ fluid material at t=2340 s and then tails off slowly.

Figure 8Evolution of the flow in space and time (Scenario A).

In the case of Lake Jircacocha, the simulated breach is clearly shifted south, compared to the observed breach. With the optimized value of the entrainment coefficient $C_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{7.15}}$, the breach volume is underestimated by 24 %, compared to the reconstruction (Fig. 7b). Also, here, this intentionally introduced discrepancy accounts for some post-event erosion. However, we note that volumes are uncertain as the reconstruction of the dam of Lake Jircacocha – in contrast to Lake Palcacocha – is a rough estimation due to lacking reference data.

Due to erosion of the dam of Lake Jircacocha, and also erosion of the valley bottom and slopes, the solid fraction of the flow increases considerably downstream. Much of the solid material, however, is deposited in the lateral parts of the flow channel, so that the flow arriving at Huaráz is fluid-dominated again (Fig. 8). The front enters the alluvial fan of Huaráz at t=2760 s, whereas the broad peak of 10 500 m³ s⁻¹ of fluid and 2000 m³ s⁻¹ of solid material (solid fraction of 16 %) is reached in the period between 3600 and 3780 s (Figs. 4 and 5c). Discharge decreases steadily afterwards. A total of 2.5 million m³ of solid and 14.0 million m³ of fluid material passes the hydrograph profile O3 until t=5400 s. Referring only to the solid, this is less material than reported by Kaser and Georges (2003). However, (i) there is still some material coming after, and (ii) pore volume has to be added to the solid volume, so that the order of magnitude of material delivered to Huaráz corresponds to the documentation in a better way. Still, the solid ratio of the hydrograph might represent an underestimation.

As prescribed by the parameter optimization, the volumes entrained along the channel are of the same order of magnitude, but lower than the reconstructed volumes summarized in Table 2: 0.7 million m³ of material is entrained upstream, 1.5 million m³ is entrained downstream of Lake Jircacocha, and 5.3 million m³ is entrained in the promontory. Figure 9a summarizes the travel times and the flow velocities of the entire process chain. Frontal velocities mostly vary between 5 and 20 m s⁻¹, with the higher values in the steeper part below Lake Jircacocha. The low and undefined velocities directly downstream of Lake Jircacocha reflect the time lag of substantial overtopping. The key numbers in terms of times, discharges, and volumes are summarized in Table 5.

Figure 9Travel times and frontal velocities for scenarios (a) A and (b) AX. Void fields in the profile graph refer to areas without clearly defined flow front.

Table 5Summary of the key results obtained with the computational experiments A–C. Refer to Tables 1 and 2 for the volumes involved, and to Table 4 for empirically expected peak discharges. Note that all entrained volumes are composed of 80 % of solid and 20 % of fluid material in terms of volume.

^* Peak of initial overtopping as response to the impact wave: 7000 m³ s⁻¹.

Download Print Version | Download XLSX

5.2 Scenario AX – event induced by overtopping; fluid with yield strength

Adding a yield strength of τ_y=5 Pa to the characteristics of the fluid substantially changes the temporal rather than the spatial evolution of the process cascade. As the fluid now behaves as fine mud instead of water and is more resistant to motion, velocities are lower, travel times are much longer, and the entrained volumes are smaller than in Scenario A (Fig. 9b; Table 5). The peak discharge at the outlet of Lake Palcacocha is reached at t=1800 s. Fluid peak discharge of 8200 m³ s⁻¹ is less than half the value yielded in Scenario A (Fig. 5d). The volume of material eroded from the dam is only slightly smaller than in Scenario A (1.4 versus 1.5 million m³). The numerically induced false positives with regard to erosion observed in Scenario A are not observed in Scenario AX, as the resistance to oscillations in the lake is higher with the added yield strength (Fig. 7c). Still, the major patterns of erosion and entrainment are the same. Interestingly, erosion is deeper in Scenario AX, reaching 76 m at the end of the simulation (Fig. 6c) and therefore the base of the entrainable material (Table 2). This is most likely a consequence of the spatially more concentrated flow and therefore higher erosion rates along the centre of the breach channel, with less lateral spreading than in Scenario A.

Consequently, Lake Jircacocha is also reached later than in Scenario A (Fig. 6d), and the peak discharge at its outlet is delayed (t=4320 s) and lower (7600 m³ s⁻¹ of fluid and 320 m³ s⁻¹ of solid material) (Fig. 5e). A total of 2.0 million m³ of material is entrained from the dam of Lake Jircacocha, with similar spatial patterns to those in Scenario A (Fig. 7d). Huaráz is reached after t=4200 s, and the peak discharge of 5000 m³ s⁻¹ of fluid and 640 m³ s⁻¹ of solid material at O3 occurs after t=6480 s (Fig. 5f). This corresponds to a solid ratio of 11 %. Interpretation of the solid ratio requires care here as the fluid is defined as fine mud, so that the water content is much lower than the remaining 89 %. The volumes entrained along the flow channel are similar in magnitude to those obtained in the simulation of Scenario A (Table 5).

5.3 Scenario B – event induced by impact wave

Scenario B is based on the assumption of an impact wave from a 3 million m³ landslide. However, due to the relatively gently sloped glacier tongue heading into Lake Palcacocha at the time of the 1941 event (Figs. 2a and 4b), only a small fraction of the initial landslide volume reaches the lake, and impact velocities and energies are reduced, compared to a direct impact from the steep slope. Approx. 1 million m³ of the landslide have entered the lake by t=120 s, an amount which only slightly increases thereafter. Most of the landslide deposits on the glacier surface. Caused by the impact wave, discharge at the outlet of Lake Palcacocha (O1) sets in at t=95 s and, due to overtopping of the impact wave, immediately reaches a relatively moderate first peak of 7000 m³ s⁻¹ of fluid discharge. The main peak of 16 900 m³ s⁻¹ of fluid and 2000 m³ s⁻¹ of solid discharge occurs at t=1200 s due do the erosion of the breach channel. Afterwards, discharge decreases relatively quickly to a low base level (Fig. 10a). The optimized value of $C_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{6.75}}$ is used also for this scenario. The depth of erosion along the main path of the breach channel is clearly less than in Scenario A (Fig. 6e). However, Table 5 shows a higher volume of eroded dam material than the other scenarios. These two contradicting patterns are explained by Fig. 11a: the overtopping due to the impact wave initiates erosion of not only the main breach, but also of a secondary breach farther north. Consequently, discharge is split among the two breaches and therefore less concentrated, explaining the lower erosion at the main channel despite a larger total amount of eroded material. The secondary drainage channel can also be deduced from observations (Fig. 3a) but has probably played a less important role than suggested by this simulation.

Figure 10Hydrographs of moraine dam failure of Lake Palcacocha (a, d), landslide dam failure of Lake Jircacocha (b, e), and the flow entering the urban area of Huaráz (c, f) for scenarios B and C. Note that, for clarity, fluid flow heights and discharges are plotted in negative direction.

Download

Figure 11Simulated versus reconstructed entrainment patterns for scenarios B and C. The total entrained height and the difference between simulated and reconstructed entrainment (error) are shown. (a) Lake Palcacocha, Scenario B. (b) Lake Jircacocha, Scenario B. (c) Lake Jircacocha, Scenario C.

The downstream results of Scenario B largely correspond to the results of Scenario A, with some delay partly related to the time from the initial landslide to the overtopping of the impact wave. Discharge at the outlet of Lake Jircacocha peaks at t=2700 s (Fig. 10b), and the alluvial fan of Huaráz is reached after 3060 s (Fig. 10c). The peak discharges at O2 and O3 are similar to those obtained in Scenario A. The erosion patterns at the dam of Lake Jircacocha (again, $C_{\mathrm{E}}={\mathrm{10}}^{-\mathrm{7.15}}$) very much resemble those yielded with scenarios A and AX (Fig. 11b), and so does the volume of entrained dam material (2.2 million m³). The same is true for the 2.5 million m³ of solid and 13.9 million m³ of fluid material entering the area of Huaráz until t=5400 s, according to this simulation.

Also in this scenario, the volumes entrained along the flow channel are very similar to those obtained in the simulation of Scenario A. The travel times and frontal velocities – resembling the patterns obtained in Scenario A, with the exception of the delay – are shown in Fig. 12a, whereas Table 5 summarizes the key numbers in terms of times, volumes, and discharges.

Figure 12Travel times and frontal velocities for scenarios (a) B and (b) C. Note that the legend of (a) also applies to (b). Void fields in the profile graph refer to areas without clearly defined flow front.

5.4 Scenario C – event induced by dam collapse

In Scenario C, we assume that the breached part of the moraine dam collapses, and the collapsed mass mixes with the water from the suddenly draining lake and flows downstream. The more sudden and powerful release, compared to the two other scenarios, leads to higher frontal velocities and shorter travel times (Fig. 12b; Table 5).

In contrast to the other scenarios, impact downstream starts earlier, as more material is released at once, instead of steadily increasing retrogressive erosion and lowering of the lake level. The fluid discharge at O1 peaks at almost 40 000 m³ s⁻¹ (Fig. 10d) rapidly after release. Consequently, Lake Jircacocha is reached already after 720 s, and the impact wave in the lake evolves more quickly than in all the other scenarios considered (Fig. 6f). The lake drains with a peak discharge of 15 400 m³ s⁻¹ of fluid and 830 m³ s⁻¹ of solid material after 1680–1740 s (Fig. 10e). In contrast to the more rapid evolution of the process chain, discharge magnitudes are largely comparable to those obtained with the other scenarios. The same is true for the hydrograph profile O3: the flow reaches the alluvial fan of Huaráz after t=2160 s, with a peak discharge slightly exceeding 10 000 m³ s⁻¹ of fluid and 2000 m³ s⁻¹ of solid material between t=2940 and 3240 s. A total of 2.7 million m³ of solid and 14.6 million m³ of fluid material enter the area of Huaráz until t=5400 s, which is slightly more than in the other scenarios, indicating the more powerful dynamics of the flow (Table 5). The fraction of solid material arriving at Huaráz remains low, with 16 % solid at peak discharge and 15 % in total. Again, the volumes entrained along the flow channel are very similar to those obtained with the simulations of the other scenarios (Table 5).

6.1 Possible trigger of the GLOF process chain

In contrast to other GLOF process chains in the Cordillera Blanca, such as the 2010 event at Laguna 513 (Schneider et al., 2014), which was clearly triggered by an ice-rock avalanche into the lake, there is disagreement upon the trigger of the 1941 multi-lake outburst flood in the Quilcay catchment. Whereas, according to contemporary reports, there is no evidence of a landslide (for example, ice avalanche) impact on the lake (Vilímek et al., 2005; Wegner, 2014), and dam rupture would have been triggered by internal erosion, some authors postulate at least a small impact starting the process chain (Portocarrero, 2014; Somos-Valenzuela et al., 2016).

Each of the three assumed initiation mechanisms of the 1941 event, represented by scenarios A/AX, B, and C, yields results which are plausible in principle. We consider a combination of all three mechanisms a likely cause of this extreme process chain. Overtopping of the moraine dam, possibly related to a minor impact wave, leads to the best correspondence of the model results with the observation, documentation, and reconstruction. In particular, the signs of minor erosion of the moraine dam north of the main breach (Fig. 3a) support this conclusion: a major impact wave, resulting in violent overtopping of the entire frontal part of the moraine dam, would supposedly also have led to more pronounced erosion in that area, as to some extent predicted by Scenario B. There is also no evidence for strong landslide–glacier interactions (massive entrainment of ice or even detachment of the glacier tongue), which would be likely scenarios in the case of a very large landslide. Anyway, the observations do not allow for substantial conclusions on the volume of a hypothetic triggering landslide: as suggested by Scenario B, even a large landslide from the slopes of Palcaraju or Pucaranra could have been partly alleviated on the rather gently sloped glacier tongue between the likely release area and Lake Palcacocha.

The minor erosional feature north of the main breach was already visible in the photo of Kinzl (Fig. 2a), possibly indicating an earlier, small GLOF. It remains unclear whether it was reactivated in 1941. Such a reactivation could only be directly explained by an impact wave, but not by retrogressive erosion only (A/AX) or internal failure of the dam (C) – so, more research is needed here. The source area of a possible impacting landslide could have been the slopes of Palcaraju or Pucaranra (Fig. 1), or the calving glacier front (Fig. 2a). Attempts to quantify the most likely release volume and material composition would be considered speculative due to the remaining difficulties in adequately simulating landslide–lake and possibly landslide–glacier–lake interactions (Westoby et al., 2014). Further research is necessary in this direction. In any case, a poor stability condition of the dam (factor of safety ∼1) could have facilitated the major retrogressive erosion of the main breach. A better understanding of the hydro-mechanical load applied by a possible overtopping wave and the mechanical strength of the moraine dam could help to resolve this issue.

The downstream patterns of the flow are largely similar for each of scenarios A, AX, B, and C, with the exception of travel times and velocities. Interaction with Lake Jircacocha disguises much of the signal of process initiation. Lag times between the impact of the flow front on Lake Jircacocha and the onset of substantial overtopping and erosion are approx. 10 in in scenarios A and B, and less than 3 min in Scenario C. This clearly reflects the slow and steady onset of those flows generated through retrogressive erosion. The moderate initial overtopping in Scenario B seems to be alleviated before reaching Lake Jircacocha. Sudden mechanical failure of the dam (Scenario C), in contrast, leads to a more sudden evolution of the flow, with more immediate downstream consequences.

6.2 Parameter uncertainties

We have tried to back-calculate the 1941 event in a way reasonably corresponding to the observation, documentation, and reconstruction, and building on physically plausible parameter sets. Earlier work on the Huascarán landslides of 1962 and 1970 has demonstrated that empirically adequate back-calculations are not necessarily plausible with regard to parameterization (Mergili et al., 2018b). This issue may be connected to equifinality issues (Beven, 1996; Beven and Freer, 2001), and in the case of the very extreme and complex Huascarán 1970 event, by the inability of the flow model and its numerical solution to adequately reproduce some of the process components (Mergili et al., 2018b). In the present work, however, reasonable levels of empirical adequacy and physical plausibility are achieved. Open questions remain with regard to the spatial differentiation of the basal friction angle required to obtain adequate results (Table 3): lower values of δ downstream from the dam of Lake Jircacocha are necessary to ensure that a certain fraction of solid passes the hydrograph profile O3 and reaches Huaráz. Still, solid fractions at O3 appear rather low in all simulations. A better understanding of the interplay between friction, drag, virtual mass, entrainment, deposition, and phase separation could help to resolve this issue (Pudasaini and Fischer, 2016a, b; Pudasaini, 2019a, b).

The empirically adequate reproduction of the documented spatial patterns is only one part of the story (Mergili et al., 2018a). The dynamic flow characteristics (velocities, travel times, hydrographs) are commonly much less well documented, particularly for events in remote areas which happened a long time ago. Therefore, direct references for evaluating the empirical adequacy of the dimension of time in the simulation results are lacking. However, travel times play a crucial role related to the planning and design of (early) warning systems and risk reduction measures (Hofflinger et al., 2019). Comparison of the results of scenarios A and AX (Fig. 9) reveals almost doubling travel times when adding a yield stress to the fluid fraction. In both scenarios, the travel times to Huaráz are of the same order of magnitude as the travel times simulated by Somos-Valenzuela et al. (2016) and therefore considered plausible, so that it is hard to decide which is the more adequate assumption. Even though the strategy of using the results of earlier simulations as reference may increase the robustness of model results, it might also reproduce errors and inaccuracies of earlier simulation attempts, and thereby confirm wrong results.

The large amount of more or less pure lake water would point towards Scenario A, whereas intense mixing and entrainment of fine material would favour Scenario AX. More work is necessary in this direction, also considering possible phase transformations (Pudasaini and Krautblatter, 2014). At the same time, the optimization and evaluation of the simulated discharges remains a challenge. Here we rely on empirical relationships gained from the analysis of comparable events (Walder and O'Connor, 1997).

6.3 Implications for predictive simulations

Considering what was said above, the findings from the back-calculation of the 1941 event can help us to better understand and constrain possible mechanisms of this extreme process chain. In principle, such an understanding can be transferred to present hazardous situations in order to inform the design of technical remediation measures. Earlier, measures were implemented not only at Lake Palcacocha (Portocarrero, 2014), but also at various other lakes such as Laguna 513: a tunnelling scheme implemented in the 1990s strongly reduced the impacts of the 2010 GLOF process chain (Reynolds, 1998; Reynolds et al., 1998; Schneider et al., 2014).

However, the findings of this study should only be applied for forward simulations in the same area or other areas with utmost care. The initial conditions and model parameters are not necessarily valid for events of different characteristics and magnitudes (Mergili et al., 2018b). In the case of Lake Palcacocha, the situation has changed substantially since 1941: the lake level is much lower and the volume larger, and the lake is directly connected to the steep glacierized slopes, so that the impact of a hypothetic landslide could be very different now. Also, the current lake is dammed by a different moraine than the pre-1941 lake, with a very different dam geometry (Somos-Valenzuela et al., 2016). In general, the mechanisms of the landslide impact into the lake, which were not the focus of the present study, would require more detailed investigations. Ideally, such work would be based on a three-phase model (Pudasaini and Mergili, 2019; considering ice as a separate phase) and consider knowledge and experience gained from comparable, well-documented events. A possible candidate for such an event would be the 2010 event at Laguna 513, which was back-calculated by Schneider et al. (2014). In general, it remains a challenge to reliably predict the outcomes of given future scenarios. The magnitude of the 1941 event was amplified by the interaction with Lake Jircacocha, whereas the 2012 GLOF process chain in the Santa Cruz Valley (Mergili et al., 2018a) alleviated due to the interaction with Lake Jatuncocha, comparable in size. While it seems clear that the result of such an interaction depends on event magnitude, topography, and the dam characteristics of the impacted lake, Mergili et al. (2018a, b) have demonstrated the high sensitivity of the behaviour of the simulated flow to the friction parameters, but also to the material involved (release mass, entrainment). A larger number of back-calculated process chains will be necessary to derive guiding parameter sets which could facilitate predictive simulations, and so will an appropriate consideration of model uncertainties and possible threshold effects (Mergili et al., 2018b). Earlier studies, considering the 2010 event at Laguna 513 (Schneider et al., 2014) and three future scenarios for Lake Palcacocha (Somos-Valenzuela et al., 2016) have followed a different strategy, using model cascades instead of integrated simulations, so that a comparison with studies based on r.avaflow is only possible to a limited extent.

Another remaining issue is the lateral spreading of the flow on the fan of Huaráz, which is overestimated in all four simulations (Figs. 8, 9, and 12): the most likely reason for this is the insufficient representation of fine-scale structures such as buildings or walls in the DEM, which would serve as obstacles confining the flow in lateral direction.