# Probabilistic tephra fallout hazard maps for Sangay volcano, Ecuador

Alessandro Tadini, Federica Pardini, Andrea Bevilacqua, Benjamin Bernard, Pablo Samaniego, Michieli Vitturi, Alvaro Aravena, Silvana Hidalgo, Olivier Roche, Nourddine Azzaoui, et al.

► To cite this version:

Alessandro Tadini, Federica Pardini, Andrea Bevilacqua, Benjamin Bernard, Pablo Samaniego, et al. Probabilistic tephra fallout hazard maps for Sangay volcano, Ecuador. Bulletin of Volcanology, 2025, 87 (10), 10.1007/s00445-025-01794-4. hal-04934013

HAL Id: hal-04934013

https://hal.science/hal-04934013v1

Submitted on 7 Feb 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 Probabilistic tephra fallout hazard maps for Sangay volcano, Ecuador

2 Alessandro Tadini¹,², Federica Pardini¹, Andrea Bevilacqua¹,*, Benjamin Bernard³, Pablo Samaniego², Mattia de'

3 Michieli Vitturi¹, Alvaro Aravena¹,⁴,⁵, Silvana Hidalgo³, Olivier Roche², Nourddine Azzaoui⁶, Augusto Neri¹, Jorge

4 Aguilar³

5 ¹Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Pisa, via Cesare Battisti 53, 56125 Pisa, Italy

6 ²Laboratoire Magmas et Volcans, CNRS, IRD, OPGC, Université Clermont Auvergne, 6 Avenue Blaise Pascal, 63178 Aubière, France

7 ³Instituto Geofisico, Escuela Politécnica Nacional, Ladrón de Guevara E11- 253 y Andalucía, Quito, Ecuador

8 ⁴Facultad de Ciencias Básicas, Universidad Catolica del Maule, Avenida San Miguel 3605, Talca, Chile

9 ⁵Millennium Institute on Volcanic Risk Research – Ckelar Volcanoes, Avenida Angamos 0610, Antofagasta, Chile

10 ⁶Laboratoire de Mathématiques Blaise Pascal, CNRS, Université Clermont Auvergne, 3 Place Vasarely, 63178 Aubière, France

11

12 * Corresponding author: andrea.bevilacqua@ingv.it

13 Abstract

14 Sangay volcano (Ecuador) shows a quasi-continuous activity at least since the 17th century and has produced several 15 eruptions which affected towns and cities at considerable distance (up to > 170 km). For this reason, despite its remote 16 location, recent efforts were aimed at reviewing its volcanic history, quantifying the occurrence probability of four 17 eruptive scenarios of different magnitude (Strong Ash Venting, Violent Strombolian, sub-Plinian and Plinian) and the 18 associated uncertainty, and, for each eruptive scenario, estimating the probability distribution of key eruptive source 19 parameters (fallout volume, average plume height and eruption duration). In this study, we utilize such information to 20 produce probabilistic hazard maps and curves. To this aim, we use coupled plume and dispersal models (PLUME-MOM- 21 TSM and HYSPLIT, respectively) with the application of a novel workflow for running an ensemble of thousands of 22 simulations following a stochastic sampling of input parameters. We produced probabilistic hazard maps for each scenario 23 by considering four ground load thresholds (i.e. 0.1, 1, 10 and 100 kg/m²) and two types of model initialization strategies, 24 based on the elicited total deposit volume and on the elicited plume height, respectively, which produced non-negligible 25 differences. In addition, we produced hazard curves for nine sites of interest from a risk perspective, corresponding to 26 towns/cities potentially affected by tephra accumulation. Finally, we also derived combined maps by merging maps of 27 single scenarios with their probability of occurrence as obtained from expert elicitation. Results indicate that in case of a 28 future eruption, even for a moderate-scale one (Violent Strombolian), probability of tephra accumulation larger than 1 29 kg/m² is relatively high (from 21to 24% considering different model initializations) in the town of Guamote, i.e. the most 30 severely affected site among those tested (43 km W of Sangay). For larger-scale events (i.e. sub-Plinian) the impact of 31 tephra accumulation results to be significant even for the city of Guayaquil (176 km W of Sangay), with probability of 32 tephra accumulation larger than 1 kg/m² from 3% to 22% considering different model initializations. For maps combining 33 single maps of historically observed scenarios, the probability (% - [5ᵗʰ-Mean-95ᵗʰ]) of having ≥ 10 kg/m² for Guamote is 34 [4-13-25] as maximum values.

35 Keywords

36 Sangay volcano; tephra fallout; probabilistic maps; uncertainty quantification; PLUME-MoM-TSM; HYSPLIT

37 Statements and declarations

38 We thank Maurizio Mulas from the ESPOL for kindly providing samples from Guayaquil and Rúben Chérrez and Pablo 39 Ruiz from the ROVE for kindly providing samples from Guaranda and Riobamba. We also thank Marjorie Encalada and 40 Edwin Telenchana for the help in collecting some field data used in this study. We finally thank Silvia Massaro for the 41 editorial handling and two anonymous reviewers for insightful comments that improved the quality of the manuscript.

42 AT was partially funded by CONVENZIONE ATTUATIVA PER IL POTENZIAMENTO DELLE ATTIVITA' DI 43 SERVIZIO A.Q. DPC-INGV 2022-2025. AA was supported by the ANID Iniciativa Milenio-ICN2021_038 "Millenium 44 Institute on Volcanic Risk Research - Ckelar Volcanoes". AT, AA, PS, NA and OR acknowledge the support received 45 from the Agence Nationale de la Recherche of the French government through the program "Investissements d'Avenir" 46 (16-IDEX-0001 CAP 20-25). This work was partially supported by the project "Pianeta Dinamico (Dynamic Planet)- 47 Working Earth"-WUNDERVUL "Toward a Wider Understanding of Volcano", funded by the Italian Ministry of 48 Education, University and Research.
49

50

The authors declare that they have no competing interests. AT, AB and AA are currently guest editor for the Special Issue "Uncertainty Quantification in Volcanology: Observations, Numerical Modelling and Hazard/Risk Assessment".

51

52

Conceptualization: All; Methodology: AT, FP, AB, BB, MdMV, AA, AN; Formal analysis and investigation: AT, FP, AB; Writing - original draft preparation: AT, FP, AB; Writing - review and editing: All; Resources: AT, FP, AB, BB.

53

# 1. Introduction

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

