# Probabilistic, scenario-based hazard assessment for pyroclastic density currents at Tungurahua volcano, Ecuador

A. Aravena, Alessandro Tadini, A. Bevilacqua, Pablo Samaniego, B. Bernard, S. Hidalgo, J. -L. Le Pennec, P. Martínez-Yáñez, J. García, Olivier Roche

<figure><img src="image_1.png" /></figure>

► To cite this version:

A. Aravena, Alessandro Tadini, A. Bevilacqua, Pablo Samaniego, B. Bernard, et al., Probabilistic, scenario-based hazard assessment for pyroclastic density currents at Tungurahua volcano, Ecuador. Bulletin of Volcanology, 2024, 86, 10.1007/s00445-024-01768-y. insu-04730597

HAL Id: insu-04730597

https://insu.hal.science/insu-04730597v1

Submitted on 10 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Probabilistic, scenario-based hazard assessment for pyroclastic density currents at Tungurahua volcano, Ecuador

A. Aravena<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> · A. Tadini<fnref n="2" />,<fnref n="3" /> · A. Bevilacqua<fnref n="3" /> · P. Samaniego<fnref n="2" />,<fnref n="4" /> · B. Bernard<fnref n="4" /> · S. Hidalgo<fnref n="4" /> · J.-L. Le Pennec<fnref n="5" /> · P. Martínez-Yáñez<fnref n="1" /> · J. García<fnref n="4" /> · O. Roche<fnref n="2" />

Received: 11 March 2024 / Accepted: 16 August 2024 / Published online: 23 September 2024
© The Author(s) 2024

## Abstract

We assess the volcanic hazard posed by pyroclastic density currents (PDCs) at Tungurahua volcano, Ecuador, using a probabilistic approach based on the analysis of calibrated numerical simulations. We address the expected variability of explosive eruptions at Tungurahua volcano by adopting a scenario-based strategy, where we consider three cases: violent Strombolian to Vulcanian eruption (VEI 2), sub-Plinian eruption (VEI 3), and sub-Plinian to Plinian eruption (VEI 4–5). PDCs are modeled using the branching energy cone model and the branching box model, considering reproducible calibration procedures based on the geological record of Tungurahua volcano. The use of different calibration procedures and reference PDC deposits allows us to define uncertainty ranges for the inundation probability of each scenario. Numerical results indicate that PDCs at Tungurahua volcano propagate preferentially toward W and NW, where a series of catchment ravines can be recognized. Two additional valleys of channelization are observed in the N and NE flanks of the volcano, which may affect the city of Baños. The mean inundation probability calculated for Baños is small (6 ± 3%) for PDCs similar to those emplaced during recent VEI 2 eruptions (July 2006, February 2008, May 2010, July 2013, February 2014, and February 2016), and on the order of 13 ± 4% for a PDC similar to that produced during the sub-Plinian phase of the August 2006 eruption (VEI 3). The highest intensity scenario (VEI 4–5), for which we present and implement a novel calibration procedure based on a few control points, produces inundation areas that nearly always include inhabited centers such as Baños, Puela, and Cotaló, among others. This calibration method is well suited for eruptive scenarios that lack detailed field information, and could be replicated for poorly known active volcanoes around the world.

**Keywords** Hazard assessment · Pyroclastic density currents · Tungurahua volcano

Editorial responsibility: Editorial responsibility: S. Massaro

<figure><img src="image_3.png" /></figure>

A. Bevilacqua
andrea.bevilacqua@ingv.it

## Introduction

Tungurahua volcano (1.47° S; 78.44° W), located ~ 8 km south of the city of Baños, is one of the most hazardous volcanoes in Ecuador and South America (Fig. 1). The eruption record of Tungurahua includes both effusive and explosive activity, with a series of events documented in historical times, such as those of 1640–1645, 1773–1782, 1885–1888, and 1916–1925 (Hall et al. 1999; Le Pennec et al. 2016) and, more recently, an eruptive episode that lasted since 1999 until 2016 (Bernard et al. 2016; Samaniego et al. 2011; Vlastélic et al. 2023). In addition, during the Late Holocene, Tungurahua experienced Plinian eruptions and sector collapses (Le Pennec et al. 2008, 2013; Bablon et al. 2018), which demonstrate the potential of Tungurahua volcano to feed long-runout distance pyroclastic density currents (PDCs), posing a permanent

<footnote n="1">¹ Facultad de Ciencias Básicas, Universidad Católica del Maule, Talca, Chile</footnote>

<footnote n="2">² Laboratoire Magmas Et Volcans, Université Clermont Auvergne, CNRS, Clermont-Ferrand, IRD, OPGC, France</footnote>

<footnote n="3">³ Istituto Nazionale Di Geofisica E Vulcanologia, Sezione Di Pisa, Pisa, Italy</footnote>

<footnote n="4">⁴ Instituto Geofísico, Escuela Politécnica Nacional, Quito, Ecuador</footnote>

<footnote n="5">⁵ Geo-Ocean, Université de Bretagne Occidentale, CNRS, Ifremer, UMR6538, IRD, IUEM, Plouzané, France</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Shaded relief map of Tungurahua volcano and its surroundings, including the location of the main inhabited zones (black labels), rivers (blue labels), and some of the main ravines (green labels). The ~3 ky BP collapse scar is indicated by a dotted line (modified from Bablon et al. 2018), while the purple line represents the limit between Tungurahua and Chimborazo Provinces. Coordinates are expressed in DD notation</figcaption></figure>

threat to the surrounding communities such as the city of Baños and other villages in the Tungurahua and Chimborazo Provinces (Fig. 1). In total, over 25,000 people live in zones that could be affected by lahars, PDCs, and other volcanic products of Tungurahua volcano (Hall et al. 1999; Samaniego et al. 2008). However, in spite of (1) the numerous contributions on Tungurahua activity that have been recently published (Anderson et al. 2018; Battaglia et al. 2019; Bernard et al. 2014; Douillet et al. 2013; Eychenne et al. 2012; Fee et al. 2010; Gaunt et al. 2020; Hall et al. 2013; 2015; Kelfoun et al. 2009; Mothes et al. 2015; Palacios et al. 2023; Parra et al. 2020; Samaniego et al. 2011); (2) the continuous effort of the Instituto Geofísico of the Escuela Politécnica Nacional (IG-EPN) in monitoring and understanding the internal dynamics of this volcano; and (3) major improvements in computational capacity, numerical models, and uncertainty quantification (Aravena et al. 2020; 2023; de' Michieli Vituri et al. 2019; 2023; Esposti Ongaro et al. 2016; Flynn and Ramsey 2020; Kelfoun 2017; Kelfoun et al. 2009; Neri et al. 2015a; Neglia et al. 2021; Sandri et al. 2016; Sobradelo and Martí 2010; Tadini et al. 2020; 2022; Tierz

et al. 2016a; 2016b), the Tungurahua volcano hazard map has not been updated since 2008 (Samaniego et al. 2008).

In this work, we present a probabilistic, scenario-based hazard assessment for PDCs produced during explosive eruptions at Tungurahua volcano. The definition of the eruptive scenarios at Tungurahua is based on the analysis of both the eruptive events preserved in the geological record and those reported in historical times (including the detailed follow-up of the 1999–2016 one), and is intended to reflect the natural variability in the activity of this volcanic system. The expected spatial distribution of volcanic products for each eruptive scenario is quantified through numerical modeling. In particular, we adopt the branching energy cone model and the branching box model (Aravena et al. 2020), using the computer programs ECMapProb 2.0 and BoxMapProb 2.0, respectively, and a set of calibration strategies described in Aravena et al. (2022) in order to sample the models’ inputs using a probabilistic approach. Thereby, this strategy allows us to define uncertainty ranges for the PDC inundation probability associated with each eruptive scenario at any site upon and around the volcano, and thus quantify the limits and strengths of the numerical estimates inferred from our hazard assessment. We do not address the evaluation of the occurrence probability of each scenario, and thus in this paper we do not include weighted combinations of the eruptive scenarios in order to generate long-term PDC hazard maps (e.g., Neri et al. 2015a; Sandri et al. 2016; Rutarindwa et al. 2019; Massaro et al. 2023) or estimates of the absolute PDC inundation probability around Tungurahua volcano for specific time windows (e.g., Bevilacqua et al. 2017; Sandri et al. 2018). Thereby, in this study we do not integrate the results in a fully probabilistic hazard assessment, which would require the construction of a frequency-magnitude model for the eruptions (Orsi et al. 2009; Marzocchi and Bebbington 2012). Our results are instead conditioned on the occurrence of the three considered eruption scenarios.

## Geological framework

### Overview

Tungurahua stratovolcano (5023 m a. s. l.; Fig. 1) is an andesitic-dacitic edifice located in the Eastern Cordillera of Ecuador, about 140 km south of Quito, constructed upon a metamorphic basement of Paleozoic and Cretaceous ages (Aspden and Litherland 1992; Litherland et al. 1993). Tungurahua is part of the Andean Northern Volcanic Zone (NVZ), a region that includes volcanoes in Ecuador and Colombia formed as a consequence of the subduction of the Nazca Plate beneath the South American Plate (Bryant et al. 2006; Nocquet et al. 2014). This volcanic edifice, which is one of the most active stratovolcanoes in the Ecuadorian

