Modelling future lahars controlled by different volcanic eruption scenarios at Cotopaxi (Ecuador) calibrated with the massively destructive 1877 lahar

Lahars are among the most hazardous mass flow processes on earth and have caused up to 23 000 casualties in single events in the recent past. The Cotopaxi volcano, 60 km southeast of Quito, has a well‐documented history of massively destructive lahars and is a hotspot for future lahars due to (i) its ~10 km2 glacier cap, (ii) its 117–147‐year return period of (Sub)‐Plinian eruptions, and (iii) the densely populated potential inundation zones (300 000 inhabitants). Previous mechanical lahar models often do not (i) capture the steep initial lahar trajectory, (ii) reproduce multiple flow paths including bifurcation and confluence, and (iii) generate appropriate key parameters like flow speed and pressure at the base as a measure of erosion capacity. Here, we back‐calculate the well‐documented 1877 lahar using the RAMMS debris flow model with an implemented entrainment algorithm, covering the entire lahar path from the volcano edifice to an extent of ~70 km from the source. To evaluate the sensitivity and to constrain the model input range, we systematically explore input parameter values, especially the Voellmy–Salm friction coefficients μ and ξ. Objective selection of the most likely parameter combinations enables a realistic and robust lahar hazard representation. Detailed historic records for flow height, flow velocity, peak discharge, travel time and inundation limits match best with a very low Coulomb‐type friction μ (0.0025–0.005) and a high turbulent friction ξ (1000–1400 m/s2). Finally, we apply the calibrated model to future eruption scenarios (Volcanic Explosivity Index = 2–3, 3–4, >4) at Cotopaxi and accordingly scaled lahars. For the first time, we anticipate a potential volume growth of 50–400% due to lahar erosivity on steep volcano flanks. Here we develop a generic Voellmy–Salm approach across different scales of high‐magnitude lahars and show how it can be used to anticipate future syneruptive lahars.

117-147-year return period of (Sub)-Plinian eruptions, and (iii) the densely populated potential inundation zones (300 000 inhabitants). Previous mechanical lahar models often do not (i) capture the steep initial lahar trajectory, (ii) reproduce multiple flow paths including bifurcation and confluence, and (iii) generate appropriate key parameters like flow speed and pressure at the base as a measure of erosion capacity. Here, we back-calculate the well-documented 1877 lahar using the RAMMS debris flow model with an implemented entrainment algorithm, covering the entire lahar path from the volcano edifice to an extent of $70 km from the source. To evaluate the sensitivity and to constrain the model input range, we systematically explore input parameter values, especially the Voellmy-Salm friction coefficients μ and ξ. Objective selection of the most likely parameter combinations enables a realistic and robust lahar hazard representation. Detailed historic records for flow height, flow velocity, peak discharge, travel time and inundation limits match best with a very low Coulomb-type friction μ (0.0025-0.005) and a high turbulent friction ξ (1000-1400 m/s 2 ). Finally, we apply the calibrated model to future eruption scenarios (Volcanic Explosivity Index = 2-3, 3-4, >4) at Cotopaxi and accordingly scaled lahars. For the first time, we anticipate a potential volume growth of 50-400% due to lahar erosivity on steep volcano flanks. Here we develop a generic Voellmy-Salm approach across different scales of high-magnitude lahars and show how it can be used to anticipate future syneruptive lahars. The Indonesian term 'lahar' refers to a rapidly flowing mixture of sediment and water that originates from a volcano (Smith & Fritz, 1989), either directly related to eruptive activity (i.e. primary/syneruptive lahar) or during post-eruptive periods (i.e. secondary lahar). These mass flows rank among the most disastrous volcanic natural hazards (Auker et al., 2013). The destructive potential of lahars arises from their sudden onset and outstanding flow characteristics (featured by high sediment-carrying capacity), speeds exceeding tens of kilometres per hour, discharges up to about 50 000 m 3 /s, and runout distances of tens to even hundreds of kilometres (Mothes et al., 2004;Thouret et al., 2020). About 56 000 lahar-related fatalities have been documented since 1500 AD (Brown et al., 2017), including about 23 000 fatalities at Nevado del Ruíz in Colombia in 1985 (Pierson et al., 1990). As population density and infrastructure development close to active volcanoes is rapidly growing, an increasing number of people are potentially exposed to volcanic hazards (Auker et al., 2013;Small & Naumann, 2001). In the last few decades, efforts have been made to reduce the risk to human life by implementing effective monitoring and early-warning systems, and by assessing probable lahar inundation areas through modelling and reconstruction of the eruptive history. The active, glacier-capped Cotopaxi volcano is among the most dangerous volcanoes worldwide, because it can produce very large eruption-triggered lahars, which may inundate densely inhabited valleys close to the volcano. In this study, we aim to contribute to the hazard assessment of syneruptive lahars by testing the applicability of a Voellmy-based debris flow model including entrainment and analysing its capability to reproduce key characteristics of lahars.
Direct collection of quantitative data on lahar events is inherently difficult, and especially rare for large-scale syneruptive flows. Understanding of the principal physical and rheological characteristics of lahars benefits from observations of experimental studies (Iverson et al., 2010) and natural debris flows (McArdell et al., 2007;McCoy et al., 2012). Key features include a pulsating progression (Doyle et al., 2011;Iverson, 1997;McArdell et al., 2007), long-lived pore-fluid pressure in excess of hydrostatic values contributing to long runout distances (Iverson, 1997;Lube et al., 2009;McArdell et al., 2007), and topographically controlled erosion and entrainment leading to a several-fold increase of discharge and flow volume after initiation (Berger et al., 2010;Berti et al., 1999;Manville, 2004;McCoy et al., 2012;Pierson et al., 1990;Scott et al., 2005).