Tephra accumulation on the ground is undoubtedly one of the major hazard sources from volcanic eruptions, due to its impact over large areas (>100 km²; Blong, 1996) and on many aspects of human life and activities (e.g., farmland/livestock, Annen and Wagner, 2003; human health, Baxter and Horwell, 2015; electrical infrastructure, Bebbington et al. 2008; roads/transportation systems, Barsotti et al. 2010; buildings stability, Macedonio and Costa, 2012; water reservoirs and vegetation, Wilson et al. 2012; crops, Ligot et al. 2024). To cope with such hazards, one of the main tools for both short-term volcanic crises management and long-term urban planning is represented by probabilistic hazard maps. Such products rely on Monte Carlo techniques (Hurst and Smith, 2004; Yang et al. 2020) to sample input parameters used to run numerical models of tephra transport and deposition, which are ultimately integrated to create hazard maps (e.g., HAZMAP, Michaud-Dubuy et al. 2021; TEPHRA2, Biass and Bonadonna, 2013; FALL3D, Vázquez et al. 2019; VOL-CALPUFF, Barsotti et al. 2018; ASH3D, Barker et al. 2019; HYSPLIT, Tadini et al. 2022). In the context of probabilistic hazard maps development, the quantification of the major sources of uncertainty is critically important in the estimation of eruptive source parameters (ESPs) and of the probabilities of occurrence of the identified scenarios, conditional on an eruption occurring within defined temporal frames (Bevilacqua et al. 2016; Sandri et al. 2016; Aravena et al. 2023). Examples for volcanic hazard studies that incorporate strategies of uncertainty quantification and propagation include comparison between available field data and existing global volcanological databases (Biass and Bonadonna, 2013), how large is the uncertainty in estimating input parameters from field data (e.g., Biass and Bonadonna, 2011; Constantinescu et al. 2022; Yang and Jenkins, 2023), how the input uncertainty propagates through the model (e.g., Woodhouse et al. 2015; de' Michieli Vitturi et al. 2016; Tadini et al. 2020) and how it is influenced by wind field variability (e.g., Macedonio et al. 2016; Madankan et al. 2014; Stefanescu et al. 2014). Among the various techniques to quantify uncertainty in the above-mentioned situations, structured expert judgment has revealed to be particularly useful in case few volcanological information is available (see e.g. Neri et al. 2008; 2015; Tadini et al. 2017b; 2021; Bevilacqua et al. 2015; 2022).

76

77

78

79

80

81

82

83

84

85

86

In this paper we present probabilistic tephra fallout hazard maps for Sangay volcano (Ecuador), located ~200 km south of Ecuador's capital city Quito and 175 km to the east of Guayaquil, the main maritime port and commercial capital of the country (Fig. 1). Because of its location, extreme humid weather conditions, and access limited to a difficult two-day trail, Sangay's geology and recent activity are still poorly known (Monzier et al. 1999). For this reason, Sangay volcano has been the topic of a companion study (Bernard et al. 2024) that includes characterization of its eruptive activity and quantification (through an expert elicitation) of the occurrence probability of different eruptive scenarios as well as the probability distribution of their ESPs. The approach described here, similarly to Tadini et al. (2022), utilizes the tephra transport and dispersal model HYSPLIT (Stein et al. 2015) to run a set of simulations whose results are then used to produce hazard maps. To improve the representation of the conditions at the volcanic source, we initialized HYSPLIT with the results of the eruption column model PLUME-MoM-TSM (de' Michieli Vitturi and Pardini, 2021) through a novel workflow that incorporates uncertainty quantification and propagation by using two different model initializations.

87

88

89

90

91

This study is organized as follows. We first summarize the results of Bernard et al. (2024) on Sangay volcano eruptive history and expert elicitation findings in Section 2. Then, in Section 3 we describe the modeling strategy, the uncertainty quantification for ESPs and occurrence probability of eruptive scenarios, and the procedure for the production of hazard maps and hazard curves. Finally, Sections 4 and 5 focus, respectively, on the presentation of the hazard maps/curves and on the discussion of their main implications, while section 6 draws some final considerations.

92

93

# 2. Background

94

## 2.1 Sangay volcano eruptive activity

95

96

97

Bernard et al. (2024) classified the activity of Sangay in quiescence, weak eruptive activity, enhanced eruptive activity, and eruptive pulses. This classification is based on the direct observation of the last ~25 years and historical accounts of the last ~400 years.
*Quiescence* is defined as a state of calm or inactivity for a relatively prolonged period and occurred six times in the 2000-2024 interval (Vasconez et al. 2022). These periods lasted between 5 and 20 months and during them few isolated ash emissions occurred, whose products impacted only the immediate vicinity (100's m) of the crater and fumarole fields.

*Weak eruptive activity* is instead the most common behavior of Sangay volcano during the 2000-2024 interval. Bernard et al. (2024) identified two types of weak eruptive activity, namely long-lasting periods (up to several years) and short episodes (few weeks to months). The formers are characterized by a low rate of seismic events (~8-18 per day), a low lava emission rate (~0.3 m³/s), and SO₂ degassing undetectable by satellite. Short episodes are instead characterized by intermediate seismic event rates (~72-88 per day), intermediate lava emission rates (~0.8-1.5 m³/s), and low SO₂ degassing (up to 224 t/day). Both are classified as weak eruptive activity because of the absence of impact outside the volcanic cone.

*Enhanced eruptive activity* has been consistently observed in the 2000-2024 interval and has manifested itself as an increase (from weak activity) in lava emission rate (average 5 m³/s), in the frequency of explosions (up to 1,000 events per day) and, more consistently, in SO₂ degassing (up to 4,500 t/day). Enhanced eruptive activity produces frequent ash fallout in local communities, as well as increased sedimentation in rivers flowing down from the volcano.

Finally, *eruptive pulses*, which occur at a much shorter time scale than the previous types of activity, are predominantly explosive and have been classified by Bernard et al. (2024) in four categories (from small to large size): Strong Ash Venting (*SAV*), Violent Strombolian (*VS*), sub-Plinian (*SPL*) and Plinian (*PL*). Among them, only the first two categories have been observed over the past 20 years. *SAV* events have frequently been associated with peaks of lava emission rate and PDCs formation (Vasconez et al. 2022; Bernard et al. 2024; Hidalgo et al. 2022) and lasted for hours to days with bent column height < 5 km above crater and sometimes higher (up to 8 km above crater) pulsatile gas plumes. During *SAV* events, ash clouds frequently reached more than 100 km from Sangay (Moran-Zuloaga et al. 2023). *VS* in the past 20 years only occurred between September 2020 and May 2021 and are easily spotted on both local and regional seismic records (Bernard et al. 2022; Moran-Zuloaga et al. 2023; Bernard et al. 2024). They have been associated with large SO₂ emissions (up to >45,000 t/day), high lava emission rate (>20 m³/s) and PDCs formation. According to the monitoring data, *VS* events lasted between 2 and 6 hours with eruptive plumes that generated column height of 5-10 km above the crater which created wide umbrella-shaped ash clouds. One *SPL* has occurred during the past 400 years (in 1628 CE) and one Holocene *PL* pulse has been recognized based on a 2-20 cm-thick dacitic tephra fallout deposit, 1 meter beneath the current soil, found up to ~20 km from the volcano (see Bernard et al. 2024 for more details).

Given the few constraints on the occurrence of *PL* eruptions at Sangay, in the following sections we first present the maps and the analyses for the *SAV*, *VS* and *SPL* scenarios only, which are the eruption scenarios historically observed over the past ~20 (SAV and VS) and ~400 (SPL) years. In the discussion we also present the results when considering the *PL* scenario.

## 2.2 Expert elicitation for eruptive scenarios at Sangay volcano

To quantify the uncertainty associated with the eruptive scenarios of Sangay volcano (described in the previous section), a performance-based expert elicitation session was conducted in 2021 (see Aspinall, 2006; Aspinall and Cooke, 2013 for generalities on expert elicitation and Bernard et al. 2024 for more details on the Sangay elicitation). Eighteen experts took part in this exercise and were asked to give their judgements on:

* the probability, for each scenario, of Sangay experiencing at least one eruption in two different timeframes, namely 10 years and 100 years in the future;

* the probability distribution, for each scenario, of three eruptive source parameters, namely the average height of the column during the pulse (in km above the crater), the total volume of the fallout deposit emitted during the pulse (in 10⁶ m³) considering an equivalent deposit density of 1000 kg/m³ and the average duration of the eruption as detected by geophysical monitoring (in hours).

For all these questions and after carefully considering all the evidence and data presented during the meeting, experts provided their judgments as the 5th/Median/95th percentile values of a probability distribution. Individual experts' answers for target questions were then pooled together by using an equal weights rule, and two independent performance-based models (Classical Model, Cooke, 1991; ERF, Flandoli et al. 2011). In the following we adopt the results based on the Classical Model; for more information on the other approaches see Bernard et al. (2024). To derive the pooled probability density functions for each question (representing the view of the group in the form of a new virtual expert, called "Decision Maker"), individual experts' answers were sampled 10000 times (de' Michieli Vitturi et al. 2024).
148 Results of the probability for the four scenarios for the next 100 years (utilized for map combination, see Section 3.3) are displayed in Fig.2. Note that we plot probability density functions of occurrence probabilities, i.e. uncertainty distributions of probability estimates, according to a doubly stochastic approach (Tadini et al. 2017a; Bevilacqua et al. 2018; 2021).
149
150 Results of the three eruptive source parameters for the four scenarios (utilized in the “Simulation setting” section) are instead displayed in Figure 3.
151
152

153

154 **3. Methods**

155 *3.1 Numerical models and workflow for map production*

156 The HYSPLIT (Hybrid Single-Particle Lagrangian Integrated Trajectory) model was used to simulate the atmospheric dispersion of tephra particles and their deposition on the ground (Stein et al. 2015). HYSPLIT has been extensively utilized for both research and operational purposes, with some Volcanic Ash Advisory Centers relying on it as a standard model (Mastin et al. 2017). In order to execute the numerical simulations, the model requires inputs on meteorological conditions and on the tephra injection rates in the atmosphere. Meteorological conditions are supplied as a dataset where the parameters describing the evolution of the atmospheric properties in space and time are listed on a 3D grid. Wind speed and direction are the main parameters controlling the transport and deposition on the ground of tephra particles. Tephra injection conditions, describing how tephra particles are released into the atmosphere (including factors such as amount, duration, and vertical distribution), also need to be defined as inputs. To this aim, we coupled HYSPLIT with the 1D steady-state eruption column model PLUME-MoM-TSM (two-size moment - de' Michieli Vitturi and Pardini, 2021). PLUME-MoM-TSM simulates the ascent of an eruption mixture through the atmosphere from the vent source until it reaches its maximum height, accounting for wind effects and simulating the umbrella cloud spreading at the neutral buoyancy height. Tephra particles lost from the column below the neutral buoyancy height are injected into the atmosphere from pointwise sources located on the plume axis. Particles reaching the neutral buoyancy enter the atmosphere through a circular area source which best fits the position and upwind spreading of the umbrella cloud. Once in the atmosphere, particle transport and deposition on the ground is simulated by HYSPLIT based on the assumed particle properties and the characteristics of the wind field. The coupling between PLUME-MoM-TSM and HYSPLIT is done through Python scripts that allow the results of PLUME-MoM-TSM to be automatically used to generate the input file for HYSPLIT. A detailed description of this procedure, together with several applications, can be found in Pardini et al. (2020) and Tadini et al. (2020). We note that PLUME-MoM-TSM can be initialized with either mass flow rate or the maximum column height, in both cases by taking into account atmospheric conditions. In the latter case, an inversion step is performed to retrieve the conditions at the vent (initial velocity and plume radius) that correspond to a mass flow rate capable of generating a column with the desired height.
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178

179 In this study, we applied the above-mentioned models to produce probabilistic maps of tephra accumulation on the ground for Sangay volcano. To achieve this, we developed a numerical workflow that automatically runs an ensemble of thousands of simulations, i.e. a Monte Carlo simulation, where each ensemble member is initialized with a specific set of ESPs and meteorological conditions. Eruption conditions include mass eruption rate (or column height), eruption duration, and physical properties of the tephra particles, such as size, shape and density. ESPs of each ensemble member can be either fixed values common to all the members or random values sampled from pre-defined probability density functions. Indeed, some eruption parameters are more difficult to constrain than others and/or have a higher impact on model results. These considerations guided the selection of the parameters that could be treated as fixed and equal among all the ensemble members, as well as those best defined by using probability density distributions to assign different values to each ensemble member. Probability density distributions were derived from a literature review and/or from a Monte Carlo sampling of the results of the expert elicitation described above. Finally, to account for the uncertainty on meteorological conditions, the eruption onset time of each ensemble member was randomly sampled in the period 2012-2022, taking into account seasonal variability by having the same number of samples within each month of the considered period. The assumption behind this choice is that the last decade of meteorological conditions is appropriate to statistically capture the main features of atmospheric circulation in the coming decades. In addition, our sampling time window fully overlaps with those of recent similar studies, which consider time windows ranging from 6 to 29 years (Bonasia et al. 2011, 2014; Titos et al. 2022).
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195

196 Once the appropriate set of initial conditions are defined for each ensemble member, the numerical simulations are executed and processed to produce the outcomes of interest which, in our case, are probabilistic hazard maps and curves (see Section 3.3).
197
198

199
## 3.2 Simulation setting

The domain considered in the tephra hazard assessment has a vertical extension of 50 km above ground and covers the entire Ecuador, spanning from 5.5°S to 1.5°N in latitude and from 81.0°W to 75.0°W in longitude, as shown in Fig. 1. The numerical simulations in HYSPLIT were performed on a grid with horizontal resolution of 0.05° (approximately 5 km). The adopted meteorological information derives from the ERA5 dataset (Hersbach et al. 2020) supplied by the European Centre for Medium-Range Weather Forecasts (ECMWF). ERA5 is a reanalysis dataset for the global weather with data available from 1940. Data is gridded hourly with a horizontal resolution of 0.25° x 0.25° (approximately 27-28 km in latitude and longitude at the equator) and 37 vertical pressure levels (from 1000 to 1 hPa). A dataset covering 10 years from 2012 to 2022 was used to perform the numerical simulations, under the assumption that such dataset is appropriate for evaluating the hazard due to future eruptions at Sangay volcano. As a preliminary analysis, wind rose diagrams (for ERA5 data) showing wind speed and direction at the Sangay volcano location were constructed for different height intervals from the topographic surface to 50 km (Fig. 4). It emerges that the wind mainly blows toward the west for all the considered height intervals. The greatest spread in wind direction is observed from 10 km to 20 km, where wind appears to blow in sectors NW-W and W-SW, with a lesser extent toward the east. A strong dominance of the winds blowing toward the west is clearly observed at the highest heights (20-35 km and 35-50 km), with a smaller extent observed in the eastward direction. For these latter two heights, a seasonal analysis for, respectively, dry season (Oct-May) and wet season (Jun-Sep) does not evidence a preferential seasonal dominance for the 20-35 km range, while for the 35-50 km range, stronger winds tend to blow more frequently toward W in the wet season, although in the dry season there is still a prevalence of W-blowing winds (see Fig. 1S from Supporting Information 2). Since typical neutral buoyancy levels for our four eruption types are located between 10 and 35-40 km (see e.g. Mastin et al. 2009), we provide wind pattern of the 1200 simulations for the SAV and PL eruption types in Figs. 2S-3S for the 10-20, 20-35 and 35-50 km ranges. From the analysis of these diagrams (results are the same also for the VS and SPL eruptions), we could conclude that our samples fully capture the variability of the wind pattern of the whole dataset.

The definition of the eruption conditions of the four investigated scenarios (SAV, VS, SPL, and PL) was based on both literature review and expert elicitation. The parameters defined through probability density distributions are: eruption duration, total mass of the deposit, plume height, and Total Grain Size Distribution (TGSD). The remaining parameters were fixed for all the ensemble members, and a complete list can be found in Supporting Information 1. The uncertainty on the TGSD was based on literature review. In particular, TGSD was assumed independent from eruption size (given the few constraints available for TGSD at Sangay, especially from proximal deposits - see Bernard et al. 2022) and as a lognormal distribution defined by a mean value μ and a standard deviation σ. Both μ and σ were considered uncertain and defined as uniform distributions in the ranges [-1; 3] φ for μ and [2; 3] φ for σ. Eruption duration, total mass of the deposit, and plume height were supplied by the elicited distributions in Bernard et al. (2024).

An initial set of 10000 samples for each eruption scenario was produced by independently sampling the distributions of eruption duration, total volume of the deposit (converted into mass considering 1000 kg/m² as deposit density), plume height, μ and σ. Since mass flow rate was not directly supplied by the elicitation, for each of the 10000 samples we computed the corresponding average mass flow rate by converting the total deposit volume into total mass (with an additional 5% mass to account for magma water content) and dividing this latter by the duration of the eruption.

As stated above, the eruption column model can be initialized with either total volume or column height. Therefore, for each scenario, we performed two sets of numerical simulations using either one of the initialization schemes. In order to have a good compromise between statistically meaningful outputs and acceptable computational times, for each eruption scenario and initialization scheme we performed 1200 simulations by running PLUME-MoM-TSM and HYSPLIT. The eruption parameters used to initialize the 1200 simulations are a subset of the original 10000 samples, and a preliminary analysis was done to ensure consistency between the parameters values used to initialize the models and the results of the elicitation. This analysis consisted in running PLUME-MoM-TSM only (without HYSPLIT) over the entire set of 10000 samples and rejecting the samples whose numerical results were inconsistent with the results of the elicitation.

In particular, for the simulations initialized with the total volume (tv), we rejected the simulations giving as output a value of plume height outside the range of values of the corresponding elicited distribution. Similarly, for the simulations initialized with plume height (ph), we rejected those producing a total volume erupted outside the range of the elicited distribution. We remark that we did this analysis for each eruption scenario and initialization scheme.

Figure 3 illustrates the comparison between the original distributions formed by the 10000 samples (red dashed lines) and the filtered distributions of the “accepted” samples (blue dashed lines). The 1200 simulations are randomly sampled from these latter distributions. For the simulations initialized with tv (left panels in Fig. 3), 2912 and 5058 samples passed the check, respectively, in the SAV and VS scenarios. In particular, in the SAV scenario the majority of samples has a mass
flow rate (and a total deposit volume) lower than those consistent with the elicited plume height. In the SPL and PL scenarios, a higher proportion of samples passed the check, with 8167 samples in the SPL scenario and 9129 in the PL scenario (for PL, see Fig. 19S from Supporting Information 2). These results indicate a good agreement between the elicited deposit volumes and plume heights in these scenarios. Regarding the simulations initialized with ph (right panels in Fig. 3), it emerged that ca. 80-90% of samples passed the check for the SAV (8964) and VS (8385) scenarios, while 5331 samples were accepted in the SPL scenario. In the PL scenario (see Fig. 19S from Supporting information 2), 8686 samples passed the filtering step. This filtering approach mitigated the independence of the 10000 samples of eruptive source parameters, especially between mass flow rate (derived from total volume) and column height and, at a lower degree, with respect to the previous two parameters and eruption duration. The only purely independent parameter is TGSD, whose effect in the final results has been assumed secondary with respect to the other inputs. Future studies with more quantitative data on TGSD for Sangay volcano could help in updating our results.

To associate meteorological conditions with each ensemble member, the start date and time of the eruptions were randomly sampled from 2012 to 2022 (the start and end years of the ERA5 dataset). The end dates of the eruptions were derived from their durations, and the numerical runs continued for additional 24 hours after the end of each eruption to achieve complete tephra deposition on the ground. To guarantee a uniform representation across all months, the meteorological dataset was sampled to initialize 100 simulations each month, totaling 1200 simulations.

269

### 3.3 Type of outputs

The outcomes that we produced by running a Monte Carlo simulation consisting of an ensemble of 1200 members are probabilistic hazard maps and curves for tephra accumulation on the ground.

Hazard maps cover the entire domain considered in the simulations and show the probability of exceeding specific thresholds of ground load. They are obtained by calculating the probability $P_{\tau}$ that the ground load at each pixel $N_{x,y}$ of the computational grid exceeds a certain ground load threshold ($\tau$). $P_{\tau}$ is determined by counting the number of times ($n_i$) $\tau$ is exceeded by the ground load resulting from each simulation of the ensemble ($GL_i$), and then dividing by the total number $m$ of simulations:

$$P_{\tau}(N_{x,y}) = \frac{\sum_{i=1}^{m} n_i}{m} \qquad [1]$$

where

1, if $[GL_i(N_{x,y}) \ge \tau]$ [2]

$n_i =$ 0, otherwise [3]

In this study, in order to obtain a set of results able to provide useful information about the potential impacts (Jenkins et al. 2015), we adopted four threshold values: 0.1, 1, 10 and 100 kg/m². These thresholds were chosen as they can be associated to various degrees of damage and disruption to communities and particularly:

* 0.1 kg/m² is sufficient enough to be detectable by residents and might cause minor vegetation damage and the covering of road markings (Magill et al. 2013);

* 1 kg/m² may lead to airport closure, contamination of water supplies, damage of electrical appliances, disruption of road, rail and air transport, and respiratory problems (Jenkins et al. 2012);

* 10 kg/m² cause minor structural damage for weaker building and major damage on crop and vegetation (Jenkins et al. 2012);

* 100 kg/m² may cause partial or total collapse of several buildings (Jenkins et al. 2012);

Hazard curves, instead, are local results computed at specific locations of interest. For each investigated location, the corresponding hazard curve indicates the exceedance probability $EP = Pr(GL > X)$, where GL is the ground load computed from the simulations and X is a specific value within a certain interval. To derive such curves, outputs of ground load
298 GLi resulting from the m simulations are sorted in ascending order (GL₀, GL₁, GL₂,..., GLₘ-₁), while the exceedance
299 probability EPi is computed as follows (following Bonadonna 2006):

$$
300 \quad EP_i = 1 - \frac{i}{m}, \quad \text{for } 0 \le i < m \hfill [4]
$$

301 The locations chosen to compute the hazard curves are those reported in Figure 1 and include the capital city of Ecuador,
302 Quito, as well as some main cities/towns located within 200 km from Sangay volcano.

303 All maps and curves are produced for either single scenarios and combined scenarios. In this latter case, we consider the
304 uncertainty distribution of the probability of having at least one eruption among each scenario in the next 100 years (Fig.
305 2); the resulting maps/curves are related to the occurrence, for each scenario, of the largest event within the next 100
306 years, excluding replication of the same event and cumulative tephra load by multiple events. The combination of maps
307 and curves was done, respectively, for the SAV and VS, for the SAV, VS and SPL and for the SAV, VS, SPL and PL
308 joint scenarios. For example, to obtain the hazard maps of the SAV-VS-SPL-PL (the others are similar in the approach),
309 we start from the probabilistic maps (H) of the four scenarios and the 10000 samples defining the distributions of
310 probabilities (P) of SAV, VS, SPL and PL for the next 100 yrs (Fig. 2). For each sample, we combine the four scenarios
311 as follows:

$$
312 \quad H_{\text{comb}} = P_{\text{SAV}} * (1-P_{\text{VS}}) * (1-P_{\text{SPL}}) * (1-P_{\text{PL}}) * H_{\text{SAV}} + P_{\text{VS}} * (1-P_{\text{SPL}}) * (1-P_{\text{PL}}) * H_{\text{VS}} + P_{\text{SPL}} * (1-P_{\text{PL}}) * H_{\text{SPL}} + P_{\text{PL}} * H_{\text{PL}} \quad [5]
$$

313 By applying equation [5] for each of the 10000 samples defining the distributions of probabilities, we obtain 10000
314 maps/curves and from these latter we show those representing the mean value and the 5th/95th percentiles to illustrate the
315 uncertainty associated with the map combination. We then applied the same approach to derive the combined hazard
316 curves (representing the mean value and the 5th/95th percentiles) for each location.

317

318 **4. Results**

319 In this section we present, in the single eruptive scenarios, the hazard maps for the 1 kg/m² threshold and the hazard
320 curves for all sites. Then we present, in the combined SAV-VS-SPL and SAV-VS-SPL-PL scenarios, hazard maps for
321 the 1 kg/m² threshold and hazard curves for Guamote, Guayaquil and Macas. The remaining maps (0.1, 10 and 100 kg/m²)
322 and curves (combined scenarios for the remaining six sites) are available in Supporting Information 2. Moreover,
323 interpolated values at 0.1, 1, 10, 100 kg/m² for all the hazard curves are available in Supporting Information 3. We also
324 provide ASCII Arc/Info files for all maps in Supporting Information 4.

325

326 **4.1 Hazard maps/curves in the single scenarios (SAV, VS, SPL)**

327 By looking at the maps (Fig. 5), the curves (Fig. 6) and the interpolated values of the hazard curves (Supporting
328 Information 3) for single scenarios, we note that, between the two initialization setups, differences tend to be higher with
329 increasing the size of the eruption (i.e., larger in the SPL than in the SAV/VS). For example, considering the 1 kg/m² and
330 the interpolated probabilities of hazard curves, we observe that differences between initialization setups at all sites are of
331 maximum ~12% in the SAV/VS, while they could reach up to ~40% for the SPL. Regarding the sites for which we
332 produced hazard curves (Fig. 6), Guamote (located ~43 km W of Sangay) is the town that, on average, would suffer more
333 tephra accumulation (probability of tephra accumulation of 1 kg/m² between 4.6% and 78% in the three scenarios, see
334 Supporting Information 3). We stress that at Guamote, even in case of a moderate-size eruption such as a VS, the
335 probability of having a ground load of 1 kg/m² is relatively high (i.e. 21-24%). For the same eruption type, the town of
336 Alausí (~60 km SW of Sangay) would also be affected with similar probabilities (16%-18%), while the two largest cities
337 in Ecuador (Quito and Guayaquil) would suffer a potential tephra accumulation of 1 kg/m² mostly for a SPL eruption
338 type with, respectively, 0.7-6% and 3-22% probability.

339

340 **4.2 Hazard maps/curves for combined scenarios (SAV, VS, SPL)**

341 Figures 7-8 show, respectively, the hazard maps of the 1 kg/m² threshold and the hazard curves (for Guamote, Guayaquil
342 and Macas towns) related to the combined SAV, VS and SPL, which are the eruptive styles actually observed over the
343 last 20 years (SAV and VS) and 400 years (SPL), respectively (for the 0.1, 10 and 100 kg/m² thresholds please refer to
344 Figs. 7S-9S from Supporting Information 2; hazard curves for the remaining 6 sites are reported in Figs. 10S-11S from
345

346

347

348

349

350

Supporting Information 2). For completeness, we also report in Supporting Information 2 (Figs. 12S-18S) hazard maps and curves of combined smallest-scale styles (SAV and VS). Even in this case (similarly to section 4.1), we note that differences between the different adopted initialization setups could be relevant (reaching up to 40% from the sites where hazard curves were produced, see Supporting Information 3).

### 4.3 Combined maps when considering the Plinian (PL) event

351

352

353

354

355

356

357

358

359

In the previous section, we have presented maps of the combined SAV, VS and SPL scenarios, as they represent eruptions that occurred at historical and recent times. However, the PL scenario has also been considered by Bernard et al. (2024) due to the presence of a tephra deposit associated with the Holocene eruptive history of Sangay that is consistent with a Plinian event. The expert elicitation outcome assigned to this scenario a median probability of occurrence within the next 100 years of ~1%, but with an upper uncertainty bound (95th percentile) of ~20%. In Figure 9 we report both the hazard maps (1 kg/m² threshold, for the other thresholds see Fig. 18S from Supporting Information 2) and the hazard curves of this scenario for the two initialization setups, where we can easily observe that differences between the two setups are sometimes high (i.e. up to ~40% when comparing interpolated probabilities from hazard curves, see Supporting Information 3).

360

361

362

363

364

365

366

367

368

369

Figures 10 and 11 display the maps (1 kg/m²) and curves (Guamote Guayaquil, Macas) of the four combined scenarios (for other thresholds and sites please refer to Figs. 21S-25S from Supporting Information 2). The maps/curves with only three scenarios (Figs. 7-8) exhibit larger differences for more distal location with respect to vent area and proximal locations (see for instance the hazard curve for the town of Guamote). Differences are also more pronounced between the maps/curves with the initialization tv (when compared with those produced with an initialization based on ph) and between maps/curves referred to the 95th percentile (as compared to the Mean and the 5th percentile). For example, for the town of Guamote and the 1 kg/m² threshold (see Supporting Information 3), the difference between the 95th percentile probabilities of the tv curves is ~18%, while differences between the Means and the 5th percentiles are, respectively, ~8% and ~3%. Conversely, for the ph curves, differences for the 5th percentiles, Means and 95th percentiles are, respectively, ~3%, ~4% and ~6%. Such differences are proportionally similar for the other sites.

370

371

372

373

374

Average differences are around ~6% for all sites, thresholds, percentiles and initialization setups, while a maximum difference of ~28% is computed for the Ambato site, tv initialization, threshold 0.1 kg/m² and 95th percentile. For comparison, average differences between the SAV-VS-SPL and the SAV-VS interpolated values is ~9% for all sites, thresholds, percentiles and initialization setups, while the maximum difference is ~40% for Macas site, ph initialization, threshold 0.1 kg/m² and 95th percentile.

375

376

## 5. Discussion

377

### 5.1 Comparison of hazard maps with historical observations and existing hazard maps

378

379

380

381

382

383

384

The results presented in the previous sections can be compared both to observed mass loadings at key locations during recent eruptions and to existing hazard maps. With respect to the former, Bernard et al. (2022) provided a detailed description of the VS eruptive event of September 20, 2022. During this pulse, ground load was measured in some sites for which we have produced hazard maps, e.g. Guamote, Alausí and Guayaquil. Values for these three locations are, respectively, 0.19, 0.175 and 0.02 kg/m²: in the same locations, hazard curves of the VS eruptive pulse (Fig. 7) showed probability of exceeding the same values of ~50%, ~45% and ~20%. Our results are therefore in accordance with these observations.

385

386

387

388

389

390

391

392

393

394

To further compare our results with actual observations, we provide some new data about ground load measured at 5 towns (Alausí, Guamote, Guaranda, Guayaquil, Riobamba, see Fig. 1) during 12 eruptive pulses (8 SAV and 4 VS) occurred between June 2020 and June 2024 (Table 1). These data, although based on a relatively small number of events, allow to count the number of eruptions which caused, in each location, a ground load larger than 0.1 and 1 kg/m², and to compare them with a binomial model using the probability of exceeding the same thresholds conditioned on the occurrence of either SAV or VS eruptions for the two initializations (mfr and ph) as interpolated from the hazard curves (see Supporting Information 3). This comparison is presented in Table 2, and the binomial likelihoods of the best performing initializations in each location are all above 25%, with peaks above 95% in several cases. These statistics are not enough for validating our approach, but we could appreciate that our modeling results associate significant likelihood to the real observations.
395
396
397
398

We remind, however, that for the smallest tephra accumulation values (0.1 and 1 kg/m²), corresponding to about 0.1, 1 mm of deposit, the ash is very easily remobilized by the local winds, producing large spatial variations. In this context, these tephra loading can be considered only as average values valid on a relatively large area, with local variations that may reach one or more orders of magnitude.

399
400
401
402
403
404
405
406
407
408
409

The current official IG-EPN hazard maps for Sangay volcano has been developed by Ordoñez et al. (2011). In this map, two areas potentially affected by tephra fallout were identified with respect to two scenarios (“small” and “moderate”, roughly corresponding to SAV and VS, respectively) and modeled with fixed ESPSs and random wind directions/velocities. Although a detailed comparison with respect to our fully probabilistic results is not straightforward, we note that in the official IG-EPN map, the area potentially affected by tephra fallout for moderate eruption at Sangay encloses the town of Guamote, which is the one which would suffer most of tephra accumulation according to our maps/curves. The area impacted by tephra fallout for small eruptions (SAV) is instead more restricted in areas proximal to the volcano excluding Guamote, while our maps/curves indicate a probability of exceeding 1 kg/m² of 9-13% even in case of SAV eruptions. In addition, we also note that, in Ordoñez et al. (2011), the towns of Alausí and Riobamba are not enclosed neither within the “moderate” eruption scenario (VS) impact area, while indeed they has been affected, respectively, three times and one time since 2020 by tephra accumulation > 0.1 kg/m² (see Table 1).

410
411

412

## 5.2 Differences between the two initialization setups

413
414
415
416
417

The differences between the two initialization setups ph and tv are already anticipated from the input parameter distributions (see section 3.2) and are more evident for larger magnitude scenarios (i.e. SPL and PL). In the SPL/PL scenarios, and especially in the PL scenario, the elicited distribution of average plume height is not always consistent to the elicited fallout total volume distribution, and in these scenarios the ph setup produces a significantly greater fallout hazard than the tv setup (see Fig. 19S from Supporting Information 2).

418
419
420
421
422
423
424
425

We remind that when we had to calculate the mass flow rate, we obtained it from the elicited total volume (converted into mass considering an equivalent density of 1000 kg/m²) and duration, and averaged it over the course of the simulated eruption. The high uncertainty in estimating tephra fallout volume/mass, as well as in eruption duration, contributed to the large uncertainty on the mass flow rate, also considering that total volume and duration have been sampled independently in our simulations. Nevertheless, we note that an average plume height that reaches high elevation, not only implies a greater total volume than a lower plume with the same duration, but also allows ash to be transported far from the vent area by the strongest directional winds (see also Fig. 4). This may have enhanced the effect of the uncertainty affecting the plume height on the fallout hazard assessments, if compared to the total volume or the eruption duration.

426
427
428
429
430
431

In particular, we note that the elicitation in Bernard et al. (2024) estimated the total deposit volume, which may underestimate the total erupted volume significantly because of the very fine ash fraction, which does not settle for long periods and can range from a negligible amount to more than 30% of the erupted mass (Rose and Durant, 2009). Conversely, the plume height elicitation was considering variable wind conditions, but the wind field in the simulation was sampled independently from the plume height value. Consequently, the association of high plumes and strong winds may have produced an over-estimation of the erupted volume in the second approach.

432
433
434
435
436
437

To quantify these differences, we have produced (Figure 12) difference maps between the combined mean maps of the two initialization setups (SAV-VS, SAV-VS-SPL and SAV-VS-SPL-PL). In Figure 12, maps initialized with ph are subtracted from the maps initialized with tv: this implies that negative values are those where ph maps have higher mean probability and vice-versa. In Figure 12, we observe that tv maps are characterized by higher hazard values in areas proximal (< 25 km to the W) with respect to Sangay location (with values never exceeding +5% in Figure 12), while ph maps imply higher probabilities on medial-distal areas (values up to 16% more than mfr maps in Figure 12).

438
439
440
441
442

Despite these differences, there is not a strong physical evidence for preferring one type of initialization over the other, because all the ESPs are significantly uncertain. However, we note that during the elicitation the experts provided a numerical value (in a scale from 1 to 10) expressing their own confidence in estimating the considered ESPs, and in this regard average plume height attained a slightly higher confidence with respect to total volume and duration (Bernard et al. 2024).

443
444
445

While we are not suggesting a preference for an initialization over the other, we stress that the choice made might have implications for vulnerability/risk assessment studies. For example, considering the SAV-VS-SPL maps, the probability (% - [5th-Mean-95th]) of having ≥ 1 kg/m² threshold (roughly corresponding to 1 mm of ash accumulation, sufficient for
airport closure) for the site of Guayaquil is, respectively, [0.6-1.3-2.3] for the tv initialization and [1.8-7.4-15.0] for the ph initialization (see Fig. 8 and Supporting Information 3). Similarly, the probability (%) - [5th-Mean-95th)] of having ≥ 10 kg/m² threshold (roughly corresponding to 10 mm of ash accumulation, sufficient for minor structural damages to edifices) for the site of Guamote is, respectively, [1.4-3.8-7.0] for the tv initialization and [4.1-13.3-25.4] for the ph initialization (see Fig. 8 and Supporting Information 3). In both examples the differences between the two model initializations is 3-5 times higher for the ph.

By looking at Figures 5-6 (and Fig. 19S from Supporting Information 2), we finally observe that, for SAV/VS initialized with tv and SPL/PL initialized with ph, our filtering approach rejected a larger number of ESPs. We have therefore produced a ‘mixed’ set of maps, in which we applied the same procedure of section 3.3 (equation [5]), but combining only the maps with the ph initialization for SAV and VS, and the maps with the tv initialization for SPL/PL. Results (maps and hazard curves) are displayed in Figures 26S-31S from Supporting Information 2.

## 6. Concluding remarks

This study provides a new tephra fallout hazard assessment for Sangay volcano through the production of probabilistic hazard maps and curves. These products have been developed using the coupled PLUME-MoM/HYSPLIT models with inputs from a stochastic sampling of atmospheric conditions and eruptive source parameters. The model itself has been initialized with either total volume and average plume height, producing two different types of maps. Four different scenarios (corresponding to eruptive pulses, Bernard et al. 2024) have been considered, and maps for each scenario have been produced for four different ground load thresholds (0.1, 1, 10 and 100 kg/m²). Hazard curves (representing exceedance probabilities) have also been produced in nine sites located both in the surroundings and at more distal locations from Sangay volcano. Maps for single eruptive scenarios have been combined according to the probability of having at least one such scenario over the next 100 years. Such probabilities (along with the eruptive source parameters for each scenario) and their uncertainties have been quantified through an expert elicitation session described in the companion paper by Bernard et al. (2024). Combined maps and curves take into account the uncertainty in probability of occurrences of different eruptive pulses in the next 100 years, and are therefore presented, for each combination setup, as mean maps and 5th/95th percentile maps, the latter two representing the uncertainty bounds.

Main results include:

* maps and curves for both single and combined scenarios indicate that the highest tephra fallout hazard is toward the W of Sangay volcano, consistently with observed volcanic activity and dominant wind directions. The town of Guamote (43 km W of Sangay) is the site, among the investigated ones, with the highest hazard, as it could be potentially affected by significant tephra accumulation (1 kg/m²) with probabilities up to 21-24% even for a moderate-scale eruption (Violent Strombolian);

* for combined scenarios, maps and curves were produced by integrating, respectively, i) the three historically observed eruptive pulses (Strong Ash Venting, Violent Strombolian and sub-Plinian), ii) the previous three and a larger Plinian event and iii) the two smallest observed pulses (Strong Ash Venting and Violent Strombolian. Average differences between combinations i) and ii) at the nine locations with hazard curves are between ~6% and ~9% for all the considered ground load thresholds and the two initialization setups. However, average differences up to 28%-40% could be observed for few sites and ground loads;