<figure><img src="image_5.png" /></figure>
Andes, presents particularly steep sides (from ~20 to 25° in the lower part up to ~40° slope in proximal domains; Bablon et al. 2018) and a complex summit morphology characterized by a series of nested structures, including an upper semi-elliptic crater elongated in the NE-SW direction and an irregularly shaped lower crater elongated in the NNE-SSW direction. Tungurahua is surrounded by three main rivers: Puela, Chambo, and Pastaza (Fig. 1). The city of Baños (~13 k inhabitants) is mainly located on a terrace overlying the riverbanks of Pastaza River (Fig. 1). The eruptive history of Tungurahua volcano includes three constructive stages with similar trends in terms of geochemistry of major and trace elements (Bablon et al. 2018), separated by major sector collapse events (Bablon et al. 2018; Hall et al. 1999; Le Pennec et al. 2013):

(a) *Tungurahua I* (>293 ± 10–79 ± 3 ka): construction of an andesitic edifice that peaked about 2 km southeast to the present summit of Tungurahua (Bablon et al. 2018). A volume of 56 ± 33 km³ was estimated for the edifice constructed during this period, with a mean eruptive rate of about 0.6 ± 0.3 km³/ka (Bablon et al. 2018). This stage was seemingly followed by ~50 ky of quiescence, and then a western sector collapse occurred at ~30–35 ky BP (Bablon et al. 2018; Le Pennec et al. 2013) leaving ~10 km³ of deposits (Bustillos 2008). The remnants of Tungurahua I can be recognized in the northern, eastern, and southern flanks of the volcano, and consist of a series of andesitic and dacitic lava flows and breccia deposits (Bablon et al. 2018; Hall et al. 1999).

(b) *Tungurahua II* (29 ± 2–~3 ka): construction of a second stratocone mainly composed of andesitic lava flows (Hall et al. 1999; Le Pennec et al. 2013), with a mean eruptive rate of about 0.9 ± 0.2 km³/ka (Bablon et al. 2018). The end of this eruptive stage was marked by a sector collapse event that resulted in ~3 km³ of deposits that came from the western flank of Tungurahua volcano (Bustillos 2008) and covered an area between 23 and 80 km² (Bablon et al. 2018). Remnants of Tungurahua II are observed on the southern upper flank and on the northern and southwestern flanks of the present volcanic edifice (Bablon et al. 2018).

(c) *Tungurahua III* (<~3 ka): construction of the present volcanic edifice by the repeated generation of lava flows, by the accumulation of deposits from pyroclastic flows and debris flows that mainly propagated through the western and northern flanks, and by the fallout of moderate volumes of tephra that were transported by the wind towards west and southwest. A mean eruption rate of 2.5 ± 1.0 km³/ka was estimated by Bablon et al. (2018) for this stage, which has been dominated by andesitic products, although three dacitic pumice fall-out deposits can be recognized as well (Le Pennec et al. 2016).

Different authors have focused on the analysis of historical data (Hall et al. 1999; Le Pennec et al. 2008, 2016), including the eruptive cycles of 1640–1645, 1773–1782, 1885–1888, and 1916–1925. PDC deposits have been associated with each of these periods (Le Pennec et al. 2016). Most were channelized through radial ravines such as Vazcún and Ulba (Fig. 1), among others, with different degrees of affectation for the city of Baños. For instance, some andesitic scoria flow deposits in the Vazcún ravine were attributed to the 1640–1645 eruptive period by Le Pennec et al. (2016). The 1773–1782 eruptive period is associated with tephra dispersion to the west, PDCs that reached the city of Baños and the Chambo and Pastaza Rivers, as well as a thick andesitic lava flow in the Juive Chico area. The 1885–1888 eruption produced extensive fall-out and PDC deposits in Juive and along other ravines of the western flank, lahars and debris flows in Vazcún and Ulba ravines (Le Pennec et al. 2016), and a thick lava flow close to the Cusua village. PDC emplacement in the Vazcún ravine is also reported for the 1916–1925 eruptive period, with little impact to Baños, but widespread impact in the western flank near Puela. Based on historical activity, a recurrence rate of about one PDC-forming eruption per century can be proposed for Tungurahua III stage (Hall et al. 1999; Le Pennec et al. 2008; 2016). Le Pennec et al. (2016) also estimated that the location of the city of Baños is impacted by PDCs on average every 350–500 years.

The last sub-Plinian eruption occurred in August 2006 (Douillet et al. 2013; Eychenne et al. 2012; Samaniego et al. 2011; Bernard et al. 2016), which generated a sustained 16–18-km-high eruptive column and multiple scoria flows that traveled along a series of ravines to the N, NW, and W from the source, and ended with the emission of a lava flow. The bulk tephra volume was of the order of 42–57 × 10⁶ m³ (Eychenne et al. 2012), while the overall volume of dense pyroclastic flow deposits was ~27 × 10⁶ m³ according to Hall et al. (2013) and 18–29 × 10⁶ m³ according to Bernard et al. (2016). This event was a paroxysmal phase of an eruptive period that started in 1999 and finished in 2016, during which volcanic activity occurred sporadically, including Strombolian, violent Strombolian, and Vulcanian events (Anderson et al. 2018; Bernard 2018; Battaglia et al. 2019; Palacios et al. 2023; Parra et al. 2016). Among the products emitted during this period, it is possible to recognize PDC deposits of a series of events that occurred in July 2006, February 2008, May 2010, July 2013, February 2014, and February 2016 (Fig. 2; Hall et al. 2015; Gaunt et al. 2020; Falasconi et al. 2023), whose deposits are here used to calibrate our numerical simulations.
<figure><img src="image_6.png" /><figcaption>**Fig. 2** Shaded relief maps of Tungurahua volcano with superposed polygons showing the footprints of a series of PDC deposits produced during recent events of this volcano (see titles). These polygons have been adopted to calibrate numerical simulations of ES1. Labels indicate the main cities (black labels) and rivers (blue labels). The ~3 ky BP collapse scar is indicated by a dotted line (modified from Bablon et al. 2018). Coordinates are expressed in DD notation</figcaption></figure>

## Definition of eruptive scenarios

Based on the eruption record of Tungurahua, we considered three eruptive scenarios of interest for PDC hazard assessment:

(a) ES1: violent Strombolian to Vulcanian eruption (VEI 2).

Events able to produce thin pyroclastic fall deposits in the volcano surroundings and to feed small-scale PDCs as a consequence of fountain collapse, low eruption column collapse or remobilization of pyroclastic material, which typically stop around a break-in-slope located at ~3000 m a.s.l. upon the Tungurahua's flank. This type of activity was frequent during the last eruption period (1999–2016), from which the deposits

<figure><img src="image_7.png" /></figure>
of six small-scale PDCs have been accurately traced (Fig. 2). These deposits, characterized though their inundation zones, were used to calibrate the numerical simulations associated with this eruptive scenario.

(b) ES2: sub-Plinian eruption (VEI 3).

This type of eruption has been common during the last ~3 ka and occurred roughly once a century (Le Pennec et al. 2008; 2016; Eychenne et al. 2012), threatening the surrounding communities. A well-documented event with these characteristics is linked to the August 2006 paroxysmal phase, where the eruptive column collapse produced a series of PDCs that reached the base of the edifice through different ravines (Hall et al. 2013; Kelfoun et al. 2009; Bernard et al. 2014, 2016). The inundation area of the PDCs produced during the August 2006 eruption was used here as a reference scenario to address this type of volcanism at Tungurahua volcano (Fig. 3).

(c) ES3: sub-Plinian to Plinian eruption (VEI 4–5).

Events able to feed long-runout PDCs related to column collapse and relatively thick fallout deposits with effects at regional scale. The recurrence of this scenario is roughly of the order of one event every 1000 years (Samaniego et al. 2008). Few field data are available to well constrain the extension of flow deposits, and thus we consider a set of control points in zones where sparse outcrops of this type of activity can be recognized (Fig. 3). In particular, we take into account two control points associated with the 1640 AD eruption ($P_1$ and $P_2$; Le Pennec et al. 2005; 2008; 2016).

Note that these scenarios differ slightly from those defined by Samaniego et al. (2008). While the lowest intensity events (scenario I) considered by Samaniego et al. (2008) are not addressed here due to the fact that no PDCs are formed during this type of eruptive activity, the intensity spectrum grouped by Samaniego et al. (2008) in the intermediate category (scenario II) is divided in two groups in this work (i.e., ES1 and ES2). The highest intensity scenarios are instead strongly consistent (i.e., scenario III and ES3).

# Methods

## The models