| Lahar modelling approaches
Incorporating the wide range of flow behaviours and the dynamic processes characteristic of lahars in runout analyses poses a challenging task, so that any theoretical model or practical simulation tool necessarily draws on physical and rheological simplifications of varying degrees. Different computer codes simulating lahars have been proposed, including the widely used empirical model LAHARZ (Iverson et al., 1998), which predicts inundation limits (Delaite et al., 2005;Huggel et al., 2008;Lee et al., 2015;Muñoz-Salinas et al., 2009;Pistolesi et al., 2014;Schneider et al., 2008;Worni et al., 2011), the cellular-automaton model LLUNPIY (Lupiano et al., 2020;Machado et al., 2015), and a variety of physically based numerical models that apply different rheology assumptions. Among these, 1-D models are today largely supplanted by multiscale continuum models with depthaveraged equations for conservation of mass and momentum. Among rheological laws, fluid dynamics-based models that propagate lahars as Newtonian flows are frequently applied. Such models (e.g. the dynamic flood-routing code Delft3D) can well reproduce dilute lakebreakout lahars from Mount Ruapehu, NZ (Carrivick et al., 2008).
However, modelling lahars as Newtonian fluids neglects that lahars differ substantially in energy dissipation compared to typical Newtonian flood flows. Fluid-dynamic models may therefore be extended by the Manning's or Chézy's coefficient to approximate internal flow resistance. The hydraulic simulation software FLO-2D (O'Brien et al., 1993), designed for floods, mudflows and debris flows, uses the viscoplastic Bingham rheology model and encompasses flow resistance terms for yield, viscous, turbulent and dispersive stresses. Flo-2D is one of the most widely used single-phase rheological models in lahar hazard analysis (Caballero and Capra, 2014;Worni et al., 2011), whereas other rheological approaches-such as the Voellmy modelhave only been rarely applied, though it has been formulated for (and proven to be successful in) simulating various types of sediment-laden flows. More sophisticated, 3-D models based on the Coulomb mixture theory (Iverson & Denlinger, 2001) describe a two-phase flow and the interactions between the solid and fluid phase. A simplified depthaveraged formulation of the mixture theory is incorporated in the simulation code Titian2D (Pitman & Le, 2005), and has been applied in a few studies for sediment-laden lahars (Procter et al., 2010;Williams et al., 2008).
With increasing model complexity and physical accuracy, a growing number of input parameters need to be defined, requiring appropriate material sampling and testing (Iverson & George, 2016).
Complex models are also computationally expensive, restricting runout modelling to a specific, small section of the flow path (Procter et al., 2010). When modelling large-magnitude lahars over long distances, it is necessary to strike a balance between robustness and physical plausibility: computationally tractable, robust models are critical for application to large-scale drainage systems, while also capturing key flow features like stage heights, flow velocity, basal shear pressure, travel times, confluence, bifurcation and run-ups. Therefore, the numerical simulation software RAMMS for debris flow is chosen in this study, as it provides a good combination of simplification and physical plausibility and allows the tracing of sediment-laden lahars at Cotopaxi along the entire trajectory in large catchments.

| Erosion in lahar models
Entrainment of eroded material into the flowing mass can have a significant impact on flow propagation (Iverson, 1997;McDougall & Hungr, 2005) and should therefore be considered in runout modelling and hazard assessment (Frank et al., 2015;Iverson & Ouyang, 2015). Although there is no common opinion on whether entrainment of substrate material reduces or enhances flow mobility (Iverson et al., 2011;Mangeney, 2011;Pudasaini & Fischer, 2016), an increase in flow volume has been observed to positively correlate with peak discharge, potential inundation area, runout distance, flow height and velocity (Rickenmann, 1999). The number of numerical models capable of simulating entrainment of mass flows is currently growing (Iverson & Ouyang, 2015), and the performance of the underlying theories is increasingly tested and validated against real mass flow events (Dietrich & Krautblatter, 2019;Frank et al., 2017;Hungr & McDougall, 2009;Hussin et al., 2012).
The simulation of entrainment processes either requires the input of user-specific growth rates or of process-based erosion rates as a function of velocity (Fagents & Baloga, 2006) or shear stress (Frank et al., 2015;Iverson, 2012). Although lahars develop a considerable erosive force related to flow height, flow density and shear stress, lahar models have only rarely addressed erosional processes or volumetric changes (Carrivick, 2007;Carrivick et al., 2010;Fagents & Baloga, 2006;Lupiano et al., 2020), largely due to the shortcomings of the models used. These models are forced to start the calculation already with the hypothetical total flow volume, either on the volcano edifice or in channelized topography further downstream. This approach may not only lead to an overprediction of flow discharge and stage height, but also means that the dynamic evolution of syneruptive lahars, characterized by initial flow acceleration, rapid volumetric growth on the steep volcano flanks, and further erosion and deposition during downstream propagation, cannot be traced in the modelling process. Previous simulations for lahars from Cotopaxi using 1-D or 2-D flood-routing models (Aguilera et al., 2004;Barberi et al., 1992;Toapaxi et al., 2019;Vera et al., 2019)

| Aims of this study
This study aims to systematically assess the hazards of massive syneruptive lahars at Cotopaxi volcano in the densely inhabited northern and southern drainage system, by back-analysing the destructive 1877 lahar event, objectively constraining best-fitting model inputs and employing these for potential future lahars. Therefore, we (i) test the applicability of a Voellmy-Salm model for massive syneruptive lahars, (ii) develop a generic model simulating the lahar from the initiation on steep volcano flanks to the distal reaches, and (iii) include erosivity and entrainment on the steep slopes of a volcano edifice and along the trajectory in the generic model.

| Geological setting and volcanic history
Cotopaxi is a large stratovolcano located in the Ecuadorian Andes of South America, $60 km southeast of Quito, Ecuador's capital ( Figure 1). The volcanic edifice stands on the metamorphic belt corresponding to the Cordillera Real, but is also partially interbedded with volcanic formations from the neighbouring Chalupas Caldera (Hall & Mothes, 2008). The glacier cap of the Cotopaxi extends from its summit (5897 m a.s.l.) down all the flanks, reaching $4900 m a.s.l. on the eastern flank and $5300 m a.s.l. on the western flank.
Similar to other tropical glaciers of South America, Cotopaxi's glacier cap has been retreating rapidly for several decades (Cáceres, 2017;Vuille et al., 2018).
The geological evolution of Cotopaxi started with an ancient complex named Cotopaxi I between 560 and 420 ka BP, followed by a long period of quiescence (Hall & Mothes, 2008). The volcano resumed its activity around 13 ka BP with the formation of the current edifice Cotopaxi II, which is characterized by andesitic and rhyolitic eruptions. Around 4.5 ka BP, the last rhyolitic eruption took place and produced the partial collapse of the northern volcano flank, with subsequent formation of a debris avalanche and a gigantic cohesive lahar (Mothes et al., 1998).
The last 2100 years at Cotopaxi II have witnessed repetitive andesitic eruptions characterized by the deposition of lava flows, regional tephra layers and large-scale syneruptive lahars (Pistolesi et al., 2011(Pistolesi et al., , 2013. Average recurrence periods between 117 and 147 years may be estimated for Cotopaxi eruptions (Barberi et al., 1995), though it has been shown that the occurrence of eruptions displays significant clustering (Pistolesi et al., 2011).
Primary lahars formed during explosive eruptions undoubtedly represent the most hazardous phenomena at the Cotopaxi Pistolesi et al., 2013;Sierra et al., 2019). These flows are formed during explosive eruptions, when pyroclastic density currents move on top of the volcano glacier and melt it, instantaneously producing large volumes of water quickly mixing with volcanic rocks (Pistolesi et al., 2013). These flows transit through the three main drainage systems descending from the volcano The most recent volcanic activity at Cotopaxi was recorded in April 2015 and identified by the monitoring system (Hidalgo et al., 2018). The anomalous records increased steadily until mid-August, when hydromagmatic activity started, producing significant ash venting from the summit crater Gaunt et al., 2016). This activity decreased progressively, waning in December 2015 (Hidalgo et al., 2018), and did not produce syneruptive lahars.

| The 1877 lahar disaster
The last syneruptive lahar from Cotopaxi occurred on 26 June 1877.
According to reports by Sodiro (1877) and Wolf (1878), Cotopaxi had been in eruption since early 1877, but in June its activity increased significantly. That day, the volcano produced a sustained explosive phase with widespread pyroclastic density currents descending all the edifice's flanks. Estimations by these authors suggest that the explosive eruption lasted at least 15 min, and that it took around 1 h for the lahar front to arrive in Latacunga (south) and the Chillos Valley (close to Quito to the north).
In the Pita drainage, the 1877 lahar experienced a significant runover, diverting $20% of its volume into the Rio Santa Clara (Mothes et al., 2004), which transits densely populated areas in the Chillos Valley. The pirating effect occurred in the site called 'La Caldera', where the Rio Pita is deeply incised and experiences a sharp bend, while the headwaters of the Rio Santa Clara are placed just above. Calculations by Mothes et al. (2004), assuming that the pirating effect was produced by overflowing resulting from super-elevation in a forced vertex, suggest that the discharge rate was >40 000 m 3 /s, but they also suggest that the flow could have been partially dammed by the very narrow canyon immediately downstream.
The lahar deposit from 1877 can still be distinguished in the landscape of the proximal north and north-eastern lower flanks of Cotopaxi: it contains a specific type of scoriaceous volcanic bomb with a unique feature of fracturing and disaggregating easily when hit lightly with a hammer (Mothes et al., 2004), which was not repeated in volcanic bombs from 18th-century lahar deposits. This enables detailed mapping of the deposits, especially in the proximal zones (<20 km from the source), as well as measurement of precise cross-sections from which physical parameters like height, velocity and discharge rate may be obtained (Aguilera et al., 2004;Mothes et al., 2004). Total lahar volumes are estimated to be about 60-75 × 10 6 m 3 in the northern drainage and about 80-100 × 10 6 m 3 in the southern drainage, respectively (Mothes et al., 2004).
The 1877 lahar was influenced by the morphology of the transited channels (Mothes et al., 2004). In the northern and southern  (Mothes et al., 2004;Pistolesi et al., 2013). Deposits from the 1877 lahar are granulometrically similar to other deposits from Cotopaxi eruptions in the 18th century and lack significant amounts of very fine-grained, clay-sized particles (Mothes et al., 2004), as is characteristic of non-cohesive flows.

| RAMMS model
We apply the physically based dynamic model RAMMS, which is widely used in runout modelling of different types of mass movements including snow avalanches, rock avalanches, debris flows and GLOFS (Frey et al., 2018;Hussin et al., 2012;Schneider et al., 2010).
The RAMMS debris flow model is formulated as a single-phase system that propagates the solid-fluid mixture as an incompressible bulk flow with mean constant density (Bartelt et al., 2017;Christen et al., 2010).
It requires careful calibration of the two governing frictional parameters of the Voellmy-Salm rheology and delivers information on key parameters for hazard assessment, such as runout distance, height and velocity of the flow, impact pressure, erosion and deposition. The core of the program is a second-order numerical solution scheme to solve the depth-averaged equations of motion in 3-D terrain (x,y,z) at time t . The mass balance is given by where _ Q x, y, t ð Þ denotes the mass production source term, H(x,y,t) the flow height and U(x,y,t) the depth-averaged velocity. The depthaveraged momentum balance in the x and y directions is given by and ∂ t HU y ð Þ+ ∂ y HU 2 y + g z k a=p respectively, where g = (g x , g y , g z ) describes the gravitational vector, k a/p is the earth pressure coefficient (normally set to 1), S g and S f denote the slope parallel gravitational acceleration and frictional deceleration, respectively.
The frictional resistance S f is described using a Voellmy-Salm approach (Salm, 1993;Voellmy, 1955), which incorporates a dry Coulomb friction μ proportional to the normal stress and a velocitysquared drag or viscous-turbulent friction ξ (similar to Chézy friction) (Bartelt et al., 1999;Christen et al., 2010;Fischer et al., 2012): where ρ is the density of the mass, φ the slope angle of the surface, and κ the terrain curvature in the flow direction. The Coulomb friction μ accounts for the resistance of the solid phase, is independent of velocity, and governs flow behaviour when moving slowly, while the turbulent friction ξ accounts for the fluid phase, scales with flow velocity, and dominates at higher velocities in acceleration zones (Bartelt et al., 2017;Christen et al., 2010).
An implemented entrainment module in RAMMS predicts the depth of erosion in terms of net decrease in elevation of the channel bed (Frank et al., 2017). It simulates the increase in flow volume due to entrainment (Frank et al., 2017), but sediment entrainment is not coupled to topographic changes in the underlying digital elevation model (DEM) or changes in flow behaviour. Field observations at the exceptionally active Illgraben catchment in Switzerland on erosion depths (Schürch et al., 2011) and erosion rates (Berger et al., 2011) of naturally occurring debris flows served as the basis for the development of the entrainment algorithm (Frank et al., 2015). The potential erosion depth e m per grid cell is described as a function of the computed shear stress acting on the channel bed, given by and is related to a critical shear stress τ c and a proportionality factor dz dτ controlling the rate of vertical erosion, given by Based on average erosion rates recorded at the Illgraben (Berger et al., 2011), a specific erosion rate dz dt is used (Frank et al., 2015):

| Model setup
We aim to set up the RAMMS lahar model in a generic way by starting the flow on the volcano flank close to the crater, enabling sediment bulking along the path and propagating the lahar until distal reaches.
The setup is designed for the purpose of back-calculating the 1877 lahar event and can be applied to forward modelling of potential future lahars.

| Model premises: DEM and limits of model domain
Lahars from Cotopaxi usually do not show coherent deposit limits but may travel up to 325 km, as reported for the 1877 lahar (Mothes et al., 2004), distally in a watery phase. Therefore, the modelling region is restricted by the largest urban centres around the volcano: to the north, the model region extends until approximately 70 km downstream of the volcano, covering the urban agglomeration of Valle de Los Chillos; to the south, it ends south of the city of Salcedo about 65 km downstream of Cotopaxi (see Figure 2). Given the large extent, a spatial resolution of 10 m is considered sufficient for capturing key flow patterns. All numerical simulations conducted in this study are based on a 10 m DEM, which we resampled from an original 4 m-resolution DEM developed by the SIGTIERRAS project in 2010 (Ministerio de Agricultura, Ganadería, Acuacultura y Pesca). No detailed information is available on valley morphology before the 1877 lahar passage. However, we can assume that the present-day topography of the partly deeply incised major channels does not differ significantly, because historical and stratigraphic evidence suggests that they confined all known historic lahars since the 18th century (Hall & Mothes, 2008).

| Model premises: Erosion and entrainment
We are aware of the high uncertainties involved with modelling lahar erosivity, as no quantitative information exists about the erosion potential of historic lahars at Cotopaxi. Alternatively, by neglecting dynamic bulking and debulking, it would be necessary to already start the model with the total lahar volume, which may lead to an overestimation of peak discharge and flow height (Frank et al., 2015). Our approach to simulate sediment entrainment along the trajectory is therefore (i) focused on modelling the overall mass balance and reproducing averaged total flow volumes reported by Mothes et al. (2004) and (ii)  abundant loose, pyroclastic sediment or glacial deposits exposed by recent glacier retreat are present, and (3) areas with limited erodibility characterized by regolith-covered bedrock, or deep bedrockflanked channels.
For each class, we specify the RAMMS entrainment parameters that control the degree of sediment entrainment into the flow (Table 1) from a set of given values (Bartelt et al., 2017;Frank et al., 2015Frank et al., , 2017. In areas with 'limited erodibility', the maximum erosion depth is constrained to 0.5 m modulated by assumed regolith thickness, whereas areas with 'average' and 'high erodibility' are considered supply-unlimited. For a rough approximation of the release hydrograph, we use an idealized three-point discharge hydrograph. The shape parameters draw upon the following empirical relation and observations:

| Estimation of initial volume and input hydrograph
1. Peak discharge. Mizuyama et al. (1992) proposed an empirical formula to estimate the peak discharge Q p of granular debris flows as a function of volume V, given by Pierson (1998) tested the applicability of this correlation for large-magnitude lahars and found that it serves as an adequate approximation for near-source peak discharges. The formula is applied to the tested initial lahar volumes ranging from 25 × 10 6 to 58 × 10 6 m 3 .
2. Time of peak discharge. Well-studied lahars at the glacier-capped volcanoes of Mount St. Helens and Nevado del Ruiz showed a relatively instantaneous release of flow (Pierson, 1995;Pierson & Scott, 1985;Pierson et al., 1990). The time of peak discharge is therefore set to 120 s after initiation.
3. Duration. We can conservatively consider eyewitness accounts of pyroclastic flow duration as a proxy for the duration of lahar release. Wolf (1878) reported that pyroclastic flows lasted for 15-20 min and caused immediate melting of snow and ice all over the crater. Application of the above-described volume-discharge correlation by Mizuyama et al. (1992) to a triangular-shaped hydrograph yields release durations of about 12 min, which is largely in accordance with eyewitness observations.
One broad release hydrograph is placed at the upstream boundary of each of the two river systems, near the present location of the glacier snouts. This enables the injected flow to drain completely into one system, and all major channels of the respective drainage system are affected simultaneously.

| Model calibration and sensitivity analysis
We back-calculate the well-documented 1877 lahar event in the northern drainage system of Cotopaxi in order to evaluate the performance of the model and to derive plausible model parameter values for future lahar scenarios. All available information on the reference event is based on eyewitness reports by Sodiro (1877)   (1) Average erodibility 0.025 0.1 1.0 - The calibration criteria act as simulation constraints and are used to determine best-fit input parameter combinations. Because input parameters may cause a significant change in the model's behaviour, we systematically analyse the sensitivity of the two Voellmy-Salm friction coefficients μ and ξ, the flow density, and the initial velocity to model outcomes for runout distance, peak discharge, travel time, flow velocity and flow height (represented by the calibration criteria).
Given the location of the release hydrographs at least 1000 m downslope of the crater, the lahar most likely has already started to evolve due to melting of glacier ice and mixing with freshly ejected volcanic products upon 'arriving' at the location of the release hydrograph. It seems likely that the initial flow has already accelerated to a certain velocity, which is approximated by exploring the initial velocity

| Evaluation of model performance
In order to objectively evaluate the quality of model predictions with RAMMS, we separately compare model outcomes (m i ) to target values (t i ) for calibration measures (i) of flow velocity, flow height, travel time and peak discharge. The performance of a simulation (ϕ) using a given set of input parameters (y) is quantified as percentage deviation, with Normalizing the differences between modelled and targeted values facilitates an easy comparison between the different rheological measures of the calibration criteria.
Similarity between simulated and mapped lahar impact area (cf. Figure 3)  For each set of input parameters (y), the values of the constraintwise performance analyses are then summed to give a single simulation performance value (ϕ total ), with This number highlights minimum model deviations within the analysed input parameter space and thereby allows us to select the most representative parameter combinations for the back-calculated event.

| Simulation of potential future lahars
Based on investigations about the eruptive history of Cotopaxi (Hall & Mothes, 2008;Mothes et al., 1998Mothes et al., , 2004 We adapt some of these simplified assumptions to the quantification of future lahar scenarios in this study by using the following approach: 1. The calibrated lahar model of the 1877 eruption with estimated release discharge and input volume (defining the input hydrograph) provides the basis for scenario modelling. The modelled event was equivalent to a VEI 3-4 eruption.
2. Due to a 50% reduction in glacier surface area, a future VEI 3-4 eruption is expected to produce lahars with release discharges half those of the 1877 lahar (assuming that the total glacier areas in 1877 and 1977 were of the same extent). This clearly represents a simplified assumption; the degree to which glacier retreat affects future lahar magnitudes may also need a closer look at eruption dynamics and triggering mechanisms (Pistolesi et al., 2013(Pistolesi et al., , 2014, as discussed later. This corresponds to a volume growth of 50 and 80% in the northern and southern drainage system, respectively. At both river systems, we can observe that the amount of eroded material varies with initial volume, revealing a general trend towards smaller erosion-attributed volume growth with increasing input volume. This trend is more pronounced in the southern drainage, where erosion and entrainment results in a volume growth ranging from 135% for an initial volume of 25 × 10 6 m 3 to 70% for an initial volume of 58 × 10 6 m 3 .

| Simulation of potential future lahars
We applied the calibrated lahar model with previously derived bestfitting input parameters (μ = 0.005 and ξ = 1400 m/s 2 ) and estimates of initial volume and discharge of the 1877 lahar to potential future lahars during three different eruption scenarios at Cotopaxi. Our approach to assess future lahar volumes considering the effect of glacier retreat yielded the following hydrograph estimates (Figure 9): for Scenarios 2, 3 and 4 in the northern drainage, we calculate peak discharges of respectively 32 000, 64 000 and 128 000 m 3 /s and corresponding input volumes of 7.8, 18.9 and 46 × 10 6 m 3 . In the southern drainage, we derive for Scenarios 2, 3 and 4 peak discharges of respectively 33 600, 67 300 and 134 500 m 3 /s and corresponding input volumes of 8.3, 20.2 and 49 × 10 6 m 3 .
The results of the scenario simulations are summarized in Table 2 and Figure 10; more detailed simulation results of flow heights, flow speed and erosion are provided in the online supporting information.
We observed that simulated lahars in the northern drainage reach total flow volumes of $20 × 10 6 m 3 in Scenario 2, $35 × 10 6 m 3 in Scenario 3, and $70 × 10 6 m 3 in Scenario 4, which corresponds to a volume growth due to entrainment of 150, 80 and 50%, respectively.  F I G U R E 9 (a) Potential future eruption scenarios at Cotopaxi (Andrade et al., 2005), with corresponding input hydrographs for lahar modelling in the northern (b) and southern (c

| Estimation of erosion capacity
To estimate erosion and entrainment by lahars, we used a forward modelling approach, which is based on estimating the erodibility of the substrate and is independent of the magnitude of an eruptive event. We demonstrated that lahars from Cotopaxi may experience substantial bulking during downstream propagation (for details, see Table 2).
Our simulations imply ( Figure 11)  clad volcanoes evolve: rapidly transitioning from an initially more dilute flow to a hyperconcentrated flow or debris flow due to coincident entrainment of pyroclastic material in proximal regions (Major & Newhall, 1989;Pierson et al., 1990).  (Table 1). The impact of topography on modelled erosive forces is also evident in the significant reduction in eroded volume as the lahars enter low-gradient areas at the foot of the volcano (see Figure 11: after 7-8 and 9-13 min in the northern and southern drainage, respectively).
In contrast to previous models at Cotopaxi, this paper calculates erosivity along the entire initial steep and the flatter channelized trajectory, whereas previous models have only anticipated the latter. To compare erosion in low-gradient reaches, we estimate that a lahar of similar magnitude like in 1877 ( Figure 11 and Table 2 Ruíz lahars (Pierson et al., 1990). Using the same empirical relation, Pistolesi et al. (2014) found for the southern 1877 lahar a bulking value of 8.5 × 10 6 m 3 (7% of assumed total volume of 120 × 10 6 m 3 ).
In terms of quantifying total eroded volumes of future lahars from Cotopaxi, it was found that they may erode and entrain between $12 × 10 6 and $43 × 10 6 m 3 of substrate material, and can potentially grow in volume by 50-400% compared to initial volume over a modelled travel distance of $70 km ( (ii) the triggering mechanism (VEI = 3 eruption), which generated pyroclastic flows interacting with the glacier cap, and (iii) the runout distance (70-100 km from the source). Moreover, (iv) recent models for lahars from Cotopaxi (Mothes et al., 2004;Ordóñez et al., 2013;Vera et al., 2019) use a bulking factor of 3 as a rule of thumb throughout all modelled eruption scenarios, which is the mean of the bulking factors of 2 to 4 reported for Nevado del Ruíz lahars (Pierson et al., 1990 (Pierson et al., 1990). This is about twice as much eroded sediment compared to our simulations at Cotopaxi, where a lahar with an initial volume of $20 × 10 6 m 3 (Scenario 3, southern drainage) may erode and entrain $36 × 10 6 m 3 . Second, we found that volumetric growth of lahars may vary remarkably among the different future eruption scenarios: in the southern drainage, initial volumes of small lahars (Scenario 2) increase by $400%, whereas lahars with large initial volumes (Scenario 4) grow by $100%. This indicates that even small initial lahars from Cotopaxi are possibly capable of entraining large volumes of sediment, favoured by the presence of easily F I G U R E 1 1 Eroded volume by lahars in the northern (blue) and southern (red) drainage as a function of simulation time (output in 60 s intervals). The northern and southern lahar started with an initial volume of 46 × 10 6 and 49 × 10 6 m 3 , respectively, and simulations of both lahars are performed with μ = 0.005 and ξ = 1400 m/s 2 . The release hydrograph in the southern sector is placed higher on the volcano flank, increasing the distance over which the flow overruns loose and easily erodible volcaniclastic material [Colour figure can be viewed at wileyonlinelibrary.com] erodible pyroclastic sediment on the volcano flanks with high erosion rates and low yield stresses ( and entrainment by large-magnitude lahars, as these erosion parameters were derived from natural debris flows in the Alps (Frank et al., 2015). The loose, pyroclastic material present on the volcano flanks, which largely controls the volumetric growth of the lahar, was characterized by erosion rates twice as high as the average observed rates in the Illgraben catchment (Frank et al., 2017), but possibly the entrainment is even more extreme. Application to lahar events with documented spatial erosion patterns is needed to quantify possible deviations. We also note that the model is limited to predict material entrainment from vertical erosion in the channel bed, but does not consider additional volume sources such as channel bank collapse, (re-) activated major landslides along the channels and in source areas on the volcano, or freshly ejected volcaniclastic material, which may attain a thickness in proximal (<10 km) zones of 1-1.5 m during a future VEI > 4 eruption (Andrade et al., 2005).

| Model performance for lahar simulation
We back-calculated the 1877 lahar event at Cotopaxi and performed a multi-criteria calibration using calibration measures such as impact area, flow height, flow velocity, travel times and peak discharge (Mothes et al., 2004;Sodiro, 1877;Wolf, 1878 (Hussin et al., 2012;Schraml et al., 2015), GLOFs (Frey et al., 2018;Schneider et al., 2014) or ice-rock avalanches (Schneider et al., 2010;Sosio et al., 2008). However, Sosio et al. (2011) used similar μ/ξ combinations (μ between 0.001 and 0.02, ξ between 700 and 1200 m/s 2 ) to simulate volcanic debris avalanches evolving from sudden failure of volcano flanks. An undisputable reason for our calibration result is that only small μ/large ξ combinations can reproduce the extremely long runout distances of lahars from Cotopaxi. A small value of μ and a larger value of ξ indicate low flow resistance (Equation (4)) and lead to a slow stopping mechanism, which in turn is the reason why the model fails to reproduce material deposition in the plains around the volcano. Large values of ξ further indicate that the turbulent friction term in Equation (4) is marginally relevant for the simulation of syneruptive lahars in this study.

| Simulation of future Cotopaxi lahars
Predicting the magnitude of future syneruptive lahars is an inherently difficult challenge, because it involves estimating the efficacy of volcano-glacier interactions during an eruption and consequent water supply, as well as the potential of material incorporation and volume growth of the lahar (Pierson, 1995). To anticipate the latter, we proposed a generic model approach and start the lahar on the steep upper volcano slopes in order to provide a first estimate of erosion capacity of lahars. As discussed previously, the scenario simulations revealed a variable degree of volume growth among the three potential future scenarios and the two drainages. This model approach can, therefore, support the estimate of volcano-specific bulking factors of lahars.
In our model, the efficacy of volcano-glacier interactions is quantified in the release hydrograph of each eruption scenario, defining the volume of slurry produced by rapid mixing of hot eruption products with ice and snow. Similar to previous lahar studies at Cotopaxi (Aguilera et al., 2004;Barberi et al., 1992;Mothes et al., 2004), our estimations of initial volumes of future lahars are directly related to the glacier extent. Because the extent of Cotopaxi's glaciers has diminished by at least $50% since the last major eruption in 1877 (Cáceres, 2017), we assumed that this would convert into a 50% reduction in initial volumes of future lahars. However, the degree to which glacier retreat influences resulting lahar volumes is an open issue. Pistolesi et al. (2013Pistolesi et al. ( , 2014 argued that lahar magnitudes expected at Cotopaxi strongly depend on eruption dynamics and the type of lahar triggering mechanism: 1877-type eruptions (VEI 3-4) with boiling-over activity mostly produce scoria flows that erode deep, linear canyons into the glacier. Glacier extent could therefore have only a minor effect on initial lahar volumes. Consequently, future eruptions will presumably generate lahars similar in magnitude to that of 1877, despite a 50% reduction in glacier surface area. In contrast, VEI > 4 eruptions produce pyroclastic surges or flows that interact with the entire glacier cap, and glacier retreat would thus exert a major control on initial lahar volumes. In a nutshell, our study shows that Voellmy-Salm models can effectively and robustly anticipate future lahars, but it remains essential to communicate uncertainties involved in estimating future lahar magnitudes and consequently simulated lahar impact areas to decision-making authorities and downstream communities.

| CONCLUSION
In this study, we simulated massive long-distance lahars from Coto- Lahar erosion is estimated by a forward model approach, that relies on the erodibility of the underlying substrate given by the current geological and topographical setting. The simulations demonstrate that the degree of volume growth decreases with release volume, ranging between 50% for large (46 × 10 6 m 3 ) and 400% for comparably smaller (8 × 10 6 m 3 ) initial volumes, respectively. Despite the uncertainties involved in this approach, the qualitative prediction of simulated lahar erosivity on the steep volcano flanks using the proposed straightforward estimation of erosion capacity is promising compared with observed spatial erosion patterns. Further information will be gained by transferring this approach to more recent lahar events, where quantitative and detailed spatial information about erosion by flow is readily available.
A reliable estimation of the magnitudes of future eruptiontriggered lahars and according hazards to downstream communities is challenging and needs to consider the degree of volcano-glacier interactions as well as the degree of solids incorporation into the flow. The generic model approach developed in this study includes lahar erosivity and can help to anticipate the volumetric growth and runout patterns of syneruptive lahars. However, critical processes determining initial lahar magnitudes (i.e. eruption magnitude, eruption dynamics, interactions between volcanic products and ice/snow, impact of glacier retreat) are not thoroughly understood and need more dedicated research.