From examination of natural events to a proposal for risk mitigation of lahars by a cellular-automata methodology: a case study for Vascún valley, Ecuador

Lahars are erosive floods, mixtures of water and pyroclastic detritus, known for being the biggest environmental disaster and causing a large number of fatalities in volcanic areas. Safety measures have been recently adopted in the threatened territories by constructing retaining dams and embankments in key positions. More disastrous events could be generated by the difficulty of maintaining these works in efficiency and for the changed risk conditions originating from their presence and the effects of their functioning. LLUNPIY/3r, a version of the cellular-automaton model LLUNPIY for lahar simulations, is presented. The growing frequency of lahars in the Vascún valley of Tungurahua volcano (Ecuador), probably due to the effects of the climatic change, has recently produced smaller and less dangerous events, sometimes favored by the collapse of ponds generated by small landslides. An investigation is performed here in order to reproduce such situations in a controlled way by the use of LLUNPIY/3r simulations. Using precise field data, points are individuated where dams by backfill, which are easy to collapse, can produce the formation of ponds; LLUNPIY/3r simulations permit projecting the triggering of small lahars by minor rainfall events or projecting, in the case of more rainfall, the anticipation of lahar detachment, avoiding simultaneous and dangerous confluence with other lahars.

1 Introduction

1.1 The problem of the lahars

Lahars are one of the most devastating phenomena in terms of the number of fatalities they cause in volcanic areas (Neall, 1976; Waythomas, 2014). They are flows other than the common stream flow and consist of pyroclastic deposits mixed with water. Their physical properties (density, viscosity, and consistency) are very similar to wet concrete that has not yet hardened (Vallance, 2000). This fluid, under steep-slope conditions, is capable of reaching speeds up to 100 km h⁻¹ and distances up to 300 km, and it becomes solid when water is gradually released in flat areas (Manville et al., 2013).

Lahars may be of a primary type (or syneruptive) if directly related to volcanic eruptions; usually when glaciers and/or snow are melted by pyroclastic or lava flows, they develop from mixing pyroclastic material with water, as in the tremendous 1985 Colombian event of Nevado del Ruiz (Pierson et al., 1990). Another case could occur when a large quantity of water is available from the breakout of a natural lake because of eruption (Manville, 2010).

Secondary lahars are produced when a large water quantity is available directly by extreme meteorological events or indirectly by the overflow of superficial water bodies. Pyroclastic deposits of previous eruptive activities are mobilized, e.g., the pyroclastic flows of Mt. Pinatubo in 1991 in the Philippines (Rodolfo et al., 1996).

Soil erosion with water inclusion along streams increases the volume of both primary and secondary lahars.

Two main triggering mechanisms are possible.

• a.

  The mobilization process is related to pyroclastic sediments sometimes mixed with some exotic material (tephra): if the superficial water amount exceeds a water height threshold, related to features of pyroclastic stratum and soil slope, then the percolation can cause a detachment in the unconsolidated stratum.

• b.

  The erosion process mainly depends on the redistribution of volcanic sediment along the slopes (Leavesley et al., 1989; Major et al., 2000): unconsolidated tephra, swept away by flows, mixes with water and gradually enlarges its volume because of the contribution of both sediments and water (Barclay et al., 2007).