For the construction of PDC hazard maps, we used the computer programs ECMapProb 2.0 and BoxMapProb 2.0 (Aravena et al. 2020, 2022). The first model is based on the energy cone assumption (Malin & Sheridan 1982; Sheridan & Malin 1983; Wadge & Isaacs 1988) and suits better to describe gravitational flows; the second follows instead the box model integral formulation (Bevilacqua et al. 2022;

<figure><img src="image_8.png" /></figure>

**Fig. 3** Shaded relief maps of Tungurahua volcano with a superposed polygon showing the footprint of a PDC produced during the sub-Plinian eruption of August 2006 (modified from Bernard et al. 2014) and a few control points where PDC deposits of the AD 1640 eruption of Tungurahua volcano have been recognized (see main text). These data have been adopted to calibrate numerical simulations of ES2 and ES3. Labels indicate the main cities (black labels) and rivers (blue labels). The ~3 ky BP collapse scar is indicated by a dotted line (modified from Bablon et al. 2018). Coordinates are expressed in DD notation

Esposti Ongaro et al. 2016; Huppert & Simpson 1980; Tadini et al. 2021) and allows describing inertial flows. Both models rely on a tree-branching approach to enhance the channelization features of the models (Aravena et al. 2020), and have been already applied for the construction of PDC hazard maps (e.g., Bevilacqua et al. 2021; Aravena et al. 2023). To address the lowest intensity scenario (i.e., ES1), probably associated with remobilization of pyroclastic material or low eruption column collapse that are deposited in proximal domains around the summit area, we used the program ECMapProb 2.0, which accounts for the strong topographic control inferred from the distribution of the deposits of the benchmark eruptions (Fig. 2). On the other side, because of the possible concomitance of different PDC generation mechanisms and transport regimes (Douillet et al. 2013; Hall et al. 2013; Kelfoun et al. 2009) during VEI 3 sub-Plinian eruptions (i.e., ES2), we tested independently both ECMapProb 2.0 and BoxMapProb 2.0 to provide conservative results for this scenario. Finally, because PDCs
during large-scale explosive eruptions at Tungurahua have been likely fed from large-scale column collapse, we used the model BoxMapProb 2.0 for ES3. To obtain conservative results, we assumed that all simulated PDCs arise from the collapse of pyroclastic material in all directions. All the simulations were performed using a 16-m resolution DEM, obtained by resampling elevation data from SigTierras (Ministry of Agriculture and Livestock, Ecuador).

## Calibration of input parameters

The inputs of ECMapProb 2.0 are collapse height ($H_{0,0}$) and energy cone slope (tan ($\varphi$)), while those of BoxMapProb 2.0 include collapsing volume ($V_{0,0}$, including pyroclasts and gas), initial particle concentration ($\phi$), Froude number ($Fr$), particle sedimentation velocity ($w_s$), pyroclast density ($\rho_p$), and air density ($\rho_a$) (Aravena et al. 2020; Esposti Ongaro et al. 2016). The definition of input parameters based on physical considerations is not straightforward and thus a calibration is necessary (Aravena et al. 2022), especially for inputs for which the models are particularly sensitive (in particular, $H_{0,0}$ and tan ($\varphi$) for ECMapProb 2.0; and $V_{0,0}$, $\phi$, and $w_s$ for BoxMapProb 2.0). In this work, we calibrated the inputs following and complementing the strategies described in Aravena et al. (2022), which are based on the development of a large set of calibration simulations with a fixed vent position and two variable input parameters ($\alpha = H_{0,0}$ and $\beta = \tan (\varphi)$ for ECMapProb 2.0; $\alpha = \log (V_{0,0})$ and $\beta = \phi$ for BoxMapProb 2.0), while the other inputs, if present, are considered constant. In the case of BoxMapProb 2.0, we applied the calibration procedures three times to test the effect of different values of $w_s$ as well.

To define a structured, reproducible calibration procedure, let us consider a set of $N \times N$ calibration simulations with fixed source position and variable input parameters within predefined ranges ($\alpha \in [\alpha_1, ..., \alpha_N]$ and $\alpha \in [\beta_1, ..., \beta_N]$). If we define $S_{m,n}$ as a non-negative similarity index between the reference scenario and the calibration simulation with inputs ($\alpha_m, \beta_n$), in order to produce a calibrated probability distribution of the model inputs, we can compute the sampling probability of this pair of inputs as $P((\alpha, \beta) \approx (\alpha_m, \beta_n)) := c_p \cdot S_{m,n}$, where $c_p$ is a normalizing constant. This strategy allows us reconstructing functions of probability distribution for the model inputs that are not based on a given parametrization of a predefined probability distribution, being instead directly controlled by the values of $S_{m,n}$ computed in the calibration simulations.

For the first two scenarios, following Aravena et al. (2022), $S_{m,n}$ was defined by comparing the inundation area of the calibration simulation with inputs ($\alpha_m, \beta_n$) with the inundation area of specific, documented PDCs (see Figs. 2 and 3), adopting the following comparison metrics:

(a) Jaccard index (JI).

Intersection area between the compared inundation polygons divided by their union area. In this case, $S_{m,n}^{(1)} := JI_{m,n}^2$.

(b) Hausdorff distance (HD).

Maximum distance between a border point of one of the inundation polygons and the other inundation polygon. In this case, $S_{m,n}$ is defined by:

$$S_{m,n}^{(2)} := \left( \frac{1}{HD_{m,n} + \varepsilon_{DEM}} \right)^2 \quad (1)$$

where $\varepsilon_{DEM}$ is the cell size of the DEM used in the calibration simulations, which is included to avoid any division by zero.

(c) Root mean squared distance (RMSD).

Root of the mean squared distance between a large set (in this work, 1000) of border points of each inundation polygon and the other inundation polygon. In this case, $S_{m,n}$ is given by:

$$S_{m,n}^{(3)} := \left( \frac{1}{RMSD_{m,n} + \varepsilon_{DEM}} \right)^2 \quad (2)$$

For the sub-Plinian to Plinian scenario (i.e., ES3), because the footprint of a benchmark deposit cannot be obtained with precision from field constraints, we implemented in the programs ECMapProb 2.0 and BoxMapProb 2.0 a new calibration metric based on a series of control points ($x_i, y_i$), with $i = 1, ..., N_c$ (Fig. 3). These points are intended to represent outcrops of documented PDCs whose traceability is not enough to precisely define an inundation polygon. Let us define $d_{m,n}(x_i, y_i)$ as the minimum distance between the i-th control point and a border point of the inundation polygon derived from a calibration simulation characterized by the inputs ($\alpha_m, \beta_n$). We define the root mean squared distance to the control points as:

$$CP_{m,n} := \frac{\sqrt{\sum_{i=1}^{N_c} d_{m,n}^2 (x_i, y_i)}}{N_c} \quad (3)$$

and the associated similarity index as:

$$S_{m,n}^{(4)} := \left( \frac{1}{CP_{m,n} + \varepsilon_{DEM}} \right)^2 \quad (4)$$

Note that the approach adopted to set the input parameters differs from deterministic sampling strategies (e.g., Ferrés et al. 2013) and from strategies based on Monte Carlo sampling methods (e.g., Clarke et al. 2020). In this work, instead of using arbitrary and independent probability distributions to set the model inputs, we incorporate
data from the geological record in a structured and reproducible calibration methodology. An underlying assumption of this approach is that topographic changes of the volcano have not dramatically changed the general behaviour of PDCs during propagation at Tungurahua volcano. Table 1 presents a summary of the inputs and assumptions used in the calibration simulations.

## Construction of hazard maps

By using the different functions of sampling probability obtained from the model calibration, derived from the use of different comparison metrics (RMSD, HD and JI), reference inundation polygons, and values of $w_s$ in the case of BoxMapProb 2.0 simulations, we performed different sets of calibrated simulations ($N = 500$ for each set). In particular, for each of the six reference inundation polygons of ES1, we constructed three maps based on three different comparison metrics (RMSD, HD, and JI) for the calculation of the sampling probability functions of the input parameters, giving rise to 18 probabilistic hazard maps. For ES2, we

constructed three maps with ECMapProb 2.0, each associated with a different comparison metric (RMSD, HD, and JI), and nine maps with BoxMapProb 2.0, derived from the use of three different comparison metrics and three different values of $w_s$ in the calibration simulations (see Table 1). Finally, for ES3, we constructed three maps with BoxMapProb 2.0, considering three different values of $w_s$ in the calibration simulations (Table 1; note that only one comparison metric was considered in this case; see the “Calibration of input parameters” section). In the calibrated simulations, in contrast to the calibration simulations, the vent position was varied uniformly in the summit zone (circle with a radius of 250 m) in order to incorporate a small uncertainty affecting the collapse position. Assuming greater variability in collapse position receives little support from geological evidence at Tungurahua volcano, which to our knowledge is characterized exclusively by volcanism from the central vent in all its recent eruptions. Results are described in terms of the inundation probability; i.e., for each pixel of the resulting map, we calculated the percentage of simulations that reach this position in order to define a probability value.