* differences between initialization setups have been quantified for the Mean maps, and indicate that maps with total volume initialization produce higher hazards (+5% maximum) for locations proximal to vent area (i.e. < 25 km to the W of Sangay) but lower hazards (-15% maximum) for medial-distal locations; possible reasons for this effect are the assumption of an erupted volume equal to the deposited volume, and the assumed independency of plume height and wind conditions.

All the outcomes summarized above clearly demonstrate the numerous sources of uncertainty that affect hazard maps and curves, related to the eruptive source parameters and atmospheric conditions (when producing maps for single scenarios), the type of model initialization (when adopting different source parameters) and the probability of occurrence of different scenarios (when combining maps from multiple scenarios). Nevertheless, a first comparison of the recent activity of Sangay is consistent with our hazard product. We believe that the set of maps and curves produced, obtained under different assumptions and with the associated uncertainty, can represent a valuable information for authorities and decision-makers aimed at mitigating the tephra fallout hazard at Sangay volcano.
498 References

499 Annen, C., and Wagner, J. J. (2003). The impact of volcanic eruptions during the 1990s. Natural Hazards Review, 44(4), 1694–175. https://doi.org/10.1061/(ASCE)1527-6988(2003)4:4(169).

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

504 Aspinall, W. P. (2006). Structured elicitation of expert judgment for probabilistic hazard and risk assessment in volcanic eruptions. Statistics in volcanology, 1, 15–30.