Different approaches were considered in literature for lahar modeling: empirical models (e.g., Schilling, 1998; Muñoz-Salinas et al., 2009) were developed, accounting mainly for some macro-observable phenomena, simplified hydrological and rheological models were developed that reduce the lahar behavior to a Newton-like behavior (e.g., O'brien, 1993; Costa, 1997), numerical methods were developed for approximating partial differential equations (PDEs; e.g., Pitman et al., 2003), and the cellular-automata (CA) alternative methodology was developed for lahar modeling, which will be shown later.

In many issues regarding complex systems, research was able to progress thanks to computer simulations, which allowed for developing multidisciplinary and transdisciplinary approaches, linked in part to the emergence of alternative computing paradigms, such as CA (Toffoli, 1984; Chopard and Droz, 1998; Iovine et al., 2007).

1.2 Multicomponent or macroscopic cellular automata

CA is both a parallel computational paradigm and an archetype for modeling “complex dynamical systems”, which are extended into space and can be described on the basis of local interactions of their constituent parts. A homogeneous CA can be seen as a d-dimensional space, partitioned in cells of uniform size, each one embedding an identical input–output computing device (a finite-state automaton). Input for each cell is given by the states of the neighboring cells, where the neighborhood conditions are determined by a time- and space-invariant pattern. At the time t=0, cells are in arbitrary states (initial conditions) and the CA evolves by changing them simultaneously at discrete times (CA step), according to the transition function σ: S^m→S, where S is the finite set of the states and m is the number of the neighboring cells (Di Gregorio and Serra, 1999).

A short exemplification is given by the CA majority: a two-dimensional space is divided into square cells so that the neighborhood of one square is given by the single element itself together with the eight surrounding it. Their states are blue (0) and red (1) and are added together within the neighborhood by the transition function. If the sum is more than 4 (the majority of neighbors are red), the next state of the cell will be red; otherwise it will be blue. Sometimes the system evolves from an initial distribution of reds and blues in a complex way, causing expansions of red and blue colors from some areas (Toffoli, 1984).

When complex macroscopic dynamical systems, like “surface flow” phenomena (lahars, debris flows, snow avalanche, lava flows, and pyroclastic flows), are modeled, the previous definitions are insufficient, and the following extensions were adopted for the multicomponent or macroscopic CA (MCA).

The abstract CA must be related univocally to the real phenomenon in its dynamics, and each cell has to correspond to a portion of the space or surface (of the territory T) where the phenomenon evolves; thus the time corresponding to a step of the transition function has to be set, and the size of the cell has to be specified, e.g., by the length of its edge. These constant values in time and space are called global parameters. P is the set of global parameters, and it includes both physical and empirical parameters. The choice of some parameters is imposed, where possible, by the precision that must be achieved in the simulation, e.g., cell dimension. The value of some other parameters is deduced instead by the physical features of the phenomenon, e.g., the parameter related to energy dissipation by turbulence. In these cases, the initial physically sounding value considered at the beginning of validation is corrected by attempts in the phase of model validation on the base of the simulation quality, depending on the comparison between the real event and simulation results. A methodology based on genetic algorithms was usually used for calibrating the parameters of our CA models of surface flows (Iovine et al., 2005).

Each characteristic, relevant to the evolution of the system and relative to the spatial portion corresponding to the cell, is individuated as a substate. The finite set Q of the states is given by the Cartesian product of the substates: $Q=Q_{\mathrm{1}}\times Q_{\mathrm{2}}\times \mathrm{\dots }\times Q_n$. Examples of a lahar model substates related to the part of territory corresponding to a cell are the average altitude (substate altitude), the thickness of the lahar (substate lahar thickness), and the depth of the erodible (unconsolidated) pyroclastic stratum (substate pyroclastic stratum depth). The dynamics of the phenomenon is expressed by the variation in the substate values in the successive steps of simulations for each cell. Considering that the features related to the third dimension may be expressed in terms of substates, it is possible to develop two-dimensional models operating three-dimensionally in reality (Avolio et al., 2012).

MCA has to account for phenomena whose dynamics involves more interacting processes, sometimes of different nature, e.g., loss of lahar energy because of erosion of the unconsolidated pyroclastic stratum of the “cell” and loss of energy of the lahar in the cell caused by its turbulence. These interacting processes compose the transition function and are called “elementary” processes of the CA. They are computed sequentially, involving the update of the MCA substates.

The last extension of MCA is the “external influences”, which account for input from the “external world” that does not depend on local interactions (that cannot be reduced to local interactions) occurring at some cells of the CA; e.g., the external influence “lava alimentation at the vents” is applied, at each step, only to those cells corresponding to the locations where the vents actually are. The value of the substate “lava quantity” is updated by adding the amount that is considered to be discharged (in the case of simulation of a real event) or that is supposed to be discharged (in the case of simulation of a conjectured event) in the cell during the time step (Di Gregorio and Serra, 1999).

Simulations of flow-like landslides have been performed by several versions of the MCA model SCIDDICA since 1987 for both subaerial and subaqueous debris, granular, and mudflows (e.g., Barca et al., 1987; Avolio et al., 2008, 2013; Mazzanti et al., 2010; Lupiano et al., 2014, 2015a, b, c, 2017). Simulations of primary and secondary lahars were performed by the MCA model LLUNPIY (Machado et al., 2014, 2015a, b; Chidichimo et al., 2016).

LLUNPIY, SCIDDICA-SS3, and SCIDDICA-SS2 are our most advanced models (in the sense that they include the features of the previous models plus other new ones) for simulating flow-like landslides and lahars. Unlike other models, which were used in lahar simulations (LAHARZ – e.g., Schilling, 1998; Muñoz-Salinas et al., 2009; TITAN2D – e.g., Sheridan et al., 2005; Williams et al., 2008; Córdoba et al., 2015), they allow the implementation of the erosion process.

1.3 Strategies of risk mitigation for lahars

Reliable simulation tools are very important in order to develop risk mitigation strategies and to test them in different conditions. Such instruments have to be used with extreme caution because the complex problem of the lahar hazard must be studied with an interdisciplinary approach (Lane et al., 2003; Leung et al., 2003); e.g., mitigation measures which involve engineered protection structures could modify hazard conditions in time and could increase the disaster risk, as better specified below.

Besides tools of early warning, which could work only partially, and aside from temporary or definitive land evacuation, which could involve a strong social impact and economic destitution, security measures have been adopted in volcanic territory by constructing retaining dams, embankments, walls, dykes, levees, and reservoirs in key positions for containing and deviating possible lahars (Scott, 1989; Verstappen, 1992; Aguilera et al., 2004; Künzler et al., 2012; Carey et al., 2012). This solution could involve a strong environmental impact: it is difficult to guarantee the constant efficiency of these works, and their presence, together with the effects of their functioning, could severely change the risk conditions (Janda et al., 1981, 1996; Scott, 1989; Procter et al., 2010; Shreve and Kelman, 2014; Wisner et al., 2012). A short paper of Kelman (2007) underlines the following:

Despite decades of evidence from research and practice demonstrating that reliance on structural approaches increases disaster risk over the long-term, structural approaches are frequently preferred without properly considering complementary or alternative measures. Examples of structural approaches are walls, dams, dykes, levees, and reservoirs. While they do provide some benefits, decisions to implement them and nothing else are usually made by emphasizing the short-term benefits and discounting the long-term costs.

The growing frequency of lahars in the area of the Vascún valley of the Tungurahua volcano, Ecuador, has recently produced smaller (shorter accumulation periods) and therefore less dangerous events (Mothes and Vallance, 2015). Moreover, small landslides, forming natural dams with temporary ponds, could easily trigger lahars by collapsing because of rainfall. This sometimes happens; e.g., the IGEPN (Instituto Geofísico Escuela Politécnica Nacional, Quito, Ecuador) reported such a case on 23 August 2008 (IGEPN, 2008a, b). These extraordinary combinations of events resulted in the idea of using the overabundant pyroclastic material, available on-site, to create easy-to-collapse artificial dams. The dam breakdown is obtained through the appropriate sizing of the cross section of the structure which is designed to fail at the achievement of a specified water level. This goal is reached through the implementation of an ad hoc numerical model, based on the finite element method (FEM), for the stability analysis of the dam slopes.

For such points, minor rainfall events can produce small lahars, and major events will anticipate the lahar detachment, avoiding simultaneous confluence with other lahars. The control of the dams collapse could permit, in various situations, many combinations for a controlled triggering of lahars in order to mitigate the risk.

The next chapter is devoted to the geological description of the Vascún valley, and the third chapter introduces LLUNPIY/3r, the MCA model for simulating primary and secondary lahars, together with its validation in simulating some significant lahars of Vascún valley. In the fourth chapter, if the building of dams that are easy to collapse is considered, then their favorable locations are hypothesized for a controlled triggering of lahars, the effects of which are simulated for possible events. At the end of the paper are conclusions and comments.

2 Geological setting of Vascún valley

Tungurahua is one of the most active and dangerous volcanoes in the Ecuadorian Andes (Cordillera Oriental) on the inside of Sangay National Park; its summit, at 5023 m a.s.l., is positioned at 01^∘28^′ S latitude and 78^∘27^′ W longitude. It is a stratovolcano, whose evolution has involved the succession of three major volcanic edifices (Tungurahua I, Tungurahua II, and Tungurahua III) since the Middle Pleistocene over a basement of metamorphic rocks. Historical eruptions have all originated from the summit crater. The main events of eruptive activity occurred during the following years: 1640–1641, 1773–1777, 1886–1888, 1916–1918, and from 1999 to the present (Hall et al., 1999; Mothes et al., 2015; Biggs et al., 2010). The average frequency of eruptions in the last 2000 years is once per century, according to the detailed studies of Le Pennec et al. (2008). They have been accompanied by strong explosions and sometimes by pyroclastic flows and lava flows that reached populated areas at the volcano's base. The formation of rain-induced lahars is also a cause of danger. Approximately 32 000 people live within the higher-risk areas, mainly in rural villages and in the touristic (thermal springs) town of Baños de Agua Santa (Mothes and Vallance, 2015). Baños de Agua Santa (1800 m a.s.l.) is only 8 km from the summit; 30 km to the NW, 30 km to the SW and 140 km to the N are the distances to the summit, respectively, from the towns of Ambato and Riobamba and from the capital Quito.

The small glacier of Tungurahua volcano is reduced both by the global warming phenomenon and by the intensification of the volcanic activity after 1999; therefore today's snow cover is completely negligible in terms of primary lahar generation. Lahar triggering, with significant frequency and magnitude, is hence subordinate to both to the intensity and duration of the rainfall and the available quantity of fresh material (cumulated unconsolidated pyroclastic material from volcanic eruption) along the slopes and within the principal canyons (Quebradas of the Vascún River, Juive Grade – La Pampa valley, and Achupashal Quebrada) and, to a lesser extent, to other factors. Almost all the lahars are confined to the canyons and converge into the Pastaza River (Mothes and Vallance, 2015).

The eruptive activity of Tungurahua volcano, during the last years, has generated a greater availability of pyroclastic debris that is periodically remobilized from atmospheric phenomena, often not particularly violent but prolonged for several days. Between 2000 and 2011, around 900 rain-induced lahars were triggered (Mothes and Vallance, 2015).

Generally, lahar magnitude is small and, consequently, causes limited damage. Precious data were supplied by the IGEPN, Quito, Ecuador, and its survey stations around the country (https://www.igepn.edu.ec/, last access: 2 December 2019); it is the ultimate authority for various volcanic hazards and earthquakes. IGEPN and its Acoustic Flow Monitor (AFM) station detect the passage of the majority of secondary lahars (e.g.: IGEPN, 2005, 2008a, b), while many others are traced by the Observatory of the Volcano Tungurahua (OVT), which is situated 13 km to the north–northwest of the crater, also with the observation contribution of local volunteers (vigias).

The 2005 and 2008 secondary lahars of Vascún valley

Vascún River gave its name to the valley on the NE flank of Tungurahua volcano; it flows into the Pastaza River (Fig. 1). The valley slopes exceed 35–40^∘ in the first (steepest) 3 km, while they range from 20 to 6^∘ in the last (less steep) 2–3 km. Furthermore, the path of Vascún River is extremely sinuous in the upper 1–2 km, where the river is characterized by a succession of tight 90^∘ bends.

Figure 1 Tungurahua volcano. The Vascún valley is inside the red box.

The Vascún valley is highly susceptible to lahar flows that threaten the nearby populated areas which were inundated several times in the past years: the town of Baños, which extends into part of the depositional area, was affected several times; the thermal structure “El Salado”, which is located along the river banks, had been repeatedly damaged by passage of lahars (Fig. 2b).

Figure 2 (a) © Google Earth view of Tungurahua volcano with indication of main localities of the study area; (b) 2005 lahar at thermal structure “El Salado”.

For LLUNPIY/3r calibration and validation, two events were taken into consideration, which took place, respectively, in 2005 and 2008. The dynamic of the events is briefly described below.

On 12 February 2005, heavy rainfall caused the remobilization of ash fall deposits, generating lahars in the Vascún River valley. The mean velocity, estimated on the base of data recorded by the alert instrumentation, varied between 3 and 7 m s⁻¹ (Williams et al., 2008) according to the morphological characteristics of valley that crossed the sector. The lahar volume, measured by the staff of the AFM station, was estimated at approximately 70 000 m³, while a subsequent investigation, carried out by IGEPN (2005), estimated it to be 55 000 m³. The lahar crossed the valley for about 10 km, flooding the El Salado baths during the passage, and then reached the Pastaza River. The chronicle of the event is taken from the work of Williams et al. (2008), and we refer to their data and simulations for comparison with our simulation.

The 2008 lahar had a different dynamic. On August 13, a small landslide produced a natural dam along Vascún River at an elevation of about 2200 m a.s.l. The dam originated a pond with a length of 100 m, a depth of 3 m, and a width of 20 m. After heavy rains on 22 August, the dam collapsed and generated a lahar. The flow velocity was estimated in 15 m s⁻¹, with a flow rate of 1120 m³ s⁻¹ and an average height of 4 m (IGEPN, 2008a, b); this flow rate is 10 times greater than that recorded in the 2005 event (IGEPN, 2005; Williams et al., 2008). The lahar reached El Salado in 5 min and devastated the pools of the thermal spa; afterwards it destroyed some houses in Las Ilusiones (a village in the Baños district) further downstream.

3 LLUNPIY/3r model for lahar simulation

LLUNPIY (Lahar modelling by Local rules based on an UNderlying PIck of Yoked processes; llunp'iy means flood in the Quechua language) is a model for simulating secondary and primary lahars according to MCA methodology applied to a complex system, whose evolution may be mainly specified in terms of local interaction. MCA features of SCIDDICA-SS3 (Avolio et al., 2013) and SCIDDICA-SS2 (Avolio et al., 2008; Lupiano et al., 2016, 2017) are inherited by LLUNPIY; LLUNPIY for secondary lahars is extensively defined in Machado et al. (2015b). Here only the features of the model, which were applied in the study cases, are reported (reduced version LLUNPIY/3r from SCIDDICA-SS2). No external influences were considered. The LLUNPIY/3r simulation starts considering data related to the altitude (value of the “altitude” substate; see Sect. 1.2), to the depth of the erodible pyroclastic stratum (value of the “soil of the cell” substate; see Sect. 1.2), and to the lahar thickness (value of the “thickness of the lahar” substate; see Sect. 1.2) for each cell.

A reliable reconstruction of the first phase of a real event of a lahar permits setting an “initial” moment, where it is possible to deduce the thickness of lahar in the territory; these data constitute the values of the substate thickness of the lahar in the first step of the simulation. In the case of the simulation of a lahar produced by the collapse of a dam holding a given water volume, the thickness of lahar is deduced by the mixing of pond water with the dam material and part of the unconsolidated pyroclastic stratum below. Note that the simulated lahar events, which occurred in the Vascún valley, do not involve the very first phase of water percolation and detachment subsequent to water inclusion (Machado, 2015), since the collapse of the temporary pond is abrupt. In the cases of past events, data permitted simulation of the phenomenon just in the phase of lahar. Furthermore, in the simulation of real and hypothesized events, all the lahars end in the Pastaza River, so the last phase of lahar deposition is omitted and the viscosity of the lahar may be considered constant for these particular cases.

3.1 Introduction to the LLUNPIY/3r version

The following quintuple defines the two-dimensional (with hexagonal cells) MCA model LLUNPIY/3r:

$$<T,P,N,Q,\mathit{\sigma }>,$$

where $T=\mathit{\left\{}\right(x,y\left)\right|x,y\in \mathbb{N}$, $\mathrm{0}\le x\le l_x$, $\mathrm{0}\le y\le l_y\mathit{\}}$ is the set of hexagonal cells which tessellate the territory where the phenomenon evolves; the cells are individuated by the points with integer coordinates (Fig. 3, left) of their centers. ℕ is the set of the natural numbers.

Figure 3 (a) The neighborhood of cell (5, 6); (b) neighborhood indices.

P is the set of both the empirical and physical global parameters. They are related to the global common features of the phenomenon (Table 1):

Table 1 Physical and empirical parameters.

$N=<(\mathrm{1},-\mathrm{1})$, $(\mathrm{0},-\mathrm{1})$, $(-\mathrm{1},\mathrm{0})$, $(-\mathrm{1},\mathrm{1})$, (0,1), (1,0), $(\mathrm{0},\mathrm{0})>$, the neighborhood, identifies the geometrical pattern of cells (Fig. 3a), which influences the state change of the “central” cell. Index 0 is assigned to the central cell, and indices $\mathrm{1},\mathrm{\dots },\mathrm{6}$ are assigned to the six adjacent cells (Fig. 3b). N=7; the sum of indices of opposite directions is always 7.

Q is the finite set of states of the finite-state automaton incorporated into the cell; it is specified, in terms of substates, as their Cartesian product (Table 2).

Table 2 Substates.

σ: Q^#N→Q is the deterministic transition function for each cell in T. The following elementary processes compose σ, and they account for the overall dynamics of the phenomenon:

• τ_mob – effects of mobilization,

• τ_lo – lahar outflows,

• τ_te – effect of turbulence,

• τ_fc – composition of flows.

Physical parameters concern physical quantities that are used in equations of the transition function and correspond to values adopted in the implementation of the model (e.g., cell apothem p_a that depends on several factors, data precision, and insuperable approximation limits related to specific features of the phenomenon) or values as the temporal correspondence of a CA step p_t, which must account for the fact that $p_{\mathrm{a}}/p_{\mathrm{t}}>v_{\mathrm{mx}}$, where v_mx is the maximum possible velocity of flows during the development of the phenomenon. In other words the shift of a flow in a CA step may not exceed the neighborhood value. All the parameters, except p_a and p_t, are empirically set in the phase of model validation by the simulation quality, and initial values of parameters were deduced by the physical features of the phenomenon; e.g., a very slow lahar would emerge unbelievably in a simulation by the largest values of p_cf, the coefficient of friction, and p_dt, the energy dissipation due to turbulence.

3.2 The elementary processes of LLUNPIY/3r

In the following, an outline of the transition function, with a description of the elementary processes updating the substates, will be provided. The complete execution of all the elementary processes accomplishes a step of the LLUNPIY/3r. The neighborhood index in square brackets, following substate specification, indicates the corresponding cell of the neighborhood. ΔQ_S indicates variation in the substate Q_S. $Q_S^{\prime }$ indicates the new value of the substate Q_S, $Q_S^{\prime }=Q_S+\mathrm{\Delta }{Q}_S$. In the case of external and internal flows, the cell to which the flow is directed is specified in the substate, inserting a subscript which precedes it, e.g., ₂Q_E[1], or the external flow of the neighbor with index 1 toward its neighbor with index 2.

3.2.1 Pyroclastic cover mobilization

Soil features together with the quantity of water content determine a value p_tm of the mobilization threshold to be compared with the kinetic head Q_KH associated to lahar debris Q_LT; when Q_KH>p_tm, the pyroclastic cover Q_D is eroded, and the lahar thickness Q_LT augments and altitude Q_A diminishes according to the following empirical formula that turned out to be valid in different models of debris flow (e.g., Avolio et al., 2008), snow avalanches (Avolio et al., 2017), and primary and secondary lahars (e.g., Machado, 2015):

$$\begin{array}{}\text{(1)}& -\mathrm{\Delta }{Q}_{\mathrm{D}}=\mathrm{\Delta }{Q}_{\mathrm{LT}}=-\mathrm{\Delta }{Q}_{\mathrm{A}}=\left(Q_{\mathrm{KH}}-p_{\mathrm{tm}}\right)\cdot p_{\mathrm{pe}}.\end{array}$$

There is, correspondingly, a dissipation of energy proportional to the depth of erosion; it is specified by a decrease in kinetic head Q_KH according to the following formula:

$$\begin{array}{}\text{(2)}& -\mathrm{\Delta }{Q}_{\mathrm{KH}}=\left(Q_{\mathrm{KH}}-p_{\mathrm{tm}}\right)\cdot p_{\mathrm{de}}.\end{array}$$

3.2.2 Effect of turbulence

A loss of kinetic head occurs by turbulence at each LLUNPIY/3r step according to the following equation:

$$\begin{array}{}\text{(3)}& -\mathrm{\Delta }{Q}_{\mathrm{KH}}=p_{\mathrm{dt}}{Q}_{\mathrm{KH}},\end{array}$$

where p_dt is an empirical parameter that accounts for the turbulence kinetic energy (Launder and Spalding, 1983). Such a parameter is referred to in LLUNPIY/3r as the substate Q_KH, which is directly related to the kinetic energy.

3.2.3 Lahar outflows

The variables f[i], $\mathrm{1}\le i\le \mathrm{6}$, specify the outflows from the central cell toward the adjacent cell i; f[0] is the part remaining in the central cell. They are computed in two steps: application of the algorithm for the minimization of differences (AMD; Avolio et al., 2012; Di Gregorio and Serra, 1999) to the “heights” in the neighborhood of the central cell and calculation of the shift of f[i], $\mathrm{1}\le i\le \mathrm{6}$ (Avolio et al., 2013).

AMD application computes the outflows f[i], $\mathrm{1}\le i\le \mathrm{6}$, which minimize the height differences in the neighborhood (Eq. 7). An alteration of height values is introduced in the central cell for taking into account the outflow run-up; furthermore the viscosity is modeled by an adherence “adh” term, the lahar quantity, which cannot leave the central cell. It varies between the two extreme values adh1 and adh2, which depend on the composition of the pyroclastic debris at the maximum and minimum water content (Machado, 2015).

This adherence method was initially used for modeling lava flows by CA in order to manage the continuous variation in viscosity by cooling of lava (e.g., Avolio et al., 2006). The approximation to account for viscosity in a CA context can be intuitively explained as follows. Without bringing a system into play in which innumerable fluid layers flow over one another, at most two layers are considered. The first layer, whose maximum thickness (adh) is determined by the coefficient of viscosity, cannot move if the thickness of fluid th exceeds adh, and a second layer with thickness th-adh is considered to slide onto the first one with a friction coefficient related to viscosity:

$$\begin{array}{}\text{(4)}& {\displaystyle }h\left[\mathrm{0}\right]=Q_{\mathrm{A}}\left[\mathrm{0}\right]+Q_{\mathrm{KH}}\left[\mathrm{0}\right]+\mathrm{adh},\text{(5)}& {\displaystyle }h\left[i\right]=Q_{\mathrm{A}}\left[i\right]+Q_{\mathrm{LT}}\left[i\right],(\mathrm{1}\le i\le \mathrm{6}),\text{(6)}& {\displaystyle }q=Q_{\mathrm{LT}}\left[\mathrm{0}\right]-\mathrm{adh}=\sum _{\mathrm{0}\le i\le \mathrm{6}}f\left[i\right]\text{(7)}& {\displaystyle }\sum _{\mathit{\{}\left(i,j\right)/\mathrm{0}\le i<j\le \mathrm{6}\mathit{\}}}\left(\mathrm{|}\left(h\left[i\right]+f\left[i\right]\right)-\left(h\left[j\right]+f\left[j\right]\right)\mathrm{|}\right).\end{array}$$

Each moving quantity f[i] (where $\mathrm{1}\le i\le \mathrm{6}$ in the following equation) may be considered to be a “cylinder”, at first entirely inside the cell, having a center of mass with coordinates Q_X[0] and Q_Y[0] and with the maximum possible radius.

The shift “sh[i]” of f[i] is calculated according to the following formulae, where the movement of the mass center is specified as the mass movement on a constant slope with a constant coefficient of friction p_cf. The movement of f[i] is directed towards the center of cell i with coordinates Q_X[i] and Q_Y[i], considering the slope angle θ[i] (Fig. 4):

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \mathrm{sh}\left[i\right]=v\cdot p_{\mathrm{t}}+g\left(\sin \mathit{\theta }\left[i\right]-p_{\mathrm{cf}}\cos \mathit{\theta }\left[i\right]\right)p_{\mathrm{t}}^{\mathrm{2}}/\mathrm{2}\\ & (\mathrm{1}\le i\le \mathrm{6}),\end{array}\end{array}$$

with “g” being the acceleration of gravity and “v” being the initial velocity,

$$\begin{array}{}\text{(9)}& v=\sqrt{\left(\mathrm{2}g\cdot Q_{\mathrm{KH}}\left[\mathrm{0}\right]\right)}.\end{array}$$

There are three possible outcomes: if the shifted cylinder is completely inside (outside) the central cell, there is only an internal (external) outflow. Otherwise two cylinders form, with the mass center corresponding to the mass center of the internal outflow and that of the external outflow. The new positions (_iQ_EX[0] and _iQ_EY[0] and _iQ_IX[0] and _iQ_IY[0]) of external and internal outflow _iQ_E[0] and _iQ_I[0] account also for the variations in kinetic head _iQ_KHE[0] and _iQ_KHI[0].

Figure 4 Outflow direction from central cell to the center of an adjacent cell in three dimensions.

3.2.4 Flow composition

Execution of the elementary process “lahar outflows” involves an updating of substates Q_LT, Q_KH, Q_X, and Q_Y by the elementary process “flow composition”. It accounts for the variation in the lahar content of the cell, i.e., variation in Q_LT and corresponding variation in Q_KH, Q_X, and Q_Y that is determined by the external outflows Q_E; these represent inflows for the cells to which they are directed. Internal outflows Q_I determine just a shift inside the cell, with variations in the substates Q_KH, Q_X, and Q_Y (Machado et al., 2015a, b).

The value of the substate Q_LT at the next step is given by its previous value minus the losses determined by the outflows _iQ_E[0] from the cell (normalized to a thickness) plus the contributions determined by the inflows _7−iQ_E[i] from the neighbors:

$$Q_{\mathrm{LT}}^{\prime }\left[\mathrm{0}\right]=Q_{\mathrm{LT}}\left[\mathrm{0}\right]+\sum _{i=\mathrm{1}}^{\mathrm{6}}\left({}_{\mathrm{7}-i}{Q}_{\mathrm{E}}\left[i\right]{-}_i{Q}_{\mathrm{E}}\left[\mathrm{0}\right]\right).$$

The other substates change correspondingly, considering a weighted average

$$\begin{array}{rl}& Q_{\mathrm{KH}}^{\prime }\left[\mathrm{0}\right]=\\ & {\displaystyle \frac{Q_{\mathrm{KH}}\left[\mathrm{0}\right]\cdot Q_{\mathrm{LT}}\left[\mathrm{0}\right]+\sum _{i=\mathrm{1}}^{\mathrm{6}}\left({}_{\mathrm{7}-i}{Q}_{\mathrm{KHE}}\left[i\right]{\cdot }_{\mathrm{7}-i}{Q}_{\mathrm{E}}\left[i\right]{-}_i{Q}_{\mathrm{KHE}}\left[\mathrm{0}\right]{\cdot }_i{Q}_{\mathrm{E}}\left[\mathrm{0}\right]\right)}{Q_{\mathrm{LT}}\left[\mathrm{0}\right]+\sum _{i=\mathrm{1}}^{\mathrm{6}}\left({}_{\mathrm{7}-i}{Q}_{\mathrm{E}}\left[i\right]{-}_i{Q}_{\mathrm{E}}\left[\mathrm{0}\right]\right)}}.\end{array}$$

The shifts, both of the external flows and internal flows, have to be considered for the new values of coordinate substates:

$$\begin{array}{rl}& Q_X^{\prime }\left[\mathrm{0}\right]=\\ & {\displaystyle \frac{Q_X\left[\mathrm{0}\right]\cdot \left(Q_{\mathrm{LT}}\left[\mathrm{0}\right]-\sum _{i=\mathrm{1}}^{\mathrm{6}}\left({}_i{Q}_{\mathrm{I}}\left[\mathrm{0}\right]\right)\right)+\sum _{i=\mathrm{1}}^{\mathrm{6}}\left({}_{\mathrm{7}-i}{Q}_{\mathrm{E}X}\left[i\right]{\cdot }_{\mathrm{7}-i}{Q}_{\mathrm{E}}\left[i\right]{+}_i{Q}_{\mathrm{I}X}\left[\mathrm{0}\right]{\cdot }_i{Q}_{\mathrm{I}}\left[\mathrm{0}\right]{-}_i{Q}_{\mathrm{E}X}\left[\mathrm{0}\right]{\cdot }_i{Q}_{\mathrm{E}}\left[\mathrm{0}\right]\right)}{Q_{\mathrm{LT}}\left[\mathrm{0}\right]+\sum _{i=\mathrm{1}}^{\mathrm{6}}\left({}_{\mathrm{7}-i}{Q}_{\mathrm{E}}\left[i\right]{-}_i{Q}_{\mathrm{E}}\left[\mathrm{0}\right]\right)}},\\ & Q_Y^{\prime }\left[\mathrm{0}\right]=\\ & {\displaystyle \frac{Q_Y\left[\mathrm{0}\right]\cdot \left(Q_{\mathrm{LT}}\left[\mathrm{0}\right]-\sum _{i=\mathrm{1}}^{\mathrm{6}}\left({}_i{Q}_{\mathrm{I}}\left[\mathrm{0}\right]\right)\right)+\sum _{i=\mathrm{1}}^{\mathrm{6}}\left({}_{\mathrm{7}-i}{Q}_{\mathrm{E}Y}\left[i\right]{\cdot }_{\mathrm{7}-i}{Q}_{\mathrm{E}}\left[i\right]{+}_i{Q}_{\mathrm{I}Y}\left[\mathrm{0}\right]{\cdot }_i{Q}_{\mathrm{I}}\left[\mathrm{0}\right]{-}_i{Q}_{\mathrm{E}Y}\left[\mathrm{0}\right]{\cdot }_i{Q}_{\mathrm{E}}\left[\mathrm{0}\right]\right)}{Q_{\mathrm{LT}}\left[\mathrm{0}\right]+\sum _{i=\mathrm{1}}^{\mathrm{6}}\left({}_{\mathrm{7}-i}{Q}_{\mathrm{E}}\left[i\right]{-}_i{Q}_{\mathrm{E}}\left[\mathrm{0}\right]\right)}}.\end{array}$$

3.3 LLUNPIY calibration and validation

We selected the 2005 and 2008 lahars of Vascún valley for LLUNPIY/3r version calibration and validation, respectively. Available data, although incomplete, of the flood phase (Machado et al., 2015b) seemed promising in obtaining reliable simulations. In fact, data of different sources were carefully compared and analyzed (Williams et al., 2008; IGEPN, 2008a, b) in order to reconstruct the two events as accurately as possible (Machado et al., 2015a, b).

The use of simulation tools (from the cellular-automata model LLUNPIY) requires detailed field data: the digital elevation model (DEM) and depth of erodible pyroclastic stratum. This implies accurate geological investigations, including subsoil tomography. Geophysical surveys allow individuating points where dams by backfill, which are easy to collapse, can enable the formation of ponds, whose breakdown can trigger a lahar (Machado, 2015; Chidichimo et al., 2016).

The simulation of the 2005 event is based on a DEM with a 1 m cell size (supplied to us by Gustavo Cordoba), while the 2008 lahar was performed with a DEM with a 5 m cell size (supplied by IGEPN). In both cases a uniform thickness of 5 m was imposed for detrital cover because detailed surveys were not available. This introduces a series of approximations that negatively influence the results of simulations. Such approximations can be reduced by an opportune survey of field data, e.g., by soil tomography, multichannel analysis of surface waves (MASW), coring, etc.

The same set of LLUNPIY/3r parameters was used in the two cases except for the parameter of progressive erosion (p_pe) because of different percentages of water in the soil. The 2005 event was triggered in a higher and very steep zone of Vascún River, when the water concentration in the soil reached critical values due to rainfall. The 2008 event was dissimilar because the breaking of a temporary pond suddenly released a larger water quantity (in comparison with 2005 case) with strong turbulence, whose effects correspond to a higher value of the parameter of progressive erosion (Machado et al., 2015b).

The results of the simulations of the 2008 event (Machado et al., 2015b) are extensively reported in this study, since this event, as previously mentioned, was caused by the same type of phenomenon whose development we want to forecast. The reliability of the simulation results, in comparison with the real event, led us to confide in the goodness of the method. The new simulations were performed with the same data precision and the same values of parameters.

Simulations of the 2005 event were limited by the partial data field and digital terrain model (DTM) information: a stretch of about 2.3 km, from an elevation of 2150 m a.s.l. (about 850 m upstream of El Salado baths) to an elevation of 1900 m a.s.l. (in correspondence of Pastaza River), was considered. The area where the simulation starts does not concern the detachment phase that occurs 8 km upriver. A kind of detachment where an initial velocity of 7 m s⁻¹ was imposed to lahar was considered in order to express the first arrival of lahar flows. An equivalent fluid approach was adopted because precise data about water flows are not available. Therefore, bulking must account not only for the erodible layer but also for water inclusion. The total mass is inclusive of the water mass. This generates a discrepancy between the lahar volume, measured on the deposit, and the “fluid” lahar volume including water to be lost in the last part of the event.

Figure 5 (a) Maximum thickness, (b) maximum velocity, and (c) erosion depth, in 2005 simulated event.

Table 3 Comparison among field data and TITAN2D and LLUNPIY simulation data.

Figure 5 shows the simulation developed with LLUNPIY in the considered sector. In particular, the maximum debris thickness values, which were reached by the lahar in simulation, are reported in Fig. 5a. Maximum velocities, reached by simulated flows (Fig. 5b), are high in steeper areas (the expected result) and gradually decrease at the downstream outlet. A velocity increase occurs south of Baños, probably because of the higher gradient of the river bed. Erosion has a trend similar to that of velocity (Fig. 5c). Table 3 synthesizes values of Fig. 5 and compares such data with IGEPN field data, reported in IGEPN 2005, and with simulation performed by TITAN2D (Williams et al., 2008). Such field data are obviously a part of the complete development of catastrophic phenomenon but extremely precious for the comparison with our simulations. Observation data are not sufficient for a precise comparison with the simulation paths. Furthermore, the lahar starts with null velocity in the simulation of (Williams et al., 2008), while LLUNPIY simulations start with 7 m s⁻¹ velocities. The difference for total eroded mass rises from the lost water volume that was not possible to be considered in measurements.

Figure 6 The 2008 simulated event. (a) Maximum thickness, (b) maximum velocity, and (c) erosion depth.

The simulation of the 2008 lahar is shown in Fig. 6: the flow speed reached values up to 20 m s⁻¹ in many areas of the valley, and the eroded material resulted in a volume of about 970 000 m³. The maximum height obtained in the simulated flow (Fig. 6a), in some sectors where the valley is particularly narrow, is 22 m, while the estimated average value by IGEPN (2008a, b) is 4 m.

Table 4 compares some data of the simulation by LLUNPIY with corresponding field data of IGEPN (2008a, b). It is possible to note that the data derived from the simulations are not very different from the known measured ones. The flow velocity of 15 m s⁻¹ represents an estimated value and not a measured one. These results demonstrate that LLUNPIY/3r is a reliable model given that the simulations are based on incomplete and sometime very approximate data concerning the pre-event and post-event. Furthermore the inevitable errors, in the records related to this event, have to be considered. Therefore an extension of LLUNPIY/3r is promising in order to introduce secondary features of the phenomenon to be tested. Simulations reproduce satisfactorily the overall dynamics of the events: there is good matching between the real and simulated lahar path and the velocity and height of detrital flow. Note that different approaches always obtain excellent results about the path because the lahar is canalized by steep faces.

Table 4 Comparison between field and LLUNPIY simulation data.

4 Lahar triggering and effects

4.1 Building rudimental dams that are easy to collapse

Ponds form along watercourses in volcanic areas, when landslides of volcanic deposits, which are originated by pyroclastic flows and lahars, act as a dam by obstructing the stream bed. The most frequent cause of a breakout of such natural ponds is the overflow of water across the newly formed dam during violent rainfall and subsequent erosion and rapid cutting down into the loose rock debris. The classification of the Glossary of Meteorology of the American Meteorological Society for rainfall intensity (Glickman, 2000) is here adopted according to the rate of precipitation R_p (measured in mm h⁻¹):

• light rain – R_p<2.5 mm h⁻¹,

• moderate rain – 2.5 mm h${}^{-\mathrm{1}}\le R_{\mathrm{p}}<\mathrm{10}$ mm h⁻¹,

• heavy rain – 10 mm h${}^{-\mathrm{1}}\le R_{\mathrm{p}}\le \mathrm{50}$ mm h⁻¹,

• violent rain – R_p>50 mm h⁻¹.

Dam collapse occurs when instable conditions arise on the downstream slope. By eroding the obstruction and flowing downstream along the river bed, the initial surge of water will incorporate a dangerous volume of sediments. This can easily produce lahars, with possible devastating effects for settlements in their path (Leung et al., 2003).

Temporary dams with a similar (but controlled) behavior can be designed and built at low cost by local backfill in order to allow the outflow of streams produced by regular rainfall events. This result is achieved by properly dimensioning the embankment according to a stability analysis. This is made by comparing the forces tending to cause movement of the mass of pyroclastic material (force of water) with those tending to resist the movement (soil strength; Lambe and Whitman, 1979). The aforementioned approach is traditionally adopted to prevent dams failure, but it will be used, in this case, to ensure their collapse at a specific water level. During rainfall events, in fact, the barred canal section fills up rather quickly, so the hypothesis behind the simulations is that the dam reaches the instability conditions for the achievement of a given hydraulic head rather than for other processes (e.g., erosion), since the first destabilizing condition is reached faster than the others which require longer times to be effective.

The FEM was applied to perform the downstream slope stability analysis. The shear strength reduction (SSR) approach, which is one of the most popular techniques to perform FEM slope analysis, was adopted (Griffiths and Lane, 1999). The SSR is simple in theory: it systematically reduces the shear strength envelope of material by a factor of safety (FS) and computes FEM models of the slope until deformations are unacceptably large or solutions do not converge. In classical soil mechanics, the factor of safety is the ratio of the shear strength at the plane of potential failure τ_f and the shear stress acting in the same plane τ, namely

$$\begin{array}{}\text{(10)}& \mathrm{FS}={\mathit{\tau }}_{\mathrm{f}}/\mathit{\tau }.\end{array}$$

For the Mohr–Coulomb criterion, the previous equation reads

$$\begin{array}{}\text{(11)}& \mathrm{FS}=\left(c+(s_{\mathrm{n}}-u)\cdot \tan \mathit{\varphi }\right)/\mathit{\tau },\end{array}$$

where c is the cohesion and ϕ is the friction angle of the material; s_n denotes the total normal stress and u the pore pressure. For the Mohr–Coulomb model, a “reduced” set of material parameters c^∗ and ϕ^∗ is computed:

$$\begin{array}{}\text{(12)}& {\displaystyle }{c}^{\ast }=c/{\mathrm{FS}}_{\mathrm{n}},\text{(13)}& {\displaystyle }\tan {\mathit{\varphi }}^{\ast }=\tan \mathit{\varphi }/{\mathrm{FS}}_{\mathrm{n}}.\end{array}$$

The problem is then solved using this set of reduced material parameters while keeping all other parameters unchanged. If convergence is obtained, the FS_n value is increased and the problem is solved again. The lowest FS_n producing non-convergence is reported as the “factor of safety” of the problem. If the resulting FS is greater than 1, for a given water level stressing the upstream slope of the dam, the structure is stable. If the iterative procedure gives back a unitary FS value, the limit equilibrium conditions have been reached. This means that the coupling of both the dam geometrical configuration and the water level situation are going to produce the collapse of the structure. This last condition is the one that must be reached for the study purposes. The pore pressure distribution inside the dam body is fundamental input for the strength reduction analysis. Such distribution is obtained as a function of the hydraulic head stressing the upstream slope of the embankment. The filtration process, implemented in the finite element model, is based on the solution of the Laplace equation (Straface et al., 2010, 2011; Molinari et al., 2014; Chidichimo et al., 2015, 2018):

$$\begin{array}{}\text{(14)}& {\mathrm{\nabla }}^{\mathrm{2}}h=\mathrm{0},\end{array}$$

where h(x,y) represents the hydraulic head distribution within the dam body, which is a function of the hydraulic conductivity of the dam material. Such dependence is defined by Darcy's law:

$$\begin{array}{}\text{(15)}& q=-K\mathrm{\nabla }h=-K{\displaystyle \frac{\mathrm{d}}{\mathrm{d}l}}\left({\displaystyle \frac{u}{\mathit{\rho }g}}+z\right),\end{array}$$

where q is Darcy's velocity, K is the hydraulic conductivity, ρ is the water density, u is the pore pressure, and z is the elevation above sea level. Several simulations were performed by taking the parameters for the numerical model from the literature. Studies performed on the geotechnical properties of the volcaniclastic formation that is found in the Andes of Ecuador and Colombia, known as Cangahua, reported that the dry unit weight of the material was found to range around 14 kN m⁻³ (Bommer et al., 2002). The strength parameters of volcanic sediments produced by recent eruptions were investigated by Orense et al. (2006), who found a value for the friction angle (ϕ) of these materials of about 40^∘, while the cohesion (c) is close to zero. The hydraulic conductivity (K) of pyroclastic beds was discovered to range between 10⁻⁴ and 10⁻⁵ m s⁻¹ (Burgisser, 2012). A middle-range value was adopted to implement the numerical models. The relatively high permeability of the pyroclastic material ensures the water outflow during the regular rainfall regimen. In case of extreme rainfall events (violent rains, typhoons, etc.), the volcanic material is no longer able to drain the water inflow, producing the water level rise which will undermine the structural stability. The dams were thought to reach a height of 3 m and to hold up to a maximum water level of 2.6 m. Assuming the aforementioned features, the first step started from the design of a stable dam configuration (Fig. 7a). This outcome was obtained by trying different dimensioning solutions in order to avoid the early collapse of the structure due, for example, to its own weight. Figure 7b shows the water velocity field in the dam cross section, while Fig. 7c shows the breaking surface causing the failure of the downstream slope, which is generated by a factor of safety of 1.28.

Figure 7 (a) Cross-sectional sketch of a stable dam with the main elements dimension, (b) water velocity field moving in the dam body and associated phreatic surface, and (c) breaking surface generated by a factor of safety of 1.28.

Once a stable dam was obtained for the chosen working conditions, the second step was to repeat the strength reduction analysis several times by slowly increasing the inclination of the downstream slope until a unitary FS was reached. The inclination was increased using the corner between the downstream slope and the dam crest as a pivot point. This resulted in a gradual increase in the θ angle and a consequent reduction of the dam base. Figure 8 shows the analysis final result, with the sizing specifications for an easy-to-collapse dam built in volcaniclastic material. To ensure a greater control over the natural dam collapse timing, a discharge channel can be arranged at the dam base. The degree of openness of this channel can be adjusted according to the flow rate values observed during the extreme rainfall events recorded over time in the area in order to delay the achievement of the triggering hydraulic head. This ploy may be necessary to avoid the undesired simultaneous collapse of different dams; the hazard could increase when different lahars are triggered at short time intervals and reach the confluence points almost at the same time.

Figure 8 (a) Cross-sectional sketch of an instable dam with the main elements dimension, (b) water velocity field moving in the dam body and associated phreatic surface, and (c) breaking surface generated by a unitary factor of safety.

4.2 Preliminary hypotheses and results of simulations

LLUNPIY was calibrated and validated for secondary lahars by simulating the two events of 12 February 2005 and 22 August 2008, which occurred in the Vascún valley of Tungurahua volcano in Ecuador (Machado et al., 2014, 2015b). In particular, the 2008 event is very important in order to confirm the value of the model parameters, tuned in the simulation even where the cause of the lahar was the breakdown of a temporary pond, generated by a small landslide. Those successful simulations permitted having confidence in the scenarios which could be realized by new simulations. Of course, a very accurate update of geological data (DEM or DTM, detrital cover depth, etc.) and sufficient geophysical surveys are indispensable for applications aimed at lahar risk mitigation.

An initial study on the potentiality of applying mitigation measures in the Vascún valley was performed by triggering lahars of a planned size (the lahar level is here considered to be the relevant datum) through the controlled collapse of rudimentary ponds.

A preliminary analysis of the principal canyons of the Vascún valley was carried out in order to individuate favorable points for positioning embankments as dams. Three points were chosen for building temporary dams: one located in Vascún River (point 1 in Figs. 9, 10, 11, and 12), the second located in one of its tributary to the right (point 2 in Figs. 9, 10, 11, and 12), and the third located in a tributary to the left (point 3 in Figs. 9, 10, 11, and 12). Vascún River in turn is a tributary of the much broader Pastaza River, where lahars of Vascún valley disperse. Simulations concern the lahars generated in the points 1, 2, and 3 (Figs. 9, 10, 11, and 12), in short, lahar 1, 2, and 3. The same initial volume of the 2008 event was selected for all the simulations except the last one. Three initial points permit the analysis of an almost-exhaustive set of possible conditions; we performed a sequence of simulations by LLUNPIY/3r, of course, with the same parameters values of the successful simulation of the 29 August 2008 event. We present here some selected simulations that appear to be interesting for many considerations, which may be deduced by their analysis.

Figure 9 Simultaneous triggering of lahars from points 1 and 2 (a) and from points 2 and 3 (b). The maximum thickness of the lahar during the conjectured events is reported in meters according the legend.

Figure 10 Single triggering events: lahar generated in point 1 (a), in point 2 (b), and in point 3 (c). The legend specifies the maximum thickness of the lahar reached during the conjectured events for each point (cell).

Figure 11 Simulations of deferred triggering of lahars generated in point 3 (a) and 1 (b). The legend specifies the maximum thickness of the lahar reached during the conjectured events for each point (cell).

Figure 12 Simulations of simultaneous triggering of lahars generated in point 1, 2, and 3: in panel (a) with total detachment volume of 13 625 m³ and in panel (b) with total detachment volume of 27 250 m³. The legend specifies the maximum thickness of the lahar reached during the conjectured events for each point (cell).

Initially, we considered two scenarios in order to investigate the effects of simultaneous and time-varying confluence of two lahars, the results of such simulations caused us to consider a larger set of cases. The former is generated by the simultaneous triggering from the points 1 with a detachment volume of 4875 m³ and 2 with a detachment volume of 4500 m³ (Fig. 9a), and the latter is generated by the simultaneous triggering from the points 2 with the same previous detachment volume and 3 with a detachment volume of 4250 m³ (Fig. 9b). The confluence of the lahars in the Vascún River (for the following the confluence point) is almost simultaneous in the latter scenario because of the similar distance between the triggering points and the confluence point. The situation is diverse for the former scenario because the distances of the triggering points from the confluence point are very different. Intuitively the former case would be less dangerous because the flood peaks of the two lahars cannot coincide at the confluence point, but the maximum thickness of the lahar in the former scenario is 27 m, while the smaller value of 21 m is detected in the latter scenario. The analysis of the two simulations showed that a much larger mass was eroded in the first part of the path from the point 1 in the former case. This unexpected result permits planning an opportune strategy according to the degree of control for triggering lahars 1, 2, and 3 from the three respective points. If triggering can be well controlled with moderate or heavy rainfall, then the best choice is to trigger lahar 3 (smaller erosion) before lahar 1 so that lahar 3 anticipates part of the erosion process in the common path of both the lahars (1 and 3) as far as the confluence point and reduces, consequently, the thickness of lahar 1. When the peak of lahar 3 goes beyond the confluence point, then lahar 2 can be triggered before lahar 1, which has to be generated as late as possible. Anyway, the last lahar to be triggered has to be surely lahar 1, but a further investigation is necessary to better understand the priority between lahar 2 and 3; the study of single lahars generated in points 1, 2, and 3 could solve the question, as it may be deduced by the following simulations.

Lahar 1 causes the maximum erosion, with a maximum thickness of 26 m because it follows the path of Vascún River, which is the main river in the valley (with a larger volume of pyroclastic cover to be eroded); lahar 3 shares a relevant part of the previous path and reaches a maximum thickness of 19 m, while lahar 2, whose path is shorter before its late confluence into Vascún River, involves the smallest erosion (maximum thickness of 16 m). Such results solve the doubt that we put forward with the first simulations.

An “cleaning” operation of pyroclastic cover could be projected by initially triggering lahar 2, with a first mobilization of the detrital cover also for the area related to the last part of Vascún River from the confluence point of lahar 2 (maximum thickness 16 m; Fig. 10b). When lahar 2 dissolves into Pastaza River, lahar 3 could be triggered with a first erosion of the detrital cover between the confluence points of lahars 2 and 3 into Vascún River so that the maximum thickness does not exceed 16 m (Fig. 10c). The most dangerous lahar 1 of Vascún River could then start at the exhaustion of lahar 3, minimizing the hazard; the maximum thickness before the confluence of lahar 3 into Vascún River does not exceed 22 m (Figs. 9a and 10a).

We tested successfully the outcomes of this strategy by simulating the triggering of the three lahars in successive times, each one immediately after the exhaustion of the previous one; the first phase concerns lahar 2 (Fig. 10b). The maximum thickness does not exceed 14 m in the last part of the path, from north of El Salado baths to south of Baños.

The erosion depth of pyroclastic cover, after lahar 2 exhaustion, prevents the maximum thickness of the successive lahar 3 to exceed 10–12 m after the confluence point with lahar 2 (Fig. 11a) because of the reduced pyroclastic cover, while 19 m are reached by triggering only lahar 3 (Fig. 10c).

Finally, the most dangerous lahar 1 does not exceed 4–8 m in the inhabited zones (Fig. 11b), while it reaches 26 m of maximum thickness in the last part of the path (Fig. 10a), when the other lahars 2 and 3 are not generated.

This last result points out the importance of cleaning the pyroclastic cover according to an opportune strategy, which can be deduced by the outcomes of simulations that explore all the possible significant cases.

The final simulations concern two cases of simultaneous triggering of all the lahars, with the same detachment volumes from points 1, 2, and 3 of the previous simulations (Fig. 12a) and the double detachment volumes from the same points (Fig. 12b) in order to understand how triggering larger volume could increase the lahar dangerousness. The results show that the maximum thickness of the lahar in the former case (Fig. 12a) is 28 m, while the maximum thickness of the lahar in the latter case (Fig. 12b) is 29 m, just a meter more. A double initial volume does not involve much larger erosion in this context, and the joint effect of a larger volume and erosion does not increase the hazard in a dramatic way.

The overall results confirm the goodness of the strategy of triggering lahars at different times according to an accurate analysis of simulations after gaining precise knowledge of the geological features of the area of application. We remember that these simulations were obtained without sufficient data about the pyroclastic cover (it is obviously overestimated) that would have led to more accurate results anyway; this issue does not compromise the reliability and validity of the proposed methodology.

5 Conclusions and comments

Our case study starts from examination of some significant natural events in the Vascún valley, an area that is heavily exposed to lahar risk. Ponds form along watercourses (rios) in the canyons of this volcanic zone, when landslides of volcanic deposits impede the normal flow. The most frequent cause of breakout of such natural ponds is the overflow of water across the newly formed dam and subsequent erosion and rapid cutting down into the loose rock debris. By eroding the blockage and flowing out watercourse downstream, the initial surge of water will incorporate part of volcanic sediments and will generate lahars. The hazard related to these lahars depends both on the features of the temporary pond and the volcanic cover along the lahar path; a larger frequency of lahars produces smaller (shorter accumulation periods) and therefore less dangerous events.

We explored the possibility of inducing lahars artificially and performed many simulations for analyzing possible different scenarios for extremely complex situations; positive results obtained from these studies permit us to settle a methodology and encourage us to continue this investigation. We propose a controlled generation of small lahars, for risk mitigation, by the collapse of temporary ponds at different times in order to avoid the superposition of lahars having the same final path. Such a proposal does not adhere to standards and is based on observations and study of favorable situations together with the usage of a robust and well-validated model of simulation in order to choose the best procedures of intervention.

The computational paradigm of cellular automata has allowed the development of a reliable model for the simulation of lahars complex dynamics. Reliable simulation tools give the opportunity to test various hypotheses and to create related scenarios to be analyzed.

LLUNPIY/3r, the model that was used for lahar simulations in the Vascún valley, is a reduced version of LLUNPIY. It does not account for the preliminary phase modeled by the fully extended LLUNPIY (Machado, 2015), when the mixing of the rain water with the unconsolidated pyroclastic stratum originates the lahar, but it directly considers a “detachment area”, the initial area where the lahar can be considered to start for simulations of both real and conjectured events. Since some data of the real event are missing at this stage, the simulation starting point corresponds to the first area crossed by the lahar, whose data have a good level of reliability.

Simulation results of lahars triggered by the collapsed dam are oversized: the field data, relative to the depth of the unconsolidated pyroclastic layer along the path of lahars, are known very approximately; a constant value of 5 m was adopted, certainly not lower than the real one in any part of the lahar path but possibly exaggerated in some other parts. Of course, it was preferred to consider, with poor data, an overestimated lahar hazard rather than an underestimated one, but in the future, better precise data of the unconsolidated pyroclastic layer can be obtained thanks to geophysical surveys. Besides this, even if the simulations of lahars triggered by the collapsed dam produce overvalued hazard scenarios, the comparison among all the cases, where the depth of the unconsolidated pyroclastic layer is overestimated in the same way, shows that the application of this methodology with accurate field data is worthy of being taken into consideration. We remember that, in Ecuador, the two most accredited models for lahar simulation, LAHARZ (Muñoz-Salinas et al., 2009; Schilling, 1998) and TITAN2D (Sheridan et al., 2005; Williams et al., 2008), omit the erosion process. They impose the total amount of the eroded pyroclastic layer at the first simulation step, while LLUNPIY/3r starts from an initial amount, adding the new eroded quantity according to a step-by-step computation of the erosive detachment.

LLUNPIY/3r is limited for application to Vascún valley (or similar cases) because all the lahars end in Pastaza River without significant variation in viscosity, so the possible last phase affecting lahars in areas with small slopes, i.e., the water loss and the resulting solidification, fails to be considered; LLUNPIY (Machado, 2015) models such a situation, but a reliable validation of the model requires simulation of opportune real cases with detailed field data.

The possibility of simulating different scenarios with reliable field data permits forecasting the thickness of lahars, their velocity, and times of their peaks to operate for the best choice as a potential hazard with more efficient and reliable alert procedures. Applications of LLUNPIY/3r need a thorough geological study of the area of interest, especially regarding morphology (DEM and DTM), pyroclastic soil cover, and the composition of the erodible layer, also specified by geophysical surveys at the strategic points. Furthermore, it is also important to conduct a hydrological study of watercourses, where, most likely, the lahars are channeled.

Feasibility studies confirmed the previous hypotheses of building weak dams with a control drainage canal with significant cost containment. Unexpected (and sometime dangerous) situations were evidenced by simulation results, which permit evaluating the hazard of possible choices. Furthermore, more efficient and efficacious early warning protocols may be produced in such a context, and the social impact for partial evacuation could be mitigated. Interventions that solve provisionally the lahar hazard but involve future risks can be avoided. The complexity of the objective presupposes a multidisciplinary (or perhaps transdisciplinary) approach, which implies an even greater effort to ensure those competencies of different types (which are reflected by the various scientific extractions of the authors): geology, physics, mathematics, engineering, and computer science can cooperate to achieve common goals. Protocols for mitigating the lahar risk can be developed in such a context, involving also social and political sciences (Leung et al., 2003; Mercer et al., 2010).

A pilot project in which the setting up of an artificial pond is planned to conduct experiments for the triggering of small lahars in safe conditions would be necessary. This is a preliminary step for standard applications of this new strategy for reducing lahars risk.

A further achievement is the extension of LLUNPIY for modeling the flows in the urbanized area, as in the last versions of SCIDDICA (Lupiano et al., 2016, 2017).