<table><caption>**Table 1** Input parameters adopted in the calibration simulations used to construct probabilistic maps of PDC inundation for Tungurahua volcano</caption><thead><tr><th colspan="3">ECMapProb 2.0</th></tr></thead><tbody><tr><th colspan="3">Input parameter</th></tr><tr><td>Root energy cone height ($H_{0,0}$)</td><td>ES1<sup>1,a</sup></td><td>ES2<sup>2,a</sup></td></tr><tr><td>Energy cone slope (tan($\phi$))</td><td>100–1000 m</td><td>100–2000 m</td></tr><tr><td>Number of calibration simulations($N \times N$)<sup>+</sup></td><td>0.2–1.0</td><td>0.2–1.0</td></tr><tr><td></td><td>400</td><td>400</td></tr></tbody><thead><tr><th colspan="3">BoxMapProb 2.0</th></tr></thead><tbody><tr><th colspan="3">Input parameter</th></tr><tr><td>Collapsing volume ($V_{0,0}$)</td><td>ES2<sup>2,b</sup></td><td>ES3<sup>3,c</sup></td></tr><tr><td>Initial particle concentration ($\phi_{0,0}$)</td><td>10<sup>7</sup>–10<sup>10</sup> m<sup>3</sup></td><td>10<sup>8</sup>–10<sup>11</sup> m<sup>3</sup></td></tr><tr><td>Froude number Fr</td><td>0.005–0.040</td><td>0.005–0.040</td></tr><tr><td>Sedimentation velocity ($w_s$)</td><td>1.0</td><td>1.0</td></tr><tr><td>Pyroclast density ($\rho_p$)</td><td>0.05, 0.3, and 1.2 m/s</td><td>0.05, 0.3, and 1.2 m/s</td></tr><tr><td>Ambient air density ($\rho_a$)</td><td>1500 kg/m<sup>3</sup></td><td>1500 kg/m<sup>3</sup></td></tr><tr><td>Number of calibration simulations ($N \times N$)<sup>+</sup></td><td>1.0 kg/m<sup>3</sup></td><td>1.0 kg/m<sup>3</sup></td></tr><tr><td></td><td>400</td><td>400</td></tr></tbody></table>

<footnote n="1">¹Reference inundation polygons used for model calibration: July 2006, February 2008, May 2010, July 2013, February 2014, and February 2016 PDCs (Fig. 2)</footnote>

<footnote n="2">²Reference inundation polygon used for model calibration: August 2006 PDCs (Fig. 3)</footnote>

<footnote n="3">³Model calibration based on a set of control points due to lack of detailed field information (Fig. 3)</footnote>

<footnote>aThree comparison metrics are used to perform the model calibration (JI, HD, and RMSD), giving rise to three different sampling probability distributions of the inputs for each calibration polygon considered</footnote>

<footnote>bThree sets of calibration simulations were performed with variable values of sedimentation velocity ($w_s$). In each case, three comparison metrics were used to calibrate the model (JI, HD, and RMSD), giving rise to nine different sampling probability distributions of the model inputs</footnote>

<footnote>cThree sets of calibration simulations were performed with variable values of sedimentation velocity ($w_s$). In each case, one comparison metric was used to calibrate the model (CP), giving rise to three different sampling probability distributions of the model inputs</footnote>

<footnote>dNote that this differs from the number of calibrated simulations (i.e., performed using the calibration simulations to sample the model inputs), which is 500 for each set of simulations</footnote>
# Results

## ES1: violent Strombolian to Vulcanian eruption (VEI 2)

The computed probabilistic inundation maps for PDCs associated with violent Strombolian to Vulcanian eruptions (ES1) are displayed in Figs. 4 and 5. Numerical results reveal a significant effect of the crater shape and the ~3 ky BP collapse scar on PDC propagation dynamics and at least five dominant channelization ravines towards W, NW, NNW, NNE, and NE, two of which may involve the city of Baños and other inhabited centers to the north (Figs. 4 and 5; the Vazcún and Ulba ravines). The consistency between the dominant propagation direction of the documented PDCs and numerical results, which were not performed with predefined collapse directions, indicate that the preferred propagation directions of recent, small-scale PDCs at Tungurahua are controlled by crater morphology, volcano topography and channelization dynamics near the summit, which is a frequent characteristic of small-volume PDCs, instead of possible directional collapse processes. Results calculated using different calibration procedures are remarkably similar (Supplementary Figs. S1-S4), even though the inundation areas of the reference events are highly variable, ranging from 1.0 km² (February 2008) to 8.6 km² (July 2006), as well as their runout distances, which range between 2.3 km (February 2016) and 7.1 km (July 2006). In almost all the cases, the highest inundation probabilities were computed when the RMSD calibration was adopted (Supplementary Figures S1-S3 and Tables S1-S3), while no clear correlations are observed between the main geometrical properties of the calibration polygons (their inundation areas and runout distances; Fig. 2) and the area enclosed by different isoprobability curves in the resulting hazard maps (Supplementary Fig. S5). The February 2016 inundation polygon is associated with the worst calibration performance in ES1 simulations (for instance, best-fit Jaccard index of 0.18), while the best calibration performances are observed for the July 2006 and February 2014 inundation polygons (best-fit Jaccard index of 0.39 and 0.36, respectively). Driven by these differences, we constructed a weighted hazard map that considers the 18 hazard maps with different weights as a function of the performance of each set of calibration simulations in reproducing the reference PDC deposit. This map is presented in Fig. 5 and is remarkably similar to the mean, equally weighted map.

Regarding some relevant locations around Tungurahua volcano from a volcanic hazard point of view (Fig. 6a and Supplementary Tables S1-S3), the inundation probabilities for ES1 are low and, in general, well-constrained. In

particular, the inundation probability in Baños is about 6 ± 3% (maximum value of 11.6%), while the villages with the largest inundation probabilities for the eruptive scenario ES1 are Palitahua (12 ± 5%), Juive Grande (9 ± 5%), Chontapamba (8 ± 4%), Cusua (7 ± 4%), and Bilbao (7 ± 4%).

## ES2: sub-Plinian eruption (VEI 3)

For an eruptive scenario similar to the August 2006 sub-Plinian event (i.e., ES2), the channelization effect of a series of ravines towards NW, N, and NNE is also evident, as well as the influence of the crater shape and the ~3 ky BP collapse scar, even though a significant number of simulations is able to propagate a few kilometers towards SE (Fig. 7 and Supplementary Figs. S6-S9). The adoption of different numerical models, as well as different values for the particle sedimentation velocity in BoxMapProb 2.0 simulations, produces strongly similar probabilistic hazard maps for this eruptive scenario, giving rise to well-constrained values of inundation probability at specific locations, as illustrated in Fig. 6b and Supplementary Tables S1-S3. Under this eruptive scenario, the inundation probabilities in Baños, Ulba, and Cotaló are 13 ± 4%, 7 ± 2%, and 7 ± 3%, respectively. The villages with the largest inundation probabilities for the eruptive scenario ES2 are Cusua (21 ± 7%), Juive Grande (19 ± 7%), Pondera (19 ± 5%), Choglontus (18 ± 5%), and Bilbao (17 ± 6%), while the inundation probability at Palitahua is 14 ± 9%. Note that, compared to the eruptive scenario ES1, the increase in the inundation probability is particularly relevant at Puela, a consequence of the reduced effect of the crater and the ~3 ky BP collapse scar in the eruptive scenario ES2 when compared to ES1 (Fig. 6b and Supplementary Tables S1-S3).

Numerical calibration of BoxMapProb 2.0 simulations for ES2, which are based on the August 2006 sub-Plinian event, allowed computation of a calibrated value for the volume of pyroclastic materials transported in the PDCs, which can be obtained by multiplying both the calibrated variables, i.e., collapsing volume ($V_{0,0}$), which includes pyroclasts and gas, and initial particle concentration ($\phi$; see the “Calibration of input parameters” section). The calibrated volumes of collapsing pyroclasts depend on the adopted value of $w_s$ (Supplementary Fig. S10), with mean values (in logarithmic scale) between $1.5 \times 10^6$ m³ (HD comparison metric, $w_s = 0.05$ m/s) and $5 \times 10^6$ m³ (JI comparison metric, $w_s = 1.2$ m/s), and 50% of data ranging between $8 \times 10^5$ m³ and $8 \times 10^6$ m³. This range is slightly smaller than the documented volumes of the individual PDCs recognized during the August 2006 event (i.e., ~8.5–17.3 $\times 10^6$ m³; Hall et al. 2013). On the other hand, note that Bernard et al. (2016)

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 4** Probabilistic hazard maps for the eruptive scenario ES1 (violent Strombolian to Vulcanian eruption), considering separately six different reference PDCs to calibrate numerical simulations. In each panel, we present the mean inundation probability computed using three different comparison metrics (see supplementary Figs. S1-S4), which are indicated by a set of is probability curves (see legend) and a rainbow color scale. Black labels indicate the main cities, while the positions of Ulba and Vazcún ravines are indicated by yellow labels. Coordinates are expressed in DD notation</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 5** Mean, weighted mean, and maximum probabilistic hazard maps for the eruptive scenario ES1 (violent Strombolian to Vulcanian eruption). Inundation probabilities are indicated by a set of isoprobability curves (see legend) and a rainbow color scale. For computing the mean map, we assign the same weight to the 18 hazard maps associated with ES1, while the weighted mean map is obtained by assigning weights controlled by the performance of each set of calibration simulations in reproducing the reference PDC deposit (see main text). Black labels indicate the main cities, while the positions of Ulba and Vazcún ravines are indicated by yellow labels. Coordinates are expressed in DD notation</figcaption></figure>

did not separately present the volume of single PDC units and, therefore, their volume estimates are not expected to be comparable with our calibration results.

## ES3: sub-Plinian to Plinian eruption (VEI 4–5)

Finally, the mean and maximum inundation probabilities computed for sub-Plinian to Plinian eruptions (i.e., ES3) are presented in Fig. 8, where it is possible to recognize the potential role of Chambo, Patate, and Pastaza Rivers to channelize PDCs towards SW, NW, and NE from the volcano, respectively, while the results associated with different values of $w_s$ are displayed in Supplementary Figure S11. We show that PDCs similar to those that reached the locations presented in Fig. 3 (i.e., $P_1$ and $P_2$), where outcrops of PDCs