506 Aspinall, W. P., and Cooke, R. M. (2013). Quantifying scientific uncertainty from expert judgment elicitation. In: Rougier J, Sparks RSJ, Hill LJ (eds) Risk and uncertainty assessment for natural hazards. Cambridge University Press, New York, pp 64–99. https://doi.org/10.1017/CBO9781139047562.

509 Barker, S. J., Van Eaton, A. R., Mastin, L. G., Wilson, C. J., Thompson, M. A., Wilson, T. M., Davis, C., and Renwick, J. A. (2019). Modeling ash dispersal from future eruptions of Taupo supervolcano. Geochemistry, Geophysics, Geosystems, 20(7), 3375-3401. https://doi.org/10.1029/2018GC008152.

512 Barsotti, S., Andronico, D., Neri, A., Del Carlo, P., Baxter, P. J., Aspinall, W. P., and Hincks, T. (2010). Quantitative assessment of volcanic ash hazards for health and infrastructure at Mt. Etna (Italy) by numerical simulation. Journal of Volcanology and Geothermal Research, 192(1-2), 85-96. https://doi.org/10.1016/j.jvolgeores.2010.02.011.

515 Barsotti, S., Di Rienzo, D. I., Thordarson, T., Björnsson, B. B., & Karlsdóttir, S. (2018). Assessing impact to infrastructures due to tephra fallout from Öræfajökull volcano (Iceland) by using a scenario-based approach and a numerical model. Frontiers in Earth Science, 6, 196. https://doi.org/10.3389/feart.2018.00196.

