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 1Tungurahua 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.1 Introduction to the LLUNPIY/3r version

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

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.

Download

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 1Physical and empirical parameters.

Download Print Version | Download XLSX

$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 2Substates.

Download Print Version | Download XLSX

σ: 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.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 4Outflow direction from central cell to the center of an adjacent cell in three dimensions.

Download

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.

Download

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

Download Print Version | Download XLSX

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 6The 2008 simulated event. (a) Maximum thickness, (b) maximum velocity, and (c) erosion depth.

Download

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 4Comparison between field and LLUNPIY simulation data.

Download Print Version | Download XLSX

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

For the Mohr–Coulomb criterion, the previous equation reads

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:

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):

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:

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.

Download

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.

Download

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 9Simultaneous 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.

Download

Figure 10Single 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).

Download

Figure 11Simulations 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).

Download

Figure 12Simulations 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).

Download

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.