fed from Tungurahua volcano during the 1640 AD eruption can be recognized, are also likely to spread out in Baños and other nearby towns, such as Puela, Cotaló, and Ulba (Supplementary Tables S1-S3). Note, however, that the adopted models simulate collapse processes that propagate in all directions and thus we are not modeling PDCs produced by directed blasts. The calibrated volume of pyroclastic material involved in these flows exhibits mean values in logarithmic scale ranging between $3 \times 10^7 \text{ m}^3$ ($w_s = 0.05 \text{ m/s}$) and $1.3 \times 10^8 \text{ m}^3$ ($w_s = 1.2 \text{ m/s}$; Supplementary Fig. S11). These volumes are significantly larger than those computed for ES2, which may explain the large difference observed between the hazard maps associated with these eruptive scenarios (see Supplementary Tables S1-S3 and Figs. 7 and 8 for comparison).
<figure><img src="image_12.png" /><figcaption>**Fig. 6** Inundation probability computed in a series of critical positions around Tungurahua volcano, considering different scenarios and calibration methods. For ES1 (panel a), data are presented in box plots, while each symbol represent a hazard map in panel b (ES2 and</figcaption></figure>

ES3, see legend). Yellow symbols are associated with ECMapProb 2.0 simulations, while green and blue symbols refer to BoxMapProb 2.0 simulations. See Supplementary Tables S1-S3 for details

## Discussion

The eruptive chronology and the lithostratigraphic record of a volcanic system are the main source of information to define the expected effects of future eruptions, construct hazard maps, and address the associated risk. This is made by assuming that the system will have similar dynamics in the future (e.g., Gurioli et al. 2010; Calder et al. 2015; Neri et al. 2015b). However, the way the eruptive record is integrated with information obtained from other information sources, such as geophysical studies and numerical modeling, is not straightforward and often hinders reproducibility of hazard assessments. The adoption and extension of some recently published calibration strategies, which are based on documented PDC deposits (Aravena et al. 2022), allowed us to integrate the volcanological record of Tungurahua in a structured and reproducible procedure to define the inputs of a set of numerical simulations, which ultimately resulted in a series of probabilistic, scenario-based PDC hazard maps for this volcano.

A key strength of our results is that the independent use of different comparison metrics, as well as different geological datasets for numerical calibration (when possible), led to uncertainty ranges for the computed inundation

probabilities. Quantitative analysis of uncertainty is typically absent in PDC hazard maps around the world (Lindsay et al. 2023), which hampers representation of the intrinsic variability of the activity observed in volcanic systems. The approach adopted in this investigation complements some recent efforts to quantify the uncertainty associated with PDC hazard assessments in high-risk volcanic systems (e.g., Aravena et al. 2023; Bevilacqua et al. 2017; Neri et al. 2015a; Rutarindwa et al. 2019; Tierz et al. 2018; 2021). In these latter case studies, uncertainty quantification derives from the inclusion of probabilistic vent opening maps from which PDC source positions are sampled (only relevant in case of distributed volcanism) and/or, as in this contribution, from the use of different assumptions to set the inputs of numerical simulations. In our case, this was performed by adopting the following strategies: (1) independent use of different subsets of field data to calibrate the models, (2) independent use of two numerical models, and (3) independent use of multiple metrics to compare field data with the results of calibration simulations. Despite all these sources of uncertainty, the inundation probabilities for a given eruptive scenario of Tungurahua volcano in points of interest from a volcanic hazard perspective are, in general, well constrained (see Supplementary Tables S1-S3),
<figure><img src="image_13.png" /><figcaption>**Fig. 7** Probabilistic hazard maps for the eruptive scenario ES2 (sub-Plinian eruption), considering numerical results of the models ECMapProb 2.0 (top panels) and BoxMapProb 2.0 (middle panels). The mean and maximum probabilistic hazard maps, considering both the models, are included in the bottom panels. Inundation probabilities are indicated by a set of isoprobability curves (see legend) and a rainbow color scale. Labels indicate the main cities. Coordinates are expressed in DD notation</figcaption></figure>

providing clear indications of inundation probability in the context of territorial planning and especially for volcanic crisis management. We stress that the probabilistic analysis of PDC simulations requires the use of numerical models of sufficiently low computational cost to produce thousands of simulations. In fact, several probabilistic PDC hazard

assessments are based on the traditional energy cone model (Clarke et al. 2020; Sandri et al. 2018; Tierz et al. 2016a, b), which presents an extremely low computational cost, but is unable to capture channelization processes that are demonstrably relevant at Tungurahua volcano. The branching formulations (Aravena et al. 2020) may represent a good

<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 8** Mean and maximum probabilistic hazard maps for the eruptive scenario ES3 (sub-Plinian to Plinian eruption), indicated by a set of isoprobability curves (see legend) and a rainbow color scale. These hazard maps were constructed considering P₁ and P₂ (see Fig. 3) as control points for calibration effects. Labels indicate the main cities. Coordinates are expressed in DD notation</figcaption></figure>

compromise between low computational cost and capacity to capture the effect of topography during PDC propagation, which supports the use of ECMapProb 2.0 and BoxMapProb 2.0 in this study.

Adoption of an appropriate way to include uncertainty in hazard assessments depends on the specific characteristics of the studied volcanic system (e.g., monogenetic fields or calderas, where probabilistic vent opening maps are needed, versus stratovolcanoes generally characterized by summit activity only) and the availability of volcanological data to calibrate numerical simulations based on reference eruptions. In the case of Tungurahua volcano, the available volcanological information allows us to characterize reasonably well PDC deposits associated with eruptions with a VEI of 2 or 3 (i.e., ES1 and ES2; Figs. 2 and 3), which is mostly explained due to the recent eruptive cycle of 1999–2016. We stress that, based on field evidence, these eruptive scenarios include most of the PDC-forming eruptions during the last three millennia. On the other hand, we note that the highest intensity scenario (i.e., ES3) was calibrated using a limited number of field data and thus a better knowledge about the eruptive history of Tungurahua is required to quantify the

hazards associated with this eruptive scenario with higher accuracy. In this sense, although the calibrated volumes of collapsing pyroclastic material for eruptive scenarios ES2 and ES3 show overlapping (see Supplementary Figs. S10 and S12), the significant differences between their mean values may suggest that an additional intermediate scenario might be considered to describe the eruptive variability of this volcano. In order to delve deeper into this topic, in Fig. 9, we present the relationship between the collapsing volume of pyroclastic material in a set of non-calibrated, complementary simulations and the resulting runout distance and inundation area of simulated PDCs. In all the cases, which consider different values of sedimentation velocity ($w_s$), an evident break in slope in the modeled inundation area can be recognized at values of collapsing volume of pyroclastic material of about $10^{7.7} \text{ m}^3$ (i.e., $\sim 5 \times 10^7 \text{ m}^3$; see Fig. 9), and a small discontinuity in the slope is observed in the modeled runout distance at collapsing volumes of pyroclasts of about $10^{7.4} \text{ m}^3$ (i.e., $\sim 2.5 \times 10^7 \text{ m}^3$). These results indicate that, in addition to the significant difference in the collapsing volumes of ES2 and ES3, the strong differences in the resulting hazard maps are also modulated

<figure><img src="image_16.png" /><figcaption>**Fig. 9** Runout distance (left-side panels) and inundation area (right-side panels) versus volume of collapsing pyroclasts for a non-calibrated set of Box-MapProb 2.0 simulations</figcaption></figure>
by a change in the behavior of the simulated PDCs above a threshold of collapsing volume of pyroclasts, which in fact coincides with the transition between the calibrated volumes of pyroclasts of ES2 and ES3. The capacity of the topography of stratovolcanoes to influence the behaviour of PDCs has been recently addressed by Aravena and Roche (2022), who classified Tungurahua volcano as a case of intense proximal channelization and moderate distal channelization based on the analysis of numerical simulations of dense PDCs. This is due to the well-defined radial ravines that favor PDC propagation in proximal domains and pronounced tangential valleys (Puela, Chambo, and Pastaza rivers) that buffer the increase of runout distance toward N, NW, and NE and of inundation area when PDCs reach the edifice base (i.e., when they reach the abovementioned tangential valleys). Our numerical results suggest that, above a volume threshold of collapsing pyroclasts, valleys beyond the base of the volcanic edifice become relevant in PDC propagation, strongly affecting the volcanic hazard around Tungurahua. This makes even more critical the need to refine our knowledge about large-scale explosive events at Tungurahua volcano. In addition, Aravena and Roche (2022) recognized the clear effect of proximal obstacles in PDC propagation at Tungurahua, which is probably due to the crater topography and the presence of the ~3 ky BP collapse scar that limit the propagation of small-scale PDCs towards NE.

Although probabilistic volcanic hazard maps as those presented here integrate a large amount of information that are more or less easily understandable by the vast majority of the volcanological community, they are not directly accessible by local communities and decision makers (see for example Thompson et al. 2015). Consequently, the definition of probability thresholds for each scenario in order to translate our probabilistic results into hazard maps with a reduced number of levels is a critical step to use our results in volcanic risk management. However, this is beyond the objectives of this work and should be ultimately defined by Decision Makers and Civil Protection authorities with the collaboration of scientific committees. In order to provide an illustrative example, in Fig. 10, we present a three-color hazard map constructed by extracting the 50% is probability curves of the maximum probabilistic hazard map of each scenario, while the equivalent figures associated with probability thresholds of 10% and 90% are displayed in the supplementary material (Figs. S13 and S14, respectively). When compared to the hazard map presented by Samaniego et al. (2008), the most significant difference, which would be present for any probability threshold adopted for the construction of the three-color hazard map, is that our simulations suggest a non-negligible inundation probability toward SE of Tungurahua volcano for the highest intensity scenario. On the other hand, regarding lower intensity scenarios, both Samaniego et al. (2008) and our results suggest