518 Baxter, P. J., and Horwell, C. J. (2015). Impacts of eruptions on human health. In H. Sigurdsson, B. F. Houghton, S. R. McNutt, H. Rymer, and J. Stix (Eds.), The encyclopedia of volcanoes (pp. 1035–1047). Elsevier. https://doi.org/10.1016/B978-0-12-385938-9.00060-2.

521 Bebbington, M. S., Cronin, S. J., Chapman, I., and Turner, M. B. (2008). Quantifying volcanic ash fall hazard to electricity infrastructure. Journal of Volcanology and Geothermal Research, 177(4), 1055–1062. https://doi.org/10.1016/j.jvolgeores.2008.07.023.

524 Bernard, B., Samaniego, P., Mastin, L., Hernandez, S., Pino, G., Kibler, J., Encalada, M., Hidalgo, S., and Vizuete, N. (2022). Forecasting and communicating the dispersion and fallout of ash during volcanic eruptions: lessons from the September 20, 2020 eruptive pulse at Sangay volcano, Ecuador. Frontiers in Earth Science 10: https://doi.org/10.3389/feart.2022.912835.

528 Bernard, B., Tadini, A., Samaniego, P., Bevilacqua, A., Vasconez, F. J., Aravena, A., de' Michieli Vitturi, M., and Hidalgo, S. (2024). Developing hazard scenarios from monitoring data, historical chronicles and expert elicitation: case study of Sangay volcano, Ecuador. https://doi.org/10.21203/rs.3.rs-3764842/v1.