<figure><img src="image_17.png" /><figcaption>**Fig. 10** Illustrative example of a three-color hazard map of Tungurahua volcano, constructed by considering the 50% is probability curve of the maximum probabilistic hazard map of each scenario. Contours of the hazard levels defined by Samaniego et al. (2008) are also included. Labels indicate the main cities. Coordinates are expressed in DD notation</figcaption></figure>

that the tangential valleys of Chambo, Puela, and Pastaza rivers represent a major limit in the zonification of PDC hazard at Tungurahua. Further comparisons are not possible due to the above-mentioned differences in the definition of the eruptive scenarios.

Regarding limitations of our methodology, we stress that the adopted numerical models, whose reduced computational cost allows for the development of a detailed probabilistic analysis of volcanic hazard, are not able to describe some phenomena that may be relevant during PDC propagation under certain conditions, such as overbanking and detachment processes (Pensa et al. 2018; Trolese et al. 2018). In fact, recent studies show that topography-controlled processes are able to increase significantly the runout distance of PDCs in zones of abrupt slope change and confining conditions (e.g., Charbonnier et al. 2023; Kubo Hutchison and Dufek 2021). The described limitations, as occurs in every numerical modeling-based assessment of natural hazards, make it necessary to adopt conservative choices in the construction of hazard map and reinforce our need to define structured and reproducible strategies to quantify uncertainty and its propagation. On the other hand, we adopted a strategy based on a discrete set of eruption scenarios and, although they were defined considering the geological record of Tungurahua volcano, they were possibly subjected to preservation biases. This strategy is useful for volcanic crisis management because it produces reference maps for different conditions but, in the absence of estimates for the occurrence probability of each scenario, does not

<figure><img src="image_18.png" /></figure>
allow production of long-term hazard scenarios integrating
the entire range of possible PDC scales, as those presented
recently for other volcanic systems (e.g., Neri et al. 2015a;
Sandri et al. 2018). Although a few estimates of the return
period associated with some of the adopted eruption scenar-
ios are present in the literature (Le Pennec et al. 2016; Sam-
aniego et al. 2008), they do not allow for the development
of an integrated analysis to produce long-term hazard maps.

Concluding remarks

We addressed the volcanic hazard associated with PDCs at Tungurahua volcano by adopting an approach based on the development of calibrated numerical simulations for three specific eruptive scenarios, which are defined from our knowledge of the eruptive record of this volcanic system. In particular, we considered violent Strombolian to Vulcanian eruptions (VEI 2), sub-Plinian eruptions (VEI 3), and sub-Plinian to Plinian eruptions (VEI 4–5). The main conclusions associated with this investigation are summarized as follows:

(a) Small-scale PDCs produced during violent Strombolian to Vulcanian eruptions are strongly controlled by at least five dominant channelization ravines towards W, NW, NNW, NNE, and NE, and by crater topography and the ~3 kyr BP collapse scar. The simulated PDCs may reach the Pastaza River through a few ravines and produce inundation probabilities at Baños of about 6±3%.

(b) PDCs generated during VEI 3 sub-Plinian eruptions are also influenced by the proximal topographic features of the volcano (crater morphology and radial ravines). Numerical results indicate that these PDCs frequently reach the Pastaza River through a series of ravines, with inundation probabilities at Baños of 13 ± 4%.

(c) Sub-Plinian to Plinian eruptions (VEI 4–5) produce PDC inundation areas that nearly always involve inhabited centers, including Baños, Puela, and/or Cotaló. However, new volcanological studies to characterize the eruptive history of Tungurahua are required for further constraining the uncertainty affecting this eruptive scenario

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00445-024-01768-y.

**Acknowledgements** Alvaro Aravena and Alessandro Tadini were financed by the French government IDEX-ISITE initiative 16-IDEX-0001 (CAP 20-25). Alessandro was also partly funded by the Cler-Volc project-Programme 1 “Detection and characterization of volcanic plumes and ash clouds” funded by the French government “Laboratory of Excellence” initiative and the CNRS Tellus programme; by

the project “Reti Multiparametriche, Vulcani A7” funded by the Ital-
ian government; and by the program “Convenzione Attuativa per il
potenziamento delle attività di servizio” in the framework of “Accordo
Quadro DPC-INGV 2022-2025”, funded by Dipartimento della Pro-
tezione Civile, Presidenza del Consiglio dei Ministri, Italy. This is
Contribution N° 663 of the ClerVolc Program of Excellence of the
International Research Centre of Disaster Science and Sustainable
Development.

**Funding** Open access funding provided by Istituto Nazionale di Geofisica e Vulcanologia within the CRUI-CARE Agreement.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

**References**

Anderson JF, Johnson JB, Steele AL, Ruiz MC, Brand BD (2018) Diverse eruptive activity revealed by acoustic and electromagnetic observations of the 14 July 2013 intense vulcanian eruption of Tungurahua Volcano, Ecuador. Geophys Res Lett 45(7):2976–2985. https://doi.org/10.1002/2017GL076419

Aravena A, Roche O (2022) Influence of the topography of stratovolcanoes on the propagation and channelization of dense pyroclastic density currents analyzed through numerical simulations. Bull Volcanol 84(7):67. https://doi.org/10.1007/s00445-022-01576-2

Aravena A, Bevilacqua A, Neri A, Gabellini P, Ferrés D, Escobar D, Aiuppa A, Cioni R (2023) Scenario-based probabilistic hazard assessment for explosive events at the San Salvador Volcanic complex, El Salvador. J Volcanol Geoth Res 438:107809. https://doi.org/10.1016/j.jvolgeores.2023.107809

Aravena A, Cioni R, Bevilacqua A, Michieli Vitturi M, Esposti Ongaro T, Neri A (2020) Tree-branching-based enhancement of kinetic energy models for reproducing channelization processes of pyroclastic density currents. Journal of Geophysical Research: Solid Earth, 125(7), e2019JB019271. https://doi.org/10.1029/2019J B019271

Aravena A, Bevilacqua A, de' Michieli Vitturi M, Esposti Ongaro T, Neri A, Cioni R (2022) Calibration strategies of PDC kinetic energy models and their application to the construction of hazard maps. Bulletin of Volcanology, 84, 29. https://doi.org/10.1007/s00445-022-01538-8

Aspden JA, Litherland M (1992) The geology and Mesozoic collisional history of the Cordillera Real. Ecuador. Tectonophysics 205(1-3):187-204. https://doi.org/10.1016/0040-1951(92)90426-7

Bablon M, Quidelleur X, Samaniego P, Le Pennec JL, Lahitte P, Liorzou C et al (2018) Eruptive chronology of Tungurahua volcano (Ecuador) revisited based on new K-Ar ages and geomorphological reconstructions. J Volcanol Geoth Res 357:378-398. https://doi.org/10.1016/j.jvolgeores.2018.05.007

Battaglia J, Hidalgo S, Bernard B, Steele A, Arellano S, Acuña K (2019) Autopsy of an eruptive phase of Tungurahua volcano (Ecuador) through coupling of seismo-acoustic and SO2
recordings with ash characteristics. Earth Planet Sci Lett 511:223-232. https://doi.org/10.1016/j.epsl.2019.01.042

Bernard J, Eychenne J, Le Pennec J-L, Narváez D (2016) Mass budget partitioning during explosive eruptions: insights from the 2006 paroxysm of Tungurahua volcano, Ecuador. Geochem Geophys Geosyst 17:3224-3240. https://doi.org/10.1002/2016GC006431

Bernard J, Kelfoun K, Le Pennec JL, Vallejo Vargas S (2014) Pyroclastic flow erosion and bulking processes: comparing field-based vs. modeling results at Tungurahua volcano, Ecuador. Bull Volcanol 76(9):1-16. https://doi.org/10.1007/s00445-014-0858-y

Bernard B (2018) Rapid hazard assessment of volcanic ballistic projectiles using long-exposure photographs: insights from the 2010 eruptions at Tungurahua volcano, Ecuador. Volcanica 1:49-61. https://doi.org/10.30909/vol.01.01.4961

Bevilacqua A, Neri A, Bisson M, Esposti Ongaro T, Flandoli F, Isaia R et al (2017) The effects of vent location, event scale, and time forecasts on pyroclastic density current hazard maps at Campi Flegrei caldera (Italy). Front Earth Sci 5:72. https://doi.org/10.3389/feart.2017.00072

Bevilacqua A, Aravena A, Neri A, Gutiérrez E, Escobar D, Schliz M, Aiuppa A, Cioni R (2021) Thematic vent opening probability maps and hazard assessment of small-scale pyroclastic density currents in the San Salvador volcanic complex (El Salvador) and Nejapa-Chiltepe volcanic complex (Nicaragua). Nat Hazards Earth Syst Sci 21(5):1639-1665. https://doi.org/10.5194/nhess-21-1639-2021

Bevilacqua A, Aravena A, Aspinall W, Costa A, Mahony S, Neri A, Sparks RSJ, Hill B (2022) Assessing minimum pyroclastic density current mass to impact critical infrastructures: example from Aso caldera (Japan). Nat Hazards Earth Syst Sci 22(10):3329-3348. https://doi.org/10.5194/nhess-22-3329-2022

Bryant JA, Yogodzinski GM, Hall ML, Lewicki JL, Bailey DG (2006) Geochemical constraints on the origin of volcanic rocks from the Andean Northern Volcanic Zone. Ecuador. J Petrol 47(6):1147-1175. https://doi.org/10.1093/petrology/egl006

Bustillos J (2008) Las Avalanchas de Escombros en el sector del volcán Tungurahua. Doctoral dissertation. Escuela Politécnica Nacional: 151

Calder E, Wagner K, Ogburn S (2015) Volcanic hazard maps. Global volcanic hazards and risk, 335-342

Charbonnier SJ, Garin F, Rodríguez LA, Ayala K, Cancel S, Escobar-Wolf R et al (2023) Unravelling the dynamics and hazards of the June 3rd, 2018, pyroclastic density currents at Fuego volcano (Guatemala). J Volcanol Geoth Res 436:107791. https://doi.org/10.1016/j.jvolgeores.2023.107791

Clarke B, Tierz P, Calder E, Yirgu G (2020) Probabilistic volcanic hazard assessment for pyroclastic density currents from pumice cone eruptions at Aluto volcano Ethiopia. Front Earth Sci 8:348. https://doi.org/10.3389/feart.2020.00348

de' Michieli Vitturi M, Esposti Ongaro T, Lari G, Aravena A (2019) IMEX_SfloW2D 1.0: a depth-averaged numerical flow model for pyroclastic avalanches. Geosci Model Dev 12(1):581-595. https://doi.org/10.5194/gmd-12-581-2019

de' Michieli Vitturi M, Esposti Ongaro T, Engwell S, (2023) IMEX_SfloW2D v2: a depth-averaged numerical flow model for volcanic gas-particle flows over complex topographies and water. Geosci Model Dev 16:6309-6336. https://doi.org/10.5194/gmd-16-6309-2023

Douillet GA, Tsang-Hin-Sun È, Kueppers U, Letort J, Pacheco DA, Goldstein F et al (2013) Sedimentology and geomorphology of the deposits from the August 2006 pyroclastic density currents at Tungurahua volcano. Ecuador. Bull Volcanol 75(11):1-21. https://doi.org/10.1007/s00445-013-0765-7

Esposti Ongaro T, Orsucci S, Cornolti F (2016) A fast, calibrated model for pyroclastic density currents kinematics and hazard.

J Volcanol Geoth Res 327:257-272. https://doi.org/10.1016/J.JVOLGEORES.2016.08.002

Eychenne J, Le Pennec JL, Troncoso L, Gouhier M, Nedelec JM (2012) Causes and consequences of bimodal grain-size distribution of tephra fall deposited during the August 2006 Tungurahua eruption (Ecuador). Bull Volcanol 74(1):187-205. https://doi.org/10.1007/S00445-011-0517-5

Falasconi A, Cioni R, Bernard B, Samaniego P, Pistolesi M, Schiavi F (2023) What controls the formation of vulcanian bombs? A case study from the 1 February 2014 eruption of Tungurahua (Ecuador). J Volcanol Geoth Res 444:107961. https://doi.org/10.1016/j.jvolgeores.2023.107961

Fee D, Garces M, Steffke A (2010) Infrasound from Tungurahua volcano 2006-2008: Strombolian to Plinian eruptive activity. J Volcanol Geoth Res 193(1-2):67-81

Ferrés D, Granados HD, Gutiérrez RE, Farraz IA, Hernández EW, Pullinger CR, Escobar CD (2013) Explosive volcanic history and hazard zonation maps of Boquerón Volcano (San Salvador volcanic complex, El Salvador). Geol Soc Am Spec Pap 4(498):201-230. https://doi.org/10.1130/2013.2498(12)

Flynn IT, Ramsey MS (2020) Pyroclastic density current hazard assessment and modeling uncertainties for Fuego Volcano. Guatemala. Remote Sens 12(17):2790. https://doi.org/10.3390/rs12172790

Gaunt HE, Burgisser A, Mothes PA, Browning J, Meredith PG, Criollo E, Bernard B (2020) Triggering of the powerful 14 July 2013 vulcanian explosion at Tungurahua volcano, Ecuador. J Volcanol Geoth Res 392:106762. https://doi.org/10.1016/j.jvolgeores.2019.106762

Gurioli L, Sulpizio R, Cioni R, Sbrana A, Santacroce R, Luperini W, Andronico D (2010) Pyroclastic flow hazard assessment at Somma-Vesuvius based on the geological record. Bull Volcanol 72:1021-1038

Hall ML, Robin C, Beate B, Mothes P, Monzier M (1999) Tungurahua Volcano, Ecuador: structure, eruptive history and hazards. J Volcanol Geoth Res 91(1):1-21. https://doi.org/10.1016/S0377-0273(99)00047-5

Hall ML, Steele AL, Mothes PA, Ruiz MC (2013) Pyroclastic density currents (PDC) of the 16-17 August 2006 eruptions of Tungurahua volcano, Ecuador: geophysical registry and characteristics. J Volcanol Geoth Res 265:78-93. https://doi.org/10.1016/J.JVOLG EORES.2013.08.011

Hall ML, Steele AL, Bernard B, Mothes PA, Vallejo SX, Douillet GA, Ramón PA, Aguaiza SX, Ruiz MC (2015) Sequential plug formation, disintegration by Vulcanian explosions, and the generation of granular Pyroclastic Density Currents at Tungurahua volcano (2013-2014), Ecuador. J Volcanol Geoth Res 306:90-103. https://doi.org/10.1016/j.jvolgeores.2015.09.009

Huppert HE, Simpson JE (1980) The slumping of gravity currents. J Fluid Mech 99(4):785-799. https://doi.org/10.1017/S0022112080000894

Kelfoun K (2017) A two-layer depth-averaged model for both the dilute and the concentrated parts of pyroclastic currents. Journal of Geophysical Research: Solid Earth 122(6):4293-4311. https://doi.org/10.1002/2017JB014013

Kelfoun K, Samaniego P, Palacios P, Barba D (2009) Testing the suitability of frictional behaviour for pyroclastic flow simulation by comparison with a well-constrained eruption at tungurahua volcano (Ecuador). Bull Volcanol 71(9):1057-1075. https://doi.org/10.1007/S00445-009-0286-6

Kubo Hutchison, A., & Dufek, J. (2021). Generation of overspill pyroclastic density currents in sinuous channels. Journal of Geophysical Research: Solid Earth, 126(10), e2021JB022442. https://doi.org/10.1029/2021JB022442

Le Pennec JL, Samaniego P, Eissen JP, Hall M, Molina I, Robin C, Mothes P, Yepes H, Ramón P, Monzier M, Egred J (2005) Los
peligros volcánicos asociados con el Tungurahua. Corporación Editorial Nacional, Quito

Le Pennec JL, Jaya D, Samaniego P, Ramón P, Yánez SM, Egred J, van Der Plicht J (2008) The AD 1300–1700 eruptive periods at Tungurahua volcano, Ecuador, revealed by historical narratives, stratigraphy and radiocarbon dating. *J Volcanol Geoth Res* 176(1):70–81. https://doi.org/10.1016/j.jvolgeores.2008.05.019

Le Pennec JL, De Saulieu G, Samaniego P, Jaya D, Gailler L (2013) A devastating Plinian eruption at Tungurahua volcano reveals formative occupation at ~1100 cal BC in Central Ecuador. Radiocarbon 55(3):1199–1214. https://doi.org/10.1017/S0033822200048116

Le Pennec JL, Ramón P, Robin C, Almeida E (2016) Combining historical and 14C data to assess pyroclastic density current hazards in Baños city near Tungurahua volcano (Ecuador). Quatern Int 394:98–114. https://doi.org/10.1016/j.quaint.2015.06.052

Lindsay JM, Charlton D, Clive MAT, Bertin D, Ogburn S, Wright H et al (2023) The diversity of volcanic hazard maps around the world: insights from map makers. *J Appl Volcanol* 12(1):8. https://doi.org/10.1186/s13617-023-00134-5

Litherland M, Zamora A, Egüez A, Aspden J, Baez N, Beate B et al (1993) Mapa geológico de la República del Ecuador, escala 1:1,000,000. In Geological Survey Publications, Keyworth, UK

Malin MC, Sheridan MF (1982) Computer-assisted mapping of pyroclastic surges. *Science* 217(4560):637–640. https://doi.org/10.1126/SCIENCE.217.4560.637

Marzocchi W, Bebbington MS (2012) Probabilistic eruption forecasting at short and long time scales. *Bull Volcanol* 74:1777–1805. https://doi.org/10.1007/s00445-012-0633-x

Massaro S, Stocchi M, Martínez Montesinos B, Sandri L, Selva J, Sulpizio R et al (2023) Assessing long-term tephra fallout hazard in southern Italy from Neapolitan volcanoes. Nat Hazard 23(6):2289–2311. https://doi.org/10.5194/nhess-23-2289-2023