531 Bevilacqua, A., Isaia, R., Neri, A., Vitale, S., Aspinall, W. P., Bisson, M., Flandoli, F., Baxter, P. J., Bertagnini, A., Esposti Ongaro, T., Iannuzzi, E., Pistolesi, M., and Rosi, M. (2015). Quantifying volcanic hazard at Campi Flegrei caldera (Italy) with uncertainty assessment: I. Vent opening maps. Journal of Geophysical Research: Solid Earth, 120, 2309–2329. https://doi.org/10.1002/2014JB011775.

535 Bevilacqua, A., Flandoli, F., Neri, A., Isaia, R., and Vitale, S. (2016). Temporal models for the episodic volcanism of Campi Flegrei caldera (Italy) with uncertainty quantification. Journal of Geophysical Research: Solid Earth, 121, 7821–7845, https://doi.org/10.1002/2016JB013171.

538 Bevilacqua, A., Bursik, M., Patra, A., Pitman, E.B., Yang, Q., Sangani, R., Kobs- Nawotniak, S., (2018). Late Quaternary eruption record and probability of future volcanic eruptions in the Long Valley volcanic region (CA, USA). J. Geophys. Res. Solid Earth 123 (7), 5466–5494. https://doi.org/10.1029/2018JB015644.

541 Bevilacqua, A., Aravena, A., Neri, A., Gutiérrez, E., Escobar, D., Schliz, M., Aiuppa, A., and 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, 1639–1665. https://doi.org/10.5194/nhess-21-1639-2021.
545 Bevilacqua, A., Aravena, A., Aspinall, W., Costa, A., Mahony, S., Neri, A., Sparks, R., and Hill, B. (2022). Assessing minimum pyroclastic density current mass to impact critical infrastructures: example from Aso caldera (Japan).
546 Natural Hazards and Earth System Sciences, 22(10), 3329-3348. https://doi.org/10.5194/nhess-22-3329-2022.

548 Biass, S., and Bonadonna, C. (2011). A quantitative uncertainty assessment of eruptive parameters derived from tephra deposits: The example of two large eruptions of Cotopaxi volcano, Ecuador. Bulletin of Volcanology, 73(1), 73-90.
550 https://doi.org/10.1007/s00445-010-0404-5.

551 Biass, S., and Bonadonna, C. (2013). A fast GIS-based risk assessment for tephra fallout: The example of Cotopaxi volcano, Ecuador. Natural Hazards, 65(1), 477-495. https://doi.org/10.1007/s11069-012-0378-z.

553 Blong, R. J. (1996). Volcanic hazards risk assessment. In R. Scarpa, & R. I. Tilling (Eds.), Monitoring and mitigation of volcano hazards (pp. 675-698). Berlin: Springer. https://doi.org/10.1007/978-3-642-80087-0_20.

555 Bonadonna, C. (2006). Probabilistic modelling of tephra dispersion. In H. M. Mader, S. G. Coles, C. B. Connor, & L. J. Connor (Eds.), Statistics in Volcanology (pp. 243-259). London: Geological Society.

557 Bonasia, R., Capra, L., Costa, A., Macedonio, G., and Saucedo, R. (2011). Tephra fallout hazard assessment for a Plinian eruption scenario at Volcán de Colima (Mexico). Journal of Volcanology and Geothermal Research, 203(1-2), 12-22.
558 https://doi.org/10.1016/j.jvolgeores.2011.03.006.
559

560 Bonasia, R., Scaini, C., Capra, L., Nathenson, M., Siebe, C., Arana-Salinas, L., and Folch, A. (2014). Long-range hazard assessment of volcanic ash dispersal for a Plinian eruptive scenario at Popocatépetl volcano (Mexico): implications for civil aviation safety. Bulletin of Volcanology, 76, 1-16. https://doi.org/10.1007/s00445-013-0789-z.

563 Constantinescu, R., White, J. T., Connor, C. B., Hopulele-Gligor, A., Charbonnier, S., Thouret, J. C., Lindsay, J. M., and Bertin, D. (2022). Uncertainty quantification of eruption source parameters estimated from tephra fall deposits. Geophysical Research Letters, 49(6), e2021GL097425. https://doi.org/10.1029/2021GL097425.

566 Cooke, R. M. (1991). Experts in uncertainty: Opinion and subjective probability in science.

567 de' Michieli Vitturi, M., Engwell, S. L., Neri, A., and Barsotti, S. (2016). Uncertainty quantification and sensitivity analysis of volcanic columns models: Results from the integral model PLUME-MoM. Journal of Volcanology and Geothermal Research, 326, 77-91. https://doi.org/10.1016/j.jvolgeores.2016.03.014.

570 de' Michieli Vitturi, M., and Pardini, F. (2021). PLUME-MoM-TSM 1.0: A volcanic column and umbrella cloud spreading model. Geoscientific Model Development, 14(3), 1345-1377. https://doi.org/10.5194/gmd-14-1345-2021.

572 de' Michieli Vitturi, M., Bevilacqua, A., Tadini, A., and Neri, A. (2024). ELICIPY 1.0: a Python online tool for expert elicitation, SoftwareX, 25. https://doi.org/10.1016/j.softx.2024.101641.

574 Flandoli F., Giorgi, E., Aspinall, W. P., and Neri, A. (2011). Comparison of a new expert elicitation model with the classical model, equal weights and single experts, using a cross-validation technique. Reliability and Enginnering System Safety, 96:1292-1310. https://doi.org/10.1016/j.ress.2011.05.012.

577 Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R.J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J. N. (2020). The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological Society, 146(730), 1999-2049.
583 https://doi.org/10.1002/qj.3803.

584 Hidalgo, S., Vasconez, F. J., Battaglia, J., Bernard, B., Espín, P., Valade, S., Naranjo, M.-F., Campion, R., Salgado, J., Córdova, M., Almeida, M., Hernández, S., Pino, G., Gaunt, E., Bell, A., Mothes, P., and Ruiz, M. (2022). Sangay volcano (Ecuador): The opening of two new vents, a drumbeat seismic sequence and a new lava flow in late 2021. Volcanica, 5(2), Article 2. https://doi.org/10.30909/vol.05.02.295311.

588 Hurst, T., and Smith, W. (2004). A Monte Carlo methodology for modelling ashfall hazards. Journal of Volcanology and Geothermal Research, 138(3-4), 393-403. https://doi.org/10.1016/j.jvolgeores.2004.08.001.

590 Jenkins, S., McAneney, J., Magill, C., and Blong, R. (2012). Regional ash fall hazard II: Asia-Pacific modelling results and implications. Bulletin of volcanology, 74, 1713-1727. https://doi.org/10.1007/s00445-012-0628-7.
592 Jenkins, S. F., Wilson, T. M., Magill, C. R., Miller, V., Stewart, C., Marzocchi, W., and Boulton, M. (2015). Volcanic ash fall hazard and risk: Technical Background Paper for the UNISDR 2015 Global Assessment Report on Disaster Risk Reduction (Global Volcano Model and IAVCEI). www.preventionweb.net/english/hyogo/gar.

595 Ligot, N., Viera, W., Peñaherrera, D., Bernard, B., Bogaert, P., and Delmelle, P. (2024). A quantitative assessment of crop vulnerability to tephra hazard at Tungurahua volcano, Ecuador: Understanding the effect of volcanic and biological factors. Bulletin of Volcanology, 86(4), 28. https://doi.org/10.1007/s00445-024-01710-2.

598 Macedonio, G., and Costa, A. (2012). Brief Communication "Rain effect on the load of tephra deposits". Natural Hazards and Earth System Sciences, 12(4), 1229-1233. https://doi.org/10.5194/nhess-12-1229-2012.

600 Macedonio, G., Costa, A., Scollo, S., and Neri, A. (2016). Effects of eruption source parameter variation and meteorological dataset on tephra fallout hazard assessment: Example from Vesuvius (Italy). Journal of Applied Volcanology, 5(1). 1. https://doi.org/10.1186/s13617-016-0045-2.

603 Madankan, R., Pouget, S., Singla, P., Bursik, M. I., Dehn, J., Jones, M., Patra, A., Pavolonis, M., Pitman, E. B., Singh, T., and Webley, P. (2014). Computation of probabilistic hazard maps and source parameter estimation for volcanic ash transport and dispersion. Journal of Computational Physics, 271, 39-59. https://doi.org/10.1016/j.jcp.2013.11.032.

607 Magill, C., Wilson, T., and Okada, T. (2013). Observations of tephra fall impacts from the 2011 Shinmoedake eruption, Japan. Earth, Planets and Space, 65, 677-698. https://doi.org/10.5047/eps.2013.05.010.

609 Mastin, L. G., Guffanti, M., Servranckx, R., Webley, P., Barsotti, S., Dean, K., Durant, A., Ewert, J.W., Neri, A., Rose, W.I., Schneider, D., Siebert, L., Stunder, B., Swanson, G., Tupper, A., Volentik, A., and Waythomas, C. F. (2009). A multidisciplinary effort to assign realistic source parameters to models of volcanic ash-cloud transport and dispersion during eruptions. Journal of volcanology and Geothermal Research, 186(1-2), 10-21. https://doi.org/10.1016/j.jvolgeores.2009.01.008.

614 Mastin, L. G., Van Eaton, A. R., Schneider, D., and Denlinger, R. P. (2017). Ongoing efforts to make ash-cloud model forecasts more accurate. NATO Science and Technology Organization: Brussels, Belgium.

616 Michaud-Dubuy, A., Carazzo, G., and Kaminski, E. (2021). Volcanic hazard assessment for tephra fallout in Martinique. Journal of Applied Volcanology, 10, 1-20. https://doi.org/10.1186/s13617-021-00106-7.

618 Monzier, M., Robin, C., Samaniego, P., Hall, M. L., Cotten, J., Mothes, P., and Arnaud, N. (1999). Sangay volcano, Ecuador: structural development, present activity and petrology. Journal of Volcanology and Geothermal Research, 90(1-2), 49-79. https://doi.org/10.1016/S0377-0273(99)00021-9.

621 Moran-Zuloaga, D., Merchan-Merchan, W., Rodriguez-Caballero, E., Mulas, M., and Hernick, P. (2023). Long-range transport and microscopy analysis of Sangay volcanic ashes in Ecuador. Air Quality, Atmosphere & Health. https://doi.org/10.1007/s11869-023-01434-w.

624 Neri, A., Aspinall, W. P., Cioni, R., Bertagnini, A., Baxter, P. J., Zuccaro, G., Andronico, D., Barsotti, S., Cole, P. D., and Esposti Ongaro, T. (2008) Developing an event tree for probabilistic hazard and risk assessment at Vesuvius. Journal of Volcanology and Geothermal Research, 178:397-415. https://doi.org/10.1016/j.jvolgeores.2008.05.014.

627 Neri, A., Bevilacqua, A., Esposti Ongaro, T., Isaia, R., Aspinall, W. P., Bisson, M., Flandoli, F., Baxter, P. J., Bertagnini, A., and Iannuzzi, E. (2015). Quantifying volcanic hazard at Campi Flegrei caldera (Italy) with uncertainty assessment: II. Pyroclastic density current invasion maps, Journal of Geophysical Research: Solid Earth, 120, 2330-2349, https://doi.org/10.1002/2014JB011776.

631 Ordoñez, J., Vallejo, S. V., Bustillos, J., Hall, M., Andrade, D., Hidalgo, S., and Samaniego, P. (2011). Volcán Sangay, peligros volcánicos potenciales. IG-EPN/IRD. Escala 1:100.000

633 Pardini, F., Corradini, S., Costa, A., Esposti Ongaro, T., Merucci, L., Neri, A., Stelitano, D., and de' Michieli Vitturi, M. (2020). Ensemble-based data assimilation of volcanic ash clouds from satellite observations: application to the 24 December 2018 Mt. Etna explosive eruption. Atmosphere, 11(4), 359. https://doi.org/10.3390/atmos11040359.

636 Rose, W. I., and Durant, A. J. (2009). Fine ash content of explosive eruptions. Journal of Volcanology and Geothermal Research, 186(1-2), 32-39. https://doi.org/10.1016/j.jvolgeores.2009.01.010.
638 Sandri, L., Costa, A., Selva, J., Tonini, R., Macedonio, G., Folch, A., and Sulpizio, R. (2016). Beyond eruptive scenarios: Assessing tephra fallout hazard from Neapolitan volcanoes. Scientific Reports, 6(1), 1–13.
639 https://doi.org/10.1038/srep24271.
640

641 Stefanescu, E. R., Patra, A. K., Bursik, M. I., Madankan, R., Pouget, S., Jones, M., Singla, P., Singh, T., Pitman, E. B., Pavolonis, M., Morton, D., Webley, P., and Dehn, J. (2014). Temporal, probabilistic mapping of ash clouds using wind field stochastic variability and uncertain eruption source parameters: Example of the 14 April 2010 Eyjafjallajökull eruption. Journal of Advances in Modeling Earth Systems, 6(4), 1173–1184.
642
643
644
645 https://doi.org/10.1002/2014MS000332.

646 Stein, A. F., Draxler, R. R., Rolph, G. D., Stunder, B. J. B., Cohen, M. D., and Ngan, F. (2015). NOAA's HYSPLIT atmospheric transport and dispersion modeling system. Bulletin of the American Meteorological Society, 96(12), 2059–2077. https://doi.org/10.1175/BAMS-D-14-00110.1.

649 Tadini, A., Bisson, M., Neri, A., Cioni, R., Bevilacqua, A., and Aspinall, W. P. (2017a). Assessing future vent opening locations at the Somma-Vesuvio volcanic complex: 1. A new information geodatabase with uncertainty characterizations, Journal of Geophysical Research: Solid Earth, 122, 4336–4356.
650
651
652 https://doi.org/10.1002/2016JB013858.

653 Tadini, A., Bevilacqua, A., Neri, A., Cioni, R., Aspinall, W. P., Bisson, M., Isaia, R., Mazzarini, F., Valentine, G. A. V., Vitale, S., Baxter, P. J., Bertagnini, A., Cerminara, M., de' Michieli Vitturi, M., Di Roberto, A., Engwell, S. L., Esposti Ongaro, T., Flandoli, F., and Pistolesi, M. (2017b). Assessing future vent opening locations at the Somma-Vesuvio volcanic complex: 2. Probability maps of the caldera for a future Plinian/sub-Plinian event with uncertainty quantification. Journal of Geophysical Research: Solid Earth, 122:4357–4376.
654
655
656
657
658 https://doi.org/10.1002/2016JB013860.

659 Tadini, A., Roche, O., Samaniego, P., Guillin, A., Azzaoui, N., Gouhier, M., de' Michieli Vitturi, M., Pardini, F., Eychenne, J., Bernard, B., Hidalgo, S., and Le Pennec, J. L. (2020). Quantifying the uncertainty of a coupled plume and tephra dispersal model: PLUME-MOM/HYSPLIT simulations applied to Andean volcanoes. Journal of Geophysical Research: Solid Earth, 125(2). https://doi.org/10.1029/2019JB018390.

663 Tadini, A., Roche, O., Samaniego, P., Azzaoui, N., Bevilacqua, A., Guillin, A., Gouhier, M., Bernard, B., Aspinall, W., Hidalgo, S., Eychenne, J., de' Michieli Vitturi, M., Neri, A., Cioni, R., Pistolesi, M., Gaunt, E., Vallejo, S., Encalada, M., Yepes, H., Proaño, A., and Pique, M. (2021). Eruption type probability and eruption source parameters at Cotopaxi and Guagua Pichincha volcanoes (Ecuador) with uncertainty quantification. Bulletin of Volcanology, 83(35). https://doi.org/10.1007/s00445-021-01458-z.

668 Tadini, A., Azzaoui, N., Roche, O., Samaniego, P., Bernard, B., Bevilacqua, A., Hidalgo, S., Guillin, A., and Gouhier, M. (2022). Tephra fallout probabilistic hazard maps for Cotopaxi and Guagua Pichincha volcanoes (Ecuador) with uncertainty quantification. Journal of Geophysical Research: Solid Earth, 127, e2021JB022780.
669
670
671 https://doi.org/10.1029/2021JB022780.

672 Titos, M., Martínez Montesinos, B., Barsotti, S., Sandri, L., Folch, A., Mingari, L., Macedonio, G., and Costa, A. (2022). Long-term hazard assessment of explosive eruptions at Jan Mayen (Norway) and implications for air traffic in the North Atlantic. Natural Hazards and Earth System Sciences, 22(1), 139–163. https://doi.org/10.5194/nhess-22-139-2022.

676 Vasconez, F. J., Hidalgo, S., Battaglia, J., Hernandez, S., Bernard, B., Coppola, D., Valade, S., Ramón, P., Arellano, S., Liorzou, C., Almeida, M., Ortíz, M., Córdova, J., and Vásconez Müller, A. (2022). Linking ground-based data and satellite monitoring to understand the last two decades of eruptive activity at Sangay volcano, Ecuador. Bulletin of Volcanology 84:49. https://doi.org/10.1007/s00445-022-01560-w.

680 Vázquez, R., Bonasia, R., Folch, A., Arce, J. L., and Macías, J. L. (2019). Tephra fallout hazard assessment at Tacaná volcano (Mexico). Journal of South American Earth Sciences, 91, 253–259.
681
682 https://doi.org/10.1016/j.jsames.2019.02.013

683 Wilson, T. M., Stewart, C., Sword-Daniels, V., Leonard, G. S., Johnston, D. M., Cole, J. W., et al. (2012). Volcanic ash impacts on critical infrastructure. Physics and Chemistry of the Earth: Parts A/B/C, 45, 5–23.
684
685 https://doi.org/10.1016/j.pce.2011.06.006.
686 Woodhouse, M.J., Hogg, A.J., Phillips, J.C. and Rougier, J. C. (2015). Uncertainty analysis of a model of wind-blown

687 volcanic plumes. Bull Volcanol 77, 83. https://doi.org/10.1007/s00445-015-0959-2.

688 Yang, Q., Pitman, E. B., Spiller, E., Bursik, M., and Bevilacqua, A. (2020). Novel statistical emulator construction for

689 volcanic ash transport model Ash3d with physically motivated measures. Proc. R. Soc. A 476:20200161.

690 http://dx.doi.org/10.1098/rspa.2020.0161.

691 Yang, Q., and Jenkins, S. F. (2023). Two sources of uncertainty in estimating tephra volumes from isopachs: perspectives

692 and quantification. Bulletin of Volcanology, 85(8), 44. https://doi.org/10.1007/s00445-023-01652-1.

693

694

695

696

697

698

699

700

701

702

703

704

705

706

707

708

709

710

711

712

713

714

715

716

717

718

719

720

721

722
<table><caption>Table 1. Measured ground loads at selected locations from eruptions at Sangay volcano (period June 2020 - June 2024).</caption><thead><tr><th rowspan="2">Measured Ground load (kg/m²) for different eruption type (date)</th><th colspan="5">Town</th></tr><tr><th>Alausí</th><th>Guamote</th><th>Guaranda</th><th>Guayaquil</th><th>Riobamba</th></tr></thead><tbody><tr><td>SAV (8-11/6/20)</td><td>0.11</td><td>-</td><td>-</td><td>0.007</td><td>-</td></tr><tr><td>VS (20/9/20)</td><td>0.176-0.187</td><td>0.19</td><td>-</td><td>0.015-0.036</td><td>-</td></tr><tr><td>VS (5-6/3/21)</td><td>2.457</td><td>0.041</td><td>-</td><td>-</td><td>-</td></tr><tr><td>VS (11/3/21)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.105-0.710</td></tr><tr><td>VS (12/4/21)</td><td>-</td><td>-</td><td>0.024-0.046</td><td>0.008</td><td>-</td></tr><tr><td>SAV (12-13/8/22)</td><td>0.051</td><td>0.026</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SAV (4-5/11/22)</td><td>0.021</td><td>0.128</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SAV (18-21/4/23)</td><td>0.002</td><td>0.207</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SAV (21/5/23)</td><td>0.02</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SAV (7/1/24)</td><td>-</td><td>0.176</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SAV (27/1/24)</td><td>0.029</td><td>0.036</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SAV (12/6/24)</td><td>-</td><td>-</td><td>-</td><td>0.003</td><td>-</td></tr></tbody></table>

<table><caption>Table 2. Comparison between probability of exceeding different ground load (GL) thresholds (0.1 and 1 kg/m²) and the percentage of eruptive pulses which have produced ground loads larger than the same thresholds at five sites. The range of some exceeding probabilities takes into account the differences between the two initialization setups (tv and ph) for</caption><thead><tr><th>Town</th><th>Eruption Type</th><th>Initialization</th><th>Modeled probability of exceeding 0.1 kg/m²</th><th>Binomial probability of the observed eruptions with GL > 0.1 kg/m²</th><th>Modeled probability of exceeding 1 kg/m²</th><th>Binomial probability of the observed eruptions with GL > 1 kg/m²</th></tr></thead><tbody><tr><td rowspan="4">Alausí</td><td rowspan="2">SAV</td><td>tv</td><td>20%</td><td>34%</td><td>3%</td><td>79%</td></tr><tr><td>ph</td><td>32%</td><td>17%</td><td>9%</td><td>47%</td></tr><tr><td rowspan="2">VS</td><td>tv</td><td>44%</td><td>37%</td><td>16%</td><td>38%</td></tr><tr><td>ph</td><td>46%</td><td>37%</td><td>18%</td><td>40%</td></tr><tr><td rowspan="4">Guamote</td><td rowspan="2">SAV</td><td>tv</td><td>25%</td><td>20%</td><td>5%</td><td>69%</td></tr><tr><td>ph</td><td>38%</td><td>28%</td><td>13%</td><td>34%</td></tr><tr><td rowspan="2">VS</td><td>tv</td><td>51%</td><td>24%</td><td>24%</td><td>33%</td></tr><tr><td>ph</td><td>48%</td><td>27%</td><td>21%</td><td>39%</td></tr><tr><td rowspan="4">Guaranda</td><td rowspan="2">SAV</td><td>tv</td><td>4%</td><td>74%</td><td>0%</td><td>98%</td></tr><tr><td>ph</td><td>7%</td><td>56%</td><td>0%</td><td>98%</td></tr><tr><td rowspan="2">VS</td><td>tv</td><td>14%</td><td>54%</td><td>2%</td><td>92%</td></tr><tr><td>ph</td><td>14%</td><td>55%</td><td>3%</td><td>88%</td></tr><tr><td rowspan="4">Guayaquil</td><td rowspan="2">SAV</td><td>tv</td><td>0%</td><td>97%</td><td>0%</td><td>99%</td></tr><tr><td>ph</td><td>3%</td><td>81%</td><td>0%</td><td>99%</td></tr><tr><td rowspan="2">VS</td><td>tv</td><td>6%</td><td>77%</td><td>1%</td><td>98%</td></tr><tr><td>ph</td><td>6%</td><td>78%</td><td>0%</td><td>99%</td></tr><tr><td rowspan="4">Riobamba</td><td rowspan="2">SAV</td><td>tv</td><td>7%</td><td>58%</td><td>1%</td><td>94%</td></tr><tr><td>ph</td><td>11%</td><td>41%</td><td>2%</td><td>88%</td></tr><tr><td rowspan="2">VS</td><td>tv</td><td>21%</td><td>41%</td><td>7%</td><td>74%</td></tr><tr><td>ph</td><td>21%</td><td>41%</td><td>8%</td><td>71%</td></tr></tbody></table>
728 the single scenarios (SAV and VS - see Supporting Information 3). Cells with '-' indicate that for this site/eruption type,

729 no eruption had a GL ≥ the considered threshold.

730

731 **Figures**

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

732

733 **Figure 1.** a) Location of Sangay volcano (red triangle) in Ecuador (the black dotted line marks national boundaries). The

734 blue square represents the location of the capital city (Quito) and the black dots indicate some towns/cities potentially

735 affected by Sangay volcano tephra fallout; b) picture of Sangay volcano on 27/12/2021 (taken from 6 km WSW from the

736 vent).
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Probability Density Functions (PDF) of the uncertainty on the probability (for each scenario - SAV, VS, SPL, PL) of having at least one eruption in the next 100 years (Bernard et al. 2024).</figcaption></figure>

737

738

739
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Probability density functions of Mass Flow Rate (expressed as log₁₀) and plume height (in km above crater) of the SAV, VS and SPL eruptive scenarios. The figure refers to model initialization with total volume (tv) and the filtering done considering the plume height (ph) (left panels) and model initialization with ph and the filtering done considering the tv (obtained from mass flow rate from the model, converted into total mass considering eruption duration and finally converted into total volume; right panels). In each panel, the red dashed line indicates the original density distribution obtained from the 10000 samples produced from the results of the elicitation while the blue dashed line indicates the filtered density distribution obtained from the samples that passed the control (by running PLUME-MoM-TSM only). The 1200 simulations were taken from this latter distribution. Details on the filtering approach are provided in section 3.2.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Wind rose diagrams at the Sangay volcano location (2.00°S; 78.32°W), constructed from the ERA5 reanalysis dataset covering the years from 2012 to 2022. Wind speed and direction were evaluated for different height intervals above the crater: surface-2 km, 2-5 km, 5-10 km, 10-20 km, 20-35 km, and 35-50 km. Note that the diagrams show the direction toward which the wind blows.</figcaption></figure>

750

751

752

753

754
<figure><img src="image_6.png" /></figure>

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

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

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

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

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /></figure>

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

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

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

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

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

<caption>**Figure 6.** Hazard curves for 9 sites (see Figure 1) for the SAV, VS and SPL scenarios. Curves in the left columns are produced by initializing PLUME-MoM-TSM with total volume (tv), whereas those in the right column by using Plume Height (ph). Interpolated values at 0.1, 1, 10 and 100 kg/m² are reported in Supporting Information 3.</caption>
<figure><img src="image_18.png" /><figcaption>**Figure 7.** Sangay tephra fallout hazard maps for the combined SAV, VS and SPL scenarios for 1 kg/m² threshold corresponding to the 5th percentile, the mean and the 95th percentile. Maps in the left column are produced initializing PLUME-MoM-TSM with total volume (tv), whereas those in the right column using plume height (ph). Contour lines, where not indicated, follow the progression 1%, 5%, 10%, 50%, 80%.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 8.** Hazard curves for 3 sites (Guamote, Guayaquil and Macas - see Figure 1) for the combined SAV, VS and SPL scenarios corresponding to the 5th percentile, the mean and the 95th percentile. Curves in the left column are produced by initializing PLUME-MoM-TSM with total volume (tv), whereas those in the right column by using plume height (ph). Interpolated values at 0.1, 1, 10 and 100 kg/m² are reported in Supporting Information 3.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 9.** Upper panels: Sangay tephra fallout hazard maps for the PL scenario for 1 kg/m² threshold. Map in the left column is produced by initializing PLUME-MoM with total volume (tv), whereas that in the right column using Plume Height (ph). Contour lines, where not indicated, follow the progression 1%, 5%, 10%, 50%, 80%. Lower panels: hazard curves for 9 sites (see Figure 1) for the PL scenario. Curves in the left column are produced by initializing PLUME-MoM-TSM with total volume (tv), whereas those in the right column by using Plume Height (ph). Interpolated values at 0.1, 1, 10 and 100 kg/m² are reported in Supporting Information 3.</figcaption></figure>
<figure><img src="image_21.png" /></figure>

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

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

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

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

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

Figure 10. Sangay tephra fallout hazard maps for the combined SAV, VS, SPL and PL scenarios for 1 kg/m² threshold

corresponding to the 5th percentile, the mean and the 95th percentile. Maps in the left column are produced initializing

PLUME-MoM-TSM with total volume (tv), whereas those in the right column using Plume Height (ph). Contour lines,

where not indicated, follow the progression 1%, 5%, 10%, 50%, 80%.
<figure><img src="image_27.png" /><figcaption>**Figure 11.** Hazard curves for 3 sites (Guamote, Guayaquil and Macas - see Figure 1) for the combined SAV, VS, SPL and PL scenarios corresponding to the 5th percentile, the mean and the 95th percentile. Curves in the left column are produced by initializing PLUME-MoM-TSM with total volume (tv), whereas those in the right column by using Plume Height (ph). Interpolated values at 0.1, 1, 10 and 100 kg/m² are reported in Supporting Information 3.</figcaption></figure>
<figure><img src="image_28.png" /></figure>

Difference map (SAV,VS,SPL,PL)

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

Difference map (SAV,VS,SPL,PL)

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