Mothes PA, Yepes HA, Hall ML, Ramón PA, Steele AL, Ruiz MC (2015) The scientific–community interface over the fifteen-year eruptive episode of Tungurahua Volcano. Ecuador. *J Appl Volcanol* 4(1):1–15. https://doi.org/10.1186/s13617-015-0025-y

Neglia F, Sulpizio R, Dioguardi F, Capra L, Sarocchi D (2021) Shallow-water models for volcanic granular flows: a review of strengths and weaknesses of TITAN2D and FLO2D numerical codes. *J Volcanol Geoth Res* 410:107146. https://doi.org/10.1016/j.jvolgeores.2020.107146

Neri A, Bevilacqua A, Esposti Ongaro T, Isaia R, Aspinall WP, Bisson M, et al (2015a) Quantifying volcanic hazard at Campi Flegrei caldera (Italy) with uncertainty assessment: 2. Pyroclastic density current invasion maps. *J Geophys Res: Solid Earth*, 120(4), 2330–2349. https://doi.org/10.1002/2014JB011776

Neri A, Ongaro TE, Voight B, Widiwijayanti C (2015b) Pyroclastic density current hazards and risk. In Volcanic hazards, risks and disasters (pp. 109–140) Elsevier

Nocquet JM, Villegas-Lanza JC, Chlieh M, Mothes PA, Rolandone F, Jarrin P, Cisneros D, Alvarado A, Audin L, Bondoux L, Martin X, Font Y, Régnier M, Vallée M, Tran T, Beauval C, Maguiña Mendoza JM, Martínez W, Tavera H, Yepes H (2014) Motion of continental slivers and creeping subduction in the northern Andes. Nat Geosci 7(4):287–291. https://doi.org/10.1038/ngeo2099

Orsi G, Di Vito MA, Selva J, Marzocchi W (2009) Long-term forecast of eruption style and size at Campi Flegrei caldera (Italy). *Earth Planet Sci Lett* 287(1–2):265–276. https://doi.org/10.1016/j.epsl.2009.08.013

Palacios PB, Mader HM, Kendall JM, Yepes HA (2023) Seismic tremor location of 10 large paroxysmal eruptions of Tungurahua volcano. Ecuador. *Geophys J Int* 233(2):1460–1483. https://doi.org/10.1093/gji/ggac523

Parra R, Bernard B, Narváez D, Le Pennec JL, Hasselle N, Folch A (2016) Eruption Source Parameters for forecasting ash dispersion

and deposition from vulcanian eruptions at Tungurahua volcano: insights from field data from the July 2013 eruption. *J Volcanol Geoth Res* 309:1–13. https://doi.org/10.1016/JJVOLGEORES.2015.11.001

Parra R, Cadena E, Paz J, Medina D (2020) Isomass and probability maps of ash fallout due to vulcanian eruptions at Tungurahua Volcano (Ecuador) deduced from historical forecasting. Atmosphere 11(8):861. https://doi.org/10.3390/atmos11080861

Pensa A, Capra L, Giordano G, Corrado S (2018) Emplacement temperature estimation of the 2015 dome collapse of Volcán de Colima as key proxy for flow dynamics of confined and unconfinement of pyroclastic density currents. *J Volcanol Geoth Res* 357:321–338. https://doi.org/10.1016/j.jvolgeores.2018.05.010

Rutarindwa R, Spiller ET, Bevilacqua A, Bursik MI, Patra AK (2019) Dynamic probabilistic hazard mapping in the Long Valley Volcanic Region CA: integrating vent opening maps and statistical surrogates of physical models of pyroclastic density currents. *J Geophys Res: Solid Earth* 124:9600–9621. https://doi.org/10.1029/2019JB017352

Samaniego P, Le Pennec JL, Robin C, Hidalgo S (2011) Petrological analysis of the pre-eruptive magmatic process prior to the 2006 explosive eruptions at Tungurahua volcano (Ecuador). *J Volcanol Geoth Res* 199(1–2):69–84. https://doi.org/10.1016/JJVOLGEORES.2010.10.010

Samaniego P, Le Pennec JL, Barba D, Hall M, Robin C, Mothes P, et al (2008) Mapa de los peligros potenciales del volcán Tungurahua (Esc. 1/50.000), 3ra Edición (IGM-IRD). Quito, Ecuador

Sandri L, Costa A, Selva J, Tonini R, Macedonio G, Folch A, Sulpizio R (2016) Beyond eruptive scenarios: assessing tephra fallout hazard from Neapolitan volcanoes. *Sci Rep* 6(1):24271. https://doi.org/10.1038/srep24271

Sandri L, Tierz P, Costa A, Marzocchi W (2018) Probabilistic hazard from pyroclastic density currents in the Neapolitan area (Southern Italy). *J Geophys Res: Solid Earth* 123:3474–3500. https://doi.org/10.1002/2017JB014890

Sobradelo R, Martí J (2010) Bayesian event tree for long-term volcanic hazard assessment: application to Teide-Pico Viejo stratovolcanoes, Tenerife, Canary Islands. *J Geophys Res: Solid Earth* 115(B5):5206. https://doi.org/10.1029/2009JB006566

Tadini A, Roche O, Samaniego P, Guillin A, Azzaoui N, Gouhier M, et al (2020) Quantifying the uncertainty of a coupled plume and tephra dispersal model: PLUME-MOM/HYSPLIT simulations applied to Andean volcanoes. *J Geophys Research: Solid Earth*, 125(2), e2019JB018390. https://doi.org/10.1029/2019JB018390

Tadini A, Bevilacqua A, Neri A, Cioni R, Biagioli G, de' Michieli Vituri M, Esposti Ongaro T, (2021) Reproducing pyroclastic density current deposits of the 79CE eruption of the Somma-Vesuvius volcano using the box-model approach. Solid Earth 12:119–139. https://doi.org/10.5194/se-12-119-2021

Tadini A, Azzaoui N, Roche O, Samaniego P, Bernard B, Bevilacqua A, Hidalgo S, Guillin A, Gouhier M (2022) Tephra fallout probabilistic hazard maps for Cotopaxi and Guagua Pichincha volcanoes (Ecuador) with uncertainty quantification. *J Geophys Res: Solid Earth* 127(2):e2021JB022780. https://doi.org/10.1029/2021JB022780

Tadini A, Azzaoui N, Roche O, Samaniego P, Bernard B, Bevilacqua A, Hidalgo S, Guillin A, Gouhier M (2022) Tephra fallout probabilistic hazard maps for Cotopaxi and Guagua Pichincha volcanoes (Ecuador) with uncertainty quantification. *J Geophys Res: Solid Earth* 127(2):e2021JB022780. https://doi.org/10.1029/2021JB022780

Thompson MA, Lindsay JM, Gaillard JC (2015) The influence of probabilistic volcanic hazard map properties on hazard communication. *J Appl Volcanol* 4:1–24. https://doi.org/10.1186/s13617-015-0023-0

Tierz P, Sandri L, Costa A, Zaccarelli L, Di Vito MA, Sulpizio R, Marzocchi W (2016b) Suitability of energy cone for probabilistic volcanic hazard assessment: validation tests at Somma-Vesuvius
and Campi Flegrei (Italy). Bull Volcanol 78:1–15. https://doi.org/
10.1007/s00445-016-1073-9

Tierz P, Stefanescu ER, Sandri L, Sulpizio R, Valentine GA, Mar-
zocchi W, Patra AK (2018) Towards quantitative volcanic risk
of pyroclastic density currents: probabilistic hazard curves and
maps around Somma-Vesuvius (Italy). J Geophys Res: Solid Earth
123:6299–6317. https://doi.org/10.1029/2017JB015383

Tierz P, Bevilacqua A, Mead S, Spiller E, Sandri L (2021) Field data, models and uncertainty in hazard assessment of pyroclastic density currents and lahars: global perspectives. Front Earth Sci 9:665725. https://doi.org/10.3389/feart.2021.665725

Tierz P, Sandri L, Costa A, Sulpizio R, Zaccarelli L, Di Vito MA, Mar-
zocchi W (2016a). Uncertainty assessment of pyroclastic density
currents at Mount Vesuvius (Italy) simulated through the energy
cone model. Natural hazard uncertainty assessment: Modeling
and Decision Support, 125–145. https://doi.org/10.1002/97811
19028116.ch9

Trolese M, Giordano G, Komorowski JC, Jenkins SF, Baxter PJ, Cho-
lik N et al (2018) Very rapid cooling of the energetic pyroclastic
density currents associated with the 5 November 2010 Merapi
eruption (Indonesia). J Volcanol Geoth Res 358:1–12. https://doi.
org/10.1016/j.jvolgeores.2018.06.004

Vlastélic I, Sainlot N, Samaniego P, Bernard B, Nauret F, Hidalgo
S, Auclair D, Gannoun A (2023) Arc volcano activity driven
by small-scale metasomatism of the magma source. Nat Geosci
16(4):363–370. https://doi.org/10.1038/s41561-023-01143-0

Wadge G, Isaacs MC (1988) Mapping the volcanic hazards from Sou-
friere Hills Volcano, Montserrat, West Indies using an image
processor. J Geol Soc 145(4):541–551. https://doi.org/10.1144/
GSJGS.145.4.0541

<figure><img src="image_19.png" /></figure>