# Quantifying the Uncertainty of a Coupled Plume and Tephra Dispersal Model: PLUME-MOM/HYSPLIT Simulations Applied to Andean Volcanoes

M. Alessandro Tadini, Olivier Roche, P. Samaniego, A. Guillin, N. Azzaoui, M. Gouhier, M. De'Michieli, F. Pardini, Julia Eychenne, B. Bernard, et al.

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

## To cite this version:

M. Alessandro Tadini, Olivier Roche, P. Samaniego, A. Guillin, N. Azzaoui, et al., 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, 2020, 125 (2), 10.1029/2019JB018390. [hal-02453781v1](https://HAL Id: hal-02453781)

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>
1

2

3

4

5

6

7

8

# Quantifying the uncertainty of a coupled plume and tephra dispersal model: PLUME-MOM/HYSPLIT simulations applied to Andean volcanoes

A. Tadini¹, O. Roche¹, P. Samaniego¹,², A. Guillin³, N. Azzaoui³, M. Gouhier¹, M. de' Michieli

Vitturi⁴, F. Pardini⁴, J. Eychenne¹, B. Bernard², S. Hidalgo², J. L. Le Pennec¹,⁵

¹Laboratoire Magmas et Volcans, Université Clermont Auvergne, CNRS, IRD, OPGC, F-63000 Clermont-Ferrand, France.

²Escuela Politécnica Nacional, Instituto Geofísico, Ladrón de Guevara E11-253 y Andalucía, Quito, Ecuador

³Laboratoire de Mathématiques Blaise Pascal, Université Clermont Auvergne, CNRS, IRD, OPGC, F-63000

Clermont-Ferrand, France.

⁴Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Pisa, Via Cesare Battisti 53, 56126, Pisa, Italy

⁵Institut de Recherche pour le Développement, Alemania N32-188 y Guayanas, Quito, Ecuador

Corresponding author: Alessandro Tadini (Alessandro.TADINI@uca.fr)

16

17

## Keypoints

18

19

* We present an uncertainty quantification for a coupled version of a plume model (PLUME-MoM) and a tephra dispersal model (HYSPLIT)

20

21

* The model has been tested against field data of 4 eruptions from Andean volcanoes (in Ecuador and Chile) of different magnitudes/styles

22

23

* The main conclusion of the uncertainty quantification is that the model is best suited for hazard studies of higher magnitude eruptions

24

25

26

27

28

29

30

31

32

33
34 **Abstract**

Numerical modelling of tephra dispersal and deposition is essential for evaluation of volcanic hazards. Many models consider reasonable physical approximations in order to reduce computational times, but this may introduce a certain degree of uncertainty in the simulation outputs. The important step of uncertainty quantification is dealt in this paper with respect to a coupled version of a plume model (PLUME-MoM) and a tephra dispersal model (HYSPLIT). The performances of this model are evaluated through simulations of four past eruptions of different magnitudes and styles from three Andean volcanoes, and the uncertainty is quantified by evaluating the differences between modeled and observed data of plume height (at different time steps above the vent) as well as mass loading and grain size at given stratigraphic sections. Different meteorological datasets were also tested and had a sensible influence on the model outputs. Other results highlight that the model tends to underestimate plume heights while overestimating mass loading values, especially for higher magnitude eruptions. Moreover, the advective part of HYSPLIT seems to work more efficiently than the diffusive part. Finally, though the coupled PLUME-MoM/HYSPLIT model generally is less efficient in reproducing deposit grain sizes, we propose it may be used for hazard maps production for higher magnitude eruptions (sub-Plinian or Plinian) for what concern mass loading.

52

53 **Index Terms and Keywords**

4314 Mathematical and computer modeling, 3275 Uncertainty quantification, 8428 Explosive volcanism, 8488 Volcanic hazards and risks

Tephra fall, tephra dispersal, numerical modelling, uncertainty quantification, Andean volcanoes

58

59 **1. Introduction**

Volcanic tephra dispersal and deposition represent a threat for many human activities since tephra may have a huge impact on aviation and can also damage edifices, infrastructures and vegetation when it accumulates on the ground, even in relatively small quantities. For this reason, numerical models have been developed over the past decades for describing both tephra rise into the eruptive column (plume models - PMs) or its transport by wind advection [tephra transport and dispersal models - TTDM; Folch, 2012]. Since describing in great detail the physics of such phenomena requires complex 3-D multiphase models, it is useful for operational purposes (e.g. volcanic ash tracking in real time or hazard maps production) to rely on simplified models, which introduce reasonable physical assumptions. In doing so, though computational times might be reduced, approximations and uncertainties are introduced in the final results of the simulations. Uncertainties need to be therefore quantified in order to facilitate decision makers in taking both real-time and long-term informed decisions. With respect to numerical models, uncertainty quantification in literature has been done: i) for PMs, by comparing modelled and observed values of maximum plume height (or level of neutral buoyancy) and/or of the mass flow rate (in kg/s), as for instance in Folch et al. [2016] or Costa et al. [2016]; ii) for TTDMs, by comparing modelled and observed ground deposit measurements (mass loadings in kg/m²) and/or ash cloud measurements
(concentrations in the atmosphere in kg/m³) [e.g., Scollo et al., 2008; Costa et al., 2009; Bonasia et al., 2010; Folch, 2012].

The aim of the present study is therefore twofold. Firstly, we present a coupled version of two different models: i) a renewed version of PLUME-MoM, a simplified 1-D plume model developed by de'Michieli Vitturi et al. [2015], and ii) the HYSPLIT model [Stein et al., 2015], a Lagrangian TTDM developed by the National Oceanic and Atmospheric Administration (NOAA) and currently used by several Volcanic Ash Advisory Centers (VAACs) to track and forecast volcanic clouds. Secondly, we provide a quantification of the uncertainty of the coupled version of these two models by testing simulations results with data of four different recent eruptions of three Andean volcanoes (Fig. 1). These eruptions were produced by Cotopaxi [2015 eruption, Bernard et al., 2016a] and Tungurahua [2006 eruption, Eychenne et al., 2012; 2013 eruption, Parra et al., 2016] volcanoes in Ecuador, and Puyehue-Cordón Caulle volcanic complex [2011 eruption, Pistolesi et al., 2015] in Chile. With this new coupled model the volcanic particles transport is simulated throughout the whole process that is within the eruptive column and through atmospheric dispersion. Furthermore, the uncertainty quantification represents an important aspect regarding hazard maps production.

In this article, after describing the eruptions chosen for the uncertainty quantification (section 2.1), we present the PLUME-MoM and HYSPLIT models as well as the coupling of these two models (section 2.2.1). Then we present the input parameters used for the simulations (Section 2.2.2) and we describe the strategy adopted for the quantification of the uncertainty of the coupled model (Section 2.3). Results presented in Section 3 serve as a basis for the discussion in Section 4 about the uncertainties related to the input parameters and the numerical models and about also the effectiveness of these models when used for producing tephra fallout hazard maps.

102

## 2. Background

### 2.1 Eruptions selected

The four eruptions chosen for testing our simulations cover a wide range of eruptive styles (sub-Plinian, violent strombolian, vulcanian, hydrovolcanic to long-lasting ash emission), durations (from few hours up to more than 3 months) and magma compositions (andesitic to rhyolitic/rhyodacitic). The criteria for selecting these eruptions were i) the location of the volcanoes in the same geodynamic context, ii) the existence of both detailed chronologies and meteorological data for the eruptions, and iii) the availability of reasonably well constrained input parameters for the models.

#### 2.1.1 Cotopaxi 2015

The 2015 eruption of Cotopaxi (C15 – Fig. 1a) started with hydromagmatic explosions on August 14th 2015, which produced a 9-10 km-high eruptive column above the crater and moderate ash fallout to the NW of the volcano. Then, it was followed by three and a half months of moderate to low ash emissions with plumes reaching on average 2 km above the crater and directed mostly to the west [Bernard et al., 2016a; Gaunt et al., 2016].

The magmatic character of the eruption increased through time as was shown by microtextural analysis [Gaunt et al., 2016] and ash/gas geochemistry [Hidalgo et al., 2018]. Through frequent sampling missions, the ash emission rate was calculated and correlated with the eruptive tremors, and it decreased during three emission phases following the conduit opening [Bernard et al., 2016a].

The fallout deposit was characterized by a very fine-grained ash with mostly blocky fragments and few vesicular scoria [Gaunt et al., 2016]. The hydrothermal components were dominant at the onset of the eruption but rapidly faded and were replaced by juvenile material
126
127

[Gaunt et al., 2016]. In total, this eruption emitted ~1.2x10⁹ kg of ash and was characterized as a VEI 1-2 [Bernard et al., 2016a].

128

## 2.1.2 Tungurahua 2013

129
130
131
132
133
134
135
136
137

According to *Hidalgo et al.* [2015], the eruptive phase XI (T13) at Tungurahua volcano (Fig. 1a) started on July 14th 2013 and lasted 23 days. A vulcanian onset, interpreted as the opening of a plugged conduit, was followed by a paroxysm which created a ~14 km-high eruptive column [*Parra et al.*, 2016]. The ash cloud created during this eruption was divided into a high cloud (~8-9 km above the crater) moving north and an intermediate cloud (~5 km above the crater) moving west and that produced most of the ash fallout [*Parra et al.*, 2016]. The eruption intensity dropped after this paroxysm but ash emission continued with a secondary increase between July 20th and 24th. Finally the eruption stopped at the beginning of August.

138
139
140

In total, this eruption emitted ~6.7x10⁸ kg of fallout deposits (~2.9x10⁸ kg for the first day) and ~5x10⁹ kg of pyroclastic flow deposits (mostly during the first day) [*García Moreno, 2016; Parra et al., 2016*].

141
142
143
144
145
146

*Parra et al.* [2016] performed numerical simulations of the vulcanian onset of this eruption, which occurred on July 14th 2013, using the coupled WRF-FALL3D models [*Michalakes et al., 2001; Folch et al., 2009*]. By comparing the mass loading between the modeled values and the observed ones at four sampling sites, the above-mentioned authors derived a set of Eruptive Source Parameters (ESPs) useful for operational purposes in case of vulcanian eruptions at Tungurahua volcano.

147

## 2.1.3 Tungurahua 2006

148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167

At Tungurahua volcano (Fig. 1a), a paroxysmal eruption (T06) occurred on August 16th 2006, which was accompanied by regional tephra fallout and many scoria flows and surges that devastated the western half of the edifice [*Douillet et al., 2013; Hall et al., 2013*]. This eruption was characterized by vigorous lava jetting and fountaining, a vent-derived eruption column reaching 16–18 km above the vent [*Steffke et al., 2010; Eychenne et al., 2012*], numerous Pyroclastic Density Currents (PDCs) descending the southern, western and northern flanks of the volcano [*Kelfoun et al., 2009; Bernard et al., 2014*], and a massive blocky lava flow emplacing on the western flank while the explosive activity waned [*Samaniego et al., 2011; Bernard et al., 2016b*]. At the climax of the eruptive event, after 3 hours of intense PDC formation, the vent-derived ash plume developed into a sub-vertical and sustained column for 50 to 60 minutes [*Hall et al., 2013*]. The plume spread over the Inter-Andean Valley, west of the volcano, and reached the Pacific Ocean, leading to substantial lapilli and ash fallout on the nearby communities and cities (e.g., Riobamba and Ambato) located to the West. The intense PDC activity generated ash-rich, 10 km-high co-PDC plumes that spread over the same areas and deposited fine ash (<90 µm) [*Eychenne et al., 2012; Bernard et al., 2016b*].

168

## 2.1.4 Puyehue-Cordón Caulle 2011

169
170

According to *Collini et al.* [2013], the Puyehue-Cordón Caulle 2011 eruption (PCC11 - Fig. 1b) started on June 4th at 14:45 LT (18:45 UTC) with the opening of a new vent 7 km
171 NNW from the main crater of the Puyehue-Cordón Caulle complex (“We Pillán” vent – Fig. 1b). The eruptive period, which involved mainly magma of rhyolitic-rhyodacitic composition [Bonadonna et al., 2015a], lasted up to June 2012 [Jay et al., 2014] and comprised both explosive and effusive activity [Tuffen et al., 2013]. The main explosive phase, which dispersed most of the tephra toward E and SE, lasted approximately 17-27 hours [Jay et al., 2014; Bonadonna et al., 2015b]. During the first three days of the eruption, the column rose approximately between 9 and 12 km above vent, then between 4 and 9 km during the following week, and less than 6 km after June 14th [Bonadonna et al., 2015a; Biondi et al., 2017].

180 During the eruption, the mass eruption rate (MER) fluctuated between 2.8x10⁷ (during the first days) and less than 5x10⁵ kg/s after June 7th [Bonadonna et al., 2015b]. Pistolesi et al. [2015] subdivided the stratigraphic record in thirteen tephra layers: among them, the first unit (Unit I, layers A-F) represented the tephra deposited between June 4th-5th. Unit I had a total erupted mass of 4.5±1.0×10¹¹ kg and was sub-Plinian with a VEI of 4 [Bonadonna et al., 2015b]. Bonadonna et al. [2015a] calculated the total grain size distribution (TGSD) of Unit I in the range -4φ/11φ, using different datasets and methods. The results indicated a bimodal distribution with the two sub-populations (with modes at -2φ and 7φ) separated by the 3φ grain size [Bonadonna et al., 2015a].

189 Collini et al. [2013] performed numerical modellings of this eruption between June 4th to June 20th using the above-mentioned WRF-FALL3D code. The authors compared both the column mass load (in ton/km²) and ground deposit measurements between modeled and observed values. With respect to deposit thickness measurements, they compared deposit thicknesses at 37 locations, resulting in a best-fit line on a computed versus observed graphs. The PCC11 eruption was furthermore modeled by Marti et al. [2017], who simulated the eruption from June 4th up to Jun 21st using the NMMB-MONARCH-ASH model and compared the same parameters as in Collini et al. [2013]. For the ground measurements, they provided comparisons between the simulated and observed isopach maps for both the Unit I and other eruptive units cited in Pistolesi et al. [2015], finding a good agreement between modeled and observed data.

## 2.2 Numerical modeling

### 2.2.1 Models used and coupling of the codes

202 For this work, the integral plume model PLUME-MoM has been coupled with HYSPLIT, one of the most extensively used atmospheric transport and dispersion models in the atmospheric sciences community.

205 Following the approach adopted in Bursik [2001], PLUME-MoM solves the equations for the conservation of mass, momentum, energy, and the variation of heat capacity and mixture gas constant. The model accounts for particle loss during the plume rise and for radial and crosswind air entrainment parameterized using two entrainment coefficients. In contrast to previous works, in which the pyroclasts are partitioned into a finite number of bins in the Krumbein scale, PLUME-MoM adopts the method of moments to describe a continuous size distribution of one or more group of particles (i.e. juveniles, lithics...). An uncertainty quantification and a sensitivity analysis of the PLUME-MoM model were done by de'Michieli Vitturi et al. [2016] by analyzing the distribution of plume heights obtained when varying a series of input parameters (i.e. air radial/wind entrainment, exit velocity, exit temperature, water fraction and wind intensity). The above-mentioned authors showed that plume height distribution was the widest when the parameters varied were the exit velocity, exit temperature, water fraction and wind intensity. With respect to the sensitivity, de'Michieli Vitturi et al. [2016] showed that initial water fraction had the strongest influence on plume
219 height determination (i.e. the plume height decreased by a factor of ~1.54 when increasing water content from 1 to 5 wt%).

221 HYSPLIT belongs to the family of Lagrangian Volcanic ash transport and dispersion models, which have been used operationally since the mid 1990's by the International Civil Aviation Organization (ICAO) to provide ash forecast guidance. The model solves the Lagrangian equations of motion for the horizontal transport of pollutants (i.e. particles), while vertical motion depends on the pollutant terminal fall velocity. The dispersion of a pollutant may be described using three main types of configuration, “3D particle” “puff” or hybrid “particle/puff”. Particularly, in the “puff” configuration, pollutants are described by packets of ash particles (“puffs”) having a horizontal Gaussian distribution of mass described by a standard deviation σ. The puffs expand with atmospheric turbulence until they exceed the size of the meteorological grid cell (either horizontally or vertically) and then split into several new puffs, each with their respective pollutant mass. In this work, the hybrid “particle/puff” configuration has been used, in which the horizontal packets of particles have a “puff” distribution, while in the vertical they move like 3D particles. This approach allows to use a limited number of puffs to properly capture both the horizontal dispersion and the vertical wind shears. Webley et al. [2009] have evaluated the sensitivity of the model with respect to the concentration of ash in the volcanic cloud when two parameters, TGSD and the vertical distribution of ash, were varied. The sensitivity analysis was done with respect to a test case eruption (Crater Peak/Mt. Spurr, Alaska, USA, 1992). They showed that three different TGSDs had little effect on the modeled ash cloud, while a uniform concentration of ash throughout the vertical eruptive column provided results more similar to satellite measurements. For this work, some modifications have been implemented in HYSPLIT and are described in Text S1 from the Supporting Information.

243 In the present study we coupled the PLUME-MoM and HYSPLIT models with an ad-hoc Python script, which computes for each grain size, from the output of the plume model, the mass rates released from the edges of the plume at intervals of fixed height, and the mass flow that reaches the neutral buoyancy level. Then, the script assembles an input file where the source locations for HYSPLIT are defined. In addition, it is employed a utility from the HYSPLIT package to extract the wind profile at the vent, in order to provide this information to the plume model. This coupled model was used for all the studied eruptions, while for some specific cases (i.e. the simulations for the PCC11 eruption) we also implemented a best-fitting inverse version of this coupling, which was based on the approach first described by Connor and Connor [2006] and applied, among others, by Bonasia et al. [2010] and Costa et al. [2009]. The parameters for which the inversion was performed and their range of variation were identified first. We considered the mass flow rate (in kg/s), the initial water mass fraction (in wt%) and the particle shape factor [Wilson and Huang, 1979; Riley et al., 2003]. We chose these parameters because their uncertainty was higher and/or the models were more sensitive to small variations of them. The procedure was aimed at minimizing the T² function

$$T^2 = \sum_{i=1}^{N} w_i [ML_{o,i} - ML_{m,i}]^2$$

258 where the sum is extended over N stratigraphic sections used in the inversion, wᵢ are weighting factors (in our case all are equal to 1), ML₀ᵢ denotes the observed mass load (in kg/m²) and MLₘᵢ are the values predicted by the model (in kg/m²). The values of T² will be then compared to the standard Chi-2 distribution of N-p degrees of freedom, with p=3 the number of free parameters.

263 264 2.2.2 Modelling features and input parameters
265

We tested four different types of meteorological data (GDAS, NCEP/NCAR, ERA-Interim, ERA-Interim refined using WRF/ARW; see Text S1 from the Supporting Information for details) with various spatial and temporal resolutions (see Table S1 in Supporting Information), which correspond to the most widely used meteo data for studies similar to ours.

266

267

268

269

270

All the HYSPLIT simulations were done using a 0.05° (~5 km) computational grid. After the end of each emission time (i.e. the actual duration of the eruption), a further amount of 12 hours was added to the simulation in order to allow finer particles to settle down. Simulations were performed in a forward way for all the four eruptions. However, a best-fitting inverse procedure (see Section 2.2.1) was performed for the PCC11 eruption because the uncertainty in the tephra fallout total mass estimation was the highest among the four chosen eruptions. A total of 600 inversions were performed, corresponding to 200 inversions for each of the three meteo data employed for a given eruption (GDAS, NCEP/NCAR and ERA-Interim).

271

272

273

274

275

276

277

278

279

Eruption source parameters (ESPs) were estimated from earlier works for the four eruptions and some of them are reported in Table 1 (the detailed list of parameters for each eruption is available in Table S2 in Supporting information). More specifically: a) the computational grid dimension (i.e. the total span of the computational domain in degrees with respect to the vent location) was defined in order to contain all or the vast majority (>95%) of the erupted mass and to reduce as much as possible the computational time; b) the initial water content was assumed as that of typical mean values for andesitic (for C15, T13 and T06) or rhyolitic (for PCC11) magmas, following Andújar et al. [2017] and Martel et al. [2018] respectively. For the inverse simulations of PCC11, the initial water content at each iteration was sampled between 6% and 8% [Martel et al., 2018]; c) Particles exit velocities from the vent were assigned two different values [following de'Michieli Vitturi et al., 2015] corresponding to a “weak plume” case (C15 and T13) or to a “strong plume” case (T06 and PCC11); d) The heat capacity of volcanic particles was assumed with a fixed value of 1600 J/kgK following Folch et al. [2016]; e) The particles shape factor was assumed with two different values for andesitic magmas (C15, T13 and T06) and for rhyolitic ones (PCC11) following the results of Riley et al. [2003]. For the inverse simulations of the PCC11 eruption, the particle shape factor values at each iteration were sampled between 0.6 and 0.8 [Riley et al., 2003]; h) the particle density values were assumed to vary linearly between two values ($\rho_1$ and $\rho_2$) specific of two grain sizes ($\phi_1$ and $\phi_2$) according to Bonadonna and Phillips [2003]. Values of $\rho_1$, $\rho_2$, $\phi_1$, and $\phi_2$ were taken from Eychenne and Le Pennec [2012] (C15/T13/T06) and Pistolesi et al. [2015] (PCC11). For each eruption, all the other most relevant features of input parameters are described below.

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

For the Cotopaxi C15 eruption, the simulations covered the whole eruption duration (14/08/2015 - 30/11/2015) for a total of 108 days and 17 hours. Plume heights values were obtained from Bernard et al. [2016a]. With respect to the TGSD calculated in Gaunt et al. [2016] we also used several unpublished data (see Table S1 from the Supporting Information). More specifically, a total of 33 samples representative of different times during the eruption and from 4 stratigraphic sections were employed. The TGSD was derived from a weighted mean (with respect to different mass loading values) of single grain size measurements. MER values used for the simulations were recalculated from Bernard et al. [2016a] to obtain hourly values (see Table S2 from the Supporting Information).

302

303

304

305

306

307

308

309

310

For the Tungurahua T13 eruption, the simulations also covered the whole eruption duration (14/07/2013 – 30/07/2013) for a total of 16 days and 12 hours. We considered observed plume height measurements from two sources: the ones by the Washington VAAC using satellite measurements, and those from observations made by the Tungurahua Volcano

311

312

313
314

Observatory (OVT). Similarly to the C15 eruption, the TGSD was obtained from a weighted mean (with respect to different mass loading values) of single grain size measurements. Hourly values of MER were obtained from unpublished data of the total mass deposited at the Choglontus sampling site at different intervals (Table S2 from the Supporting Information).

315

316

317

318

319

320

321

322

323

324

325

326

For the Tungurahua T06 eruption, the simulations covered 4 hours corresponding to the climatic phases I and II described in Hall et al. [2013]. Plume heights were derived from Steffke et al. [2010]. An average value of the MER was initially derived from the total mass deposited over this period (see Text S1 from the Supporting Information); successively, hourly values of MER were determined after an iterative procedure aimed at obtaining modeled output values of plume heights as close as possible to observed data. This iteration was done separately for each meteo data. The TGSD was recalculated from that of Eychenne et al. [2012] by removing the mass contribution of the co-PDC part (see Text S1 from the Supporting Information).

327

328

329

330

331

332

333

Finally, for the Puyehue-Cordón Caulle PCC11 eruption, the simulations covered the initial part of the eruption corresponding to the emplacement of Unit I [Pistolesi et al., 2015] for a total of 24 hours. Daily average plume heights and MERs from Bonadonna et al. [2015b] were employed along with a TGSD calculation from Bonadonna et al. [2015a]. For the inverse simulations, the MER was sampled between two values ($10^{6.75}$ and $10^{6.95}$ kg/s), which gave the minimum and maximum total mass values provided by Bonadonna et al. [2015b] and reported also in Table S2 (Supporting Information).

334

335

## 2.3 Uncertainty quantification procedure

336

337

We quantified the uncertainty of the coupled numerical model by comparing modeled and observed values of key parameters of both the PM and the TTDM.

338

339

340

341

With respect to the PM, we compared the plume height (in meters above vent) observed against the corresponding value at the same time (or at the closest measurement available) given by the model. In this case it is important to remember that plume height in PLUME-MoM is obtained as output value using a fixed MER.

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

For the TTDM, we compared ground deposit measurements and we adopted a specific approach in order to properly address uncertainty quantification. The results of the simulations were used to compare, at each stratigraphic section, observed and modeled values of mass loading and grain size, the latter one characterized by $Md\phi$ and $\sigma\phi$ [Folk and Ward, 1957]. For mass loading we use hereafter the notation "$\Delta$ mass loading", which corresponds to the difference between the computed and the observed values of mass loading (in kg/m$^2$). In the corresponding graphs (Figs. 3b, 4b, 5b and 6b) $\Delta$ mass loading values (for each simulation) and observed mass loadings are reported for each section. A complete list of the stratigraphic sections employed is available in Table S3 from the Supporting information. We considered also the direction of the main elongation axis of the deposit by comparing isomass maps constructed from field data and those given by the model. With respect to mass loading values, additional parameters were also calculated to quantify the uncertainty of the model, which were: 1) the above-mentioned $T^2$ function (see Section 2.2.1), which was normalized (for each eruption) by dividing it with the mean values of mass loading measured in the field (MML); 2) the percentage of sections for which there was an overestimation and an underestimation; 3) the mean overestimation (MO) and the mean underestimation (MU),
$$ egin{cases} MO = \frac{\sum_{i=1}^{N_o} \Delta_i}{N_o} & \text{for } \Delta_i > 0 \\ MU = \frac{\sum_{i=1}^{N_u} \Delta_i}{N_u} & \text{for } \Delta_i < 0 \end{cases} $$

where N_o and N_u are the number of sections with overestimation and underestimation, respectively; 4) the respective ratios of MO and MU with the mean mass loading values (MML) measured in the field.

With these four parameters the aim was to define, for each eruption and each meteo data, 1) the discrepancy between the observed data and the model (T²/MML – the normalization allows to compare T² from different eruptions), 2) whether the model tends mostly to overestimate or underestimate the observed data (% of sections under or overestimated), 3) the quantification of, respectively, the absolute model mean underestimation (MU) and mean overestimation (MO) and, 4) how important are MO and MU with respect to the mean values of mass loading measured in the field (MO/MML and MU/MML ratios). Regarding the grain size data, instead, the modeled values of Mdφ or σφ were plotted as a function of the observed values at specific stratigraphic sections, and the distribution of the data relative to a perfect fit line was discussed.

## 3. Results

For all the eruptions, Fig. 2 describes the stratigraphic sections used for uncertainty quantification, Figs. 3-6 provide the results of each comparison, while Tables 2-5 summarize the values calculated for each uncertainty quantification. Complementary data given in the Supporting Information are: the output values (plume heights, mass loadings, Mdφ and σφ values, Tables S4-S7) and the simulation outputs in PDF (Figures S1-S16).

### 3.1 Cotopaxi 2015

For the C15 eruption, a total of 35 mass loading measurements [from Bernard et al., 2016a] and 4 grain size analyses [unpublished and from Gaunt et al., 2016] were used for comparison with our model (Fig. 2a).

For each meteo condition and for the values of MER considered, plume heights comparison (Fig. 3a) shows that PLUME-MoM results are generally lower than those obtained by inverting seismic signal or from satellite/video camera images, though the model data mimic the patterns of observations. The difference between observed and modeled values (Table 2) is ~435-480 m for the seismic signal and video camera images while it is ~1300-1400 m for the satellite measurements. We note, however, a few exceptions. For the seismic-derived data, exceptions are the days around the 23rd of September, where modeled plume heights are systematically higher than the inferred ones. In contrast, Fig. 3a shows that there is a very good correlation between modeled and observed plume heights estimated from video recordings for the first phase of the eruption (August and beginning of September).

Ground deposits data show a difference of about 15°-20° between the directions of modeled and observed main dispersal axes (Figs. 2a and 3b). Notice that the deposits simulated, despite in extremely low quantities (i.e. 10⁻¹⁰-10⁻¹¹ kg/m²) at more distal locations, are spread all over the computational domain (see Figs. S1 to S3). Mass loading data show that the simulations underestimate field observations at locations in the main dispersal axes (Fig. 3b). Notice that the two sections along the main dispersal axes with the highest underestimations (BNAS and PNC 4 sections, see Table S3 from the Supporting Information) have observed mass loading values of, respectively, 18 and 15 kg/m²; for these two sections,
which are very proximal (~5 and ~7 km from the vent respectively) the model predicts very
low deposition (<1 kg/m² for all the simulations). The T²/MML values (Table 2) show that the
differences between model and observed values are relatively low, and the model generally
underestimates the observed values (57% to 77% of the field sections are underestimated). An
area of model underestimation might be recognized close to the vent area along the main
dispersal axes for all the simulations (see Figure S17 from the Supporting Information). The
MO and MU values (and also the MO/MML and MU/MML ratios) are similar for the
different meteo data, and for all the cases with a higher value of MU and MU/MML (for
simulations done using the NCEP/NCAR and the ERA-Interim meteo data).

The grain size data are scarce but we note that the computed Mdφ values are almost
always shifted toward coarser sizes (Fig. 3c) and that the σφ values show that the sorting of
the computed deposit is much smaller with respect to reality (Fig. 3d). Both computed Mdφ
and σφ show nearly constant values for a given section but with different meteo data.

414

415

## 3.2 Tungurahua 2013

For the T13 eruption, a total of 48 mass loading measurements [unpublished and from
Parra *et al.*, 2016] and 29 grain size analyses [unpublished and from Parra *et al.*, 2016] were
used for the comparison (Fig. 2b).

The plume heights comparison (Fig. 4a) shows that all the simulations markly
underestimate the observations reported from both sources. The mean difference is about -2.1
km to -2.2 km (Table 3). The difference of deposit main dispersal axes is small since the
simulations done using GDAS and ERA-Interim data are almost coincident with respect to
field data while the NCAR simulation is only 8° shifted toward the SW (Figs. 2b and 4b).

The observed values of mass loading (Fig. 4b and Table S3 from the Supporting
Information) are all <3 kg/m², similarly with respect to the C15 eruption for the two sections
along the main dispersal axes (San Pedro de Sabanag and 12 de Octubre, Table S3 from the
Supporting Information). Mass loading differences have a small spread highlighted by low
T²/MML values (Table 3). This is also shown by the absolute differences (MO and MU),
which are also almost identical despite the model tends to underestimate field data at most
sections. For the T13 eruption, the distribution of sections with overestimation and
underestimation does not highlight homogeneous areas of model overestimation or
underestimation (see Fig. S18 from the Supporting Information). In Table 3 the MO/MML
and MU/MML ratios have both values <1, indicating that the difference in mass loading value
is less important than the average deposit value of mass loading. In Fig. 4b, the mass loading
differences with respect to the observed data are equally positive (overestimation) or negative
(underestimation) in proximity of the main dispersal axes, without a clear prevalence.

Grain size comparison highlights that, similarly to the C15 eruption, most of the
computed grain sizes are shifted toward constant coarser grained values (Mdφ, see Fig. 4c)
with a smaller and fairly constant sorting for much of the sections (σφ, see Fig. 4d). Notice,
however, that some simulation sorting values are along the perfect fit line (mostly
NCEP/NCAR simulation) or are even larger than the observed ones (GDAS and the ERA-
Interim simulations).

443

444

## 3.3 Tungurahua 2006

For the T06 eruption, a total of 48 mass loading measurements [Eychenne *et al.*, 2012]
and 22 grain size analyses [recalculated from Eychenne *et al.*, 2012, see also Text S1 from the
Supporting Information] were used for the comparison (Fig. 2c).
448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

Fig. 5a shows that the plume heights simulated are close to observed data, except for the NCEP/NCAR model. The ERA-Interim/WRF model, in particular, provides a low mean overestimation of about 400 m (Table 4). Notice that this simulation was characterized by a fairly low $T^2$ value, although higher with respect to the parent ERA-Interim simulation (Table 4). This difference is due to the iterative procedure described in Section 2.4, which allowed finding the hourly values of MERs that minimized the differences in plume heights. Another combination of MERs was instead used for the other three meteorological datasets. Differences in deposit main dispersal axes are the highest of the four studied eruptions and are up to about 40° toward South (see ERA-Interim meteo in Figs. 2c and 5b).

With respect to mass loading, the $T^2$/MML values (Table 4) highlight a relatively high spread of the data, which is also reflected in the MO and MU values. In this case, it could be considered that most of the sections with underestimation are concentrated in proximity of the main dispersal axis highlighted by field data (Fig. 5b). Notice that the NCEP/NCAR provides the highest values of overestimation (MO = 62.57, MO/MML = 7.68). Moreover, the T06 eruption is one of the two cases, among the studied ones, where one simulation gives more sections with overestimation than sections with underestimation (ERA-Interim/WRF, see Table 4). Considering the spatial distribution of sections with overestimation and underestimation (see Fig. S19 from the Supporting Information), then a homogeneous area of model overestimation might be identified in the proximity of the vent area along the main dispersal axes (see Fig. S19 from the Supporting Information). Figure 5b highlights an interesting pattern for all the sections since the difference in mass loading tends to increase approaching the main dispersal axis, which is particularly evident for the GDAS and the NCEP/NCAR simulations.

The grain size data show a fairly well defined trend of $M d\phi$ values, which are close to the perfect fit line (Fig. 5c). The model sorting values are instead mostly shifted toward lower values but define trends mimicking that of the perfect fit line (Fig. 5d).

## 3.4 Puyehue-Cordón Caulle 2011

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

For the PCC11 eruption, a total of 75 mass loading measurements and 24 grain size analyses [Bonadonna et al., 2015a; Pistolesi et al., 2015; unpublished] were used for the comparison (Fig. 2d). For the mass loadings, the thickness data of Pistolesi et al. [2015] were multiplied by the bulk deposit density value of 560 kg/m³ reported in Bonadonna et al. [2015a] for Unit I, in order to obtain kg/m² values. Daily average plume heights a.s.l. reported in Bonadonna et al. [2015b] have been converted into “above vent” values by subtracting the vent elevation reported in Bonadonna et al. [2015b] (1470 m a.s.l.).

For this eruption, the simulations generally overestimate the plume heights observed, which are lowered with the inverse procedure (see Table 5, Fig. 6a). The simulated deposit main dispersal axes are all shifted toward the South by 5-10° with respect to the field data (Figs. 2d and 6b).

For the mass loading, most of the $T^2$/MML values are the highest among all the simulations, with values up to 22.12 (ERA-Interim) (Table 5). MO and MU values are respectively >100 kg/m² and from -18 kg/m² up to -54 kg/m². The MO/MML and MU/MML ratios indicate anyway that mean overestimation is 3 to 6 times higher than MML and that mean underestimation is 0.3 to 1 times higher than MML. As for the other eruptions, the percentage of sections with overestimation is lower than that with underestimation, except for the simulation done with the GDAS meteo data (Table 5). From Fig. S20 from the Supporting Information, the distribution of the sections with overestimation or underestimation highlights a homogeneous area of model overestimation located 30-40 km from vent area along the main dispersal axes. The correlation between high values of mass loading overestimation and the
position of the main dispersal axis (Figure 6b) is evident only for the simulation done with the ERA-Interim meteorological data. For the other simulations instead, the sections with the highest differences are uncorrelated with respect to the position of the main dispersal axis given by the model. It is also important to underline that in this case also, sections with highest values of observed mass loadings are not correlated with the deposit main dispersal axis given by field data, a pattern that is confirmed also by the simulations (see Fig. 6b). This latter feature might be correlated with the progressive anticlockwise rotation of the ash cloud, a pattern already discussed by *Pistolesi et al.* [2015] and *Bonadonna et al.* [2015b]. To confirm this, we have also performed a more detailed analysis using satellite images to track the evolution of the ash cloud during the 04-05/06/2011: details about this method are reported in Text S1 from the Supporting Information. The sequence of images derived (Figure S21 from the Supporting Information) show that at the onset of the eruption the cloud drifted southwestwardly (130°), but as time passed, the cloud rapidly moved towards the east, reaching 105°. This compares with the main dispersal axis assessed from the field deposits integrated over the whole Unit I (layers A-F) and yielding a mean direction of 117°. However, the maximum mass loading of deposits have been recorded at much higher angles, lying between 130-135° (Figure 6b). This actually correlates with ash emissions occurring at the onset of the eruption, where the ash-rich plume might have produced rapid and en masse fallouts along the main ash cloud dispersal axis centered at 130° (Figure S21 from the Supporting Information). This is supported by mass loading values of the deposits, which are very high on the dispersal axis (green dots in Figure S21), ranging from 481.6kg/m² close to the vent (section n° 57, Table S3 from the Supporting Information) to 160 kg/m² at a greater distance. By contrast, the mass loading of samples located away from the dispersal axis (red dots in Figure S21), shows much lower values of about 5.6 kg/m², although being close to the vent. Interestingly, section n° 57 is also the one that tends to have the highest value of underestimations (up to -400 kg/m²).

Regarding the grain size data, the Mdφ values are spread on both sides of the perfect fit line (Fig. 6c). The NCEP/NCAR simulations (both direct and inverse) tend to give finer grained values with respect to the observed data. The sorting data tend to define two trends of constant values of σφ ~0.5 and ~2, and some model sorting values are higher than the observed ones (Fig. 6d). An important remark for the modeled grain sizes of the PCC11 eruption is that none of them show any bimodal distribution in contrast to the observed data. This is particularly evident for the above-mentioned section n° 57, which does not have bimodality and which has an Mdφ shifted toward more coarser-grained values.

531

532

## 4. Discussion

### 4.1 Uncertainty in the input parameters

A significant amount of uncertainty in the simulations may derive from the meteorological data employed. As also shown by other studies [e.g. *Devenish et al.*, 2012; *Webster et al.*, 2012], even small errors in the wind field can lead to large errors in the ash concentration, making therefore a point-by-point comparison of modelled with observed data a challenging task. The datasets we considered are among the most widely used in similar numerical modellings [e.g., *Webley et al.*, 2009; *Bonasia et al.*, 2012; *Folch*, 2012; *Costa et al.*, 2016]: moreover, it has also been used the mesoscale meteorological model WRF/ARW, which has been coupled with other TTDMs in similar works [e.g. *FALL3D*, *Poret et al.*, 2017]. From our results, it is not evident that a particular meteorological dataset provides systematically the best results. For instance, the GDAS dataset provides the worst results (in terms of both the T² and the MO-MU values) for the lower magnitude C15 and T13 eruptions,
546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

578

579

580

581

582

583

584

585

586

587

588

589

590

591

592

593

594

595

while it provides the best results for the T06 and PCC11 eruptions. The NCEP/NCAR dataset shows the opposite as the results are better for the C15 and T13 eruptions with respect to T06 and PCC11. The employment of the WRF/ARW model (see also Text S1 from the Supporting Information) did not result in a significant improvement of the results as it gave instead higher T²/MML values with respect to the parent ERA-Interim meteorological file (see Table 4), although for some other models the employment of the WRF/ARW model gave better results [Parra et al., 2016]. Given the high computational times necessary to process original meteo data, the refinement procedure using WRF/ARW was not applied to other longer eruptions. The meteorological data have a considerable effect on the direction of main advection of the volcanic particles, which controls the deposit main dispersal axis direction. This is particularly evident for the T06 eruption, where differences with respect to the observed axis are up to 40°. Two main reasons for such differences may be invoked: i) the meteorological data are built in a way such that their parameters remain constant for relatively long periods (3 to 6 hours) and for quite large areas (0.75°x0.75° up to 2.5°x2.5°), and within such temporal frames and spatial domains it is not possible to capture the variability of natural phenomena; ii) 4-dimensions meteorological files (especially Reanalysis products) might be less accurate over complex terrains (e.g. the Andes), for which the details of the atmospheric flow are less likely captured and there are not a lot of observations available. This could be the case for the T06 and T13 eruptions, where the rugged topography of the area surrounding the Tungurahua volcano could have caused secondary atmospheric effects not recorded in the meteorological files.

A common problem with eruption source parameters is the measurements of plume height. For instance, for the C15 eruption Bernard et al. [2016a] used three different methodologies for plume height estimates (inversion of seismic signals, video cameras observations, and satellite measurements), which gave sometimes very different values (see Fig. 3a). For the T06 eruption, Steffke et al. [2010] used two different methods of satellite observations. Therefore, it is not surprising that differences in measurements at the same time can be important. The uncertainty in plume height is also high for the T13 eruption, for which two different methods (satellite measurements and visual observations) have been employed, and for the PCC11 eruption as well, for which only daily mean values of plume height have been reported.

Mass loading values for the C15, T13 and T06 eruptions have been actually measured for each section (with various methods), but for the PCC11 they have been determined by multiplying the deposit thickness by a mean bulk deposit density value (see Section 3.4). This latter aspect is critical since density of tephra fall deposits may vary considerably owing to drastic density change between different particle sizes [e.g., Bonadonna and Phillips, 2003; Eychenne and Le Pennec, 2012; Pistolesi et al., 2015]. This is particularly important for the PCC11 eruption that has the highest T²/MML values (see Table 5), which might also be related to an uncertainty in the observed mass loading data. We also stress that the assumption of a linear variation of particle density with grain size (employed in PLUME-MoM) is a simplification since the density variation may be more complex [i.e sigmoidal rather than linear as for the T06 eruption, Eychenne and Le Pennec, 2012]. Compared to other sources of uncertainty, however, the simplification used in the simulations is expected to have a minor effect on the final results.

Finally, it is important to remark that there are also uncertainties in estimations of the initial water mass fraction in magmas. This is due primarily to the use of different methods [e.g., by direct measurements, geological inference, thermodynamic calculation or experimental approaches, see Clemens, 1984], among which the direct measurement from melt inclusions in crystals are the most used [see for example Plank et al., 2013]. As a comparison, for this study we relied on estimates made both using direct measurements from
melt inclusions and experimental approaches [Martel et al., 2018] or considering only
experimental approaches [Andújar et al., 2017]: results gave H₂O wt. % ranging between 4-6
wt. % and 6-8 wt. % for andesites and rhyolites respectively. As the water mass fraction has a
strong influence on the plume height simulated with PLUME-MoM [see section 2.2.1 and
also de'Michieli Vitturi et al., 2016], its careful estimation is therefore of primary importance.

## 4.2 Uncertainty in the numerical modelling

When MER values obtained from total deposit measurements are used as input
parameters, PLUME-MoM underestimates the plume height measurements for three out of
four eruptions tested, and there may be two main reasons for that. First, as already discussed
in the previous section, the measurements are in some specific cases uncertain. Second, the
mass eruption rate, assumed to be equal to the total mass of deposit divided with the eruption
duration, may be underestimated in some cases (e.g. the T06 eruption) since deposits of
pyroclastic density currents are neglected, hence giving lower plume heights. We note,
however, that the mean underestimations (and mean overestimations as well) of the model for
each eruption are lower with respect to the uncertainty in observed data among different
methods, and that in some cases (e.g. the T06 eruption) the refinement of the meteorological
data using the WRF/ARW model can sensibly reduce the difference in plume height with
respect to observed data.

The PLUME-MoM/HYSPLIT model tends generally to have more points
underestimating the mass loading data (see Tables 2 to 5). However, if the absolute mean
differences (MO and MU) and their ratios with mean values of mass loading (MO/MML and
MU/MML) are considered, then model overestimation is systematically higher with respect to
underestimation. For example, for the PCC11 eruption and for the simulation done using the
ERA-Interim data, MO is almost 10 times higher than MU (Table 5). The high values of MO
or MU and of their ratios with MML tend also to be higher for higher magnitude eruptions
(e.g. T06 and PCC11): in this regard the inverse procedure reduces considerably the
discrepancy between modeled and observed data as indicated for instance by the T²/MML
value for the PCC11 eruption.

The problem of model uncertainty is further illustrated by the difference in mass
loading with respect to the orientation of the stratigraphic section (Figs. 3b, 4b, 5b and 6b).
There are two opposite situations since the deposit main dispersal axis coincides either with
the lowest values of Δ mass loading (highest underestimation, e.g. C15 eruption, Fig. 3b) or
with the highest values of Δ mass loading (highest overestimation, T06 eruption, Fig. 7b, and
to a lesser extent T13 and PCC11 eruptions). This may be explained considering the advective
and diffusive parts of the transport equation used [Folch, 2012]. While the mass seems to be
correctly advected in the simulations (although with some deviation with respect to observed
data), the equations of HYSPLIT related to turbulent diffusion do not appear to work
efficiently, underestimating the horizontal diffusion and concentrating the mass close to the
main dispersal axis of advection. A similar issue has been also encountered by Hurst and
Davis [2017]. This may explain the above-mentioned mass loading underestimation or
overestimation, which are possibly increased by the fact that the HYSPLIT model does not
account for complex collective settling mechanisms of volcanic ash caused by aggregation,
gravitational instabilities, diffusive convection, particle-particle interactions and wake-capture
effects [Del Bello et al., 2017; Gouhier et al., 2019]. However, the problem of the effect of
diffusion on volcanic plumes dispersal and therefore on particle sedimentation is complex
[see for example Devenish et al., 2012]: a more rigorous study is therefore needed for
HYSPLIT to investigate the influence of different available diffusion equations on final
results.
645

646

647

648

649

650

651

652

653

654

655

656

657

The failure to take into account such mechanisms implies that the simulated finest-
grained particles are transported much further that in reality. For instance, the C15 eruption
has a particularly fine-grained TGSD [due also to its hydrovolcanic nature, Bernard *et al.*,
2016a] (see Table S2 from the Supporting Information) so that the mass is transported all over
the computational domain (see Figs. S1 to S3 from the Supporting Information). The case of
the PCC11 eruption is similar since the TGSD is up to 12φ, and an estimated amount of ~5%
of the erupted mass is transported out of the computational domain. While for this eruption
the finest fraction of the volcanic clouds circumvented the Southern hemisphere and passed
over the South of Australia [Collini *et al.*, 2013], it is possible that part of the fine ash did not
deposit (see also the issue of grain size analyses in the following paragraph). In this context,
the transport of material could have been at its maximum along the main dispersal axes, and
therefore the degree of underestimation of mass loading at proximal-medial sites along
dispersal axes is maximized as well.

658

659

660

661

662

663

664

665

666

667

668

669

Regarding the simulated grain size data, the Mdφ values are systematically coarser-
grained for the C15 and T13 low magnitude eruptions while they are either coarser-grained or
finer-grained for the PCC11 eruption. The shifting toward coarser-grained Mdφ values can be
explained by the fact that HYSPLIT neglects the above-mentioned collective settling
mechanisms of volcanic ash. For the eruptions where the amount of fine ash is higher (the
C15, T13 eruptions and partially the PCC11 one), the fine ash is transported distally, hence
causing coarser grain sizes in proximal to medial sections. Moreover, the model is not capable
of reproducing the bimodality of grain size distribution observed, as for instance in the PCC11
eruption. The σφ comparisons show that, instead, for most cases the modeled data tend to
have a lower sorting value with respect to the observed ones. These results show that the
employment of grain size data for model validation is less reliable with respect to mass
loading data.

670

671

672

673

674

675

676

677

678

679

680

681

682

683

684

685

686

687

688

689

690

691

692

693

694

Four important issues should be considered to improve the coupled PLUME-MoM/HYSPLIT model in the context of tephra fallout hazard assessments and probabilistic hazard maps production. First, the meteorological dataset must be considered carefully since it controls strongly the plume height. Second, the amount of fine ash and the duration of the eruption seem to be more critical than the magnitude of the eruption for mass loading calculations, since the simulations of higher magnitude eruptions of short duration with lower wt% of fine particles (i.e. T06 eruption) are more accurate than simulations of lower magnitude eruptions with longer durations and a higher amount of fines (i.e. the C15 and T13). If the magnitude, the amount of fine particles and the duration of the eruption are high (i.e. the PCC11 eruption), then the model tends to overestimate the natural data. Third, for the above-mentioned reasons, we recommend to employ PLUME-MoM/HYSPLIT in its present configuration for the production of hazard maps related to higher magnitude eruptions (i.e. sub-Plinian or Plinian). This is supported by our simulations of such eruptions (i.e. T06 and PCC11), for which overestimation is much higher (in terms of mean absolute values) with respect to underestimation. This latter point is important in a context of hazard assessment since underestimation may be considered as less acceptable than overestimation. Moreover, it is also important to remind that: a) specifically for our test eruptions, the lower magnitude ones tend to have longer durations and are more difficult to model due to the very high variability of both the eruptions parameters and atmospheric conditions, which are less likely to be captured; b) the T06 and PCC11 eruptions are those for which modeled and observed plume heights are more similar. Fourth, the MO/MML and MU/MML ratios may be used to account for model uncertainty and to serve as a basis for calculating coefficients that allow the creation of probabilistic maps (from the point of view of mass loading) that quantify the model mean overestimation and underestimation. For this purpose, statistical techniques might be employed to correct the model by estimating its deviance from the observed data.
695 **5. Conclusions**

696 This paper presents the coupling of the PLUME-MoM model with a renewed version
697 of the HYSPLIT tephra dispersal model. These two coupled models have been tested against
698 four eruptions of different magnitudes and styles from three Andean volcanoes. A procedure
699 of uncertainty quantification has been applied by computing the differences between modeled
700 and observed data of plume height, mass loading and grain size (in terms of Mdφ and σφ).
701 Four different meteorological datasets (GDAS, NCAR/NCEP, ERA-Interim, ERA-
702 Interim+WRF) have been tested as well. The main conclusions and future perspectives of this
703 uncertainty quantification are:

704 • None of the meteorological datasets tested produced systematically the best results
705 for all the eruptions. This implies that if a specific dataset is employed for
706 numerical modelling, its uncertainty (as quantified here) should be considered.

707 • The PLUME-MoM model tends to underestimate measured plume heights, except
708 for the eruption with the highest magnitude tested (i.e. PCC11). Though
709 underestimation might be related to uncertainties in plume height measurements
710 and input data (e.g. mass eruption rate and initial water mass fraction), we note
711 that for most of the cases we investigated the mean underestimations of the model
712 for each eruption were lower with respect to the uncertainty in observed data.

713 • The PLUME-MoM/HYSPLIT model tends generally to have more points
714 underestimating the mass loading data. If absolute mean differences are
715 considered, however, then overestimation is almost always higher than
716 underestimation. The distribution of sections with overestimation and
717 underestimation does not highlight systematically homogeneous areas of either
718 overestimation or underestimation.

719 • The advective part of the HYSPLIT model appears to work more efficiently than
720 the diffusive part. Moreover, the failure to take into account any collective settling
721 mechanisms of volcanic ash in HYSPLIT might cause important discrepancies
722 between observed and modeled data of mass loading and, above all, grain size
723 distributions.

724 • For the above-mentioned reasons, high amounts of fine particles might reduce the
725 accuracy of the model when simulating mass loading and grain size data.

726 • If the PLUME-MoM/HYSPLIT model is meant to be employed for hazard
727 assessment purposes, we recommend considering high magnitude eruptions (i.e.
728 sub-Plinian or Plinian) as target cases, and mass loading as primary parameter.

729 • Future developments of this project should consider the comparison of simulations
730 outputs with those from other models, in order to identify which model is best
731 suited for a specific eruption type

732 **Acknowledgements**

733 This research was financed by the French government IDEX-ISITE initiative 16-IDEX-0001
734 (CAP 20-25), the French Research Institute for Sustainable Development (IRD) in the context
735 of the Laboratoire Mixte International "Séismes et Volcans dans les Andes du Nord" and the
736 CNRS Tellus programme. Marco Pistolesi is warmly acknowledged for providing additional
737 unpublished data for the Puyehue-Cordón Caulle 2011 eruption. The authors are grateful to
738 Sandra Banson and Cèline Planche for the help in setting up the simulations with WRF. We
739 thank two anonymous reviewers for useful comments that improved the quality of the
740 manuscript, and we acknowledge as well the editorial handling provided by Nico Fournier
741 and Martha Savage.
742 Data supporting the analyses and conclusions presented in this study could be found in a
743 repository in the Figshare community with the following DOI:
744 https://doi.org/10.6084/m9 figshare.10298036.v1.

745

746 References

747 Andújar, J., C. Martel, M. Pichavant, P. Samaniego, B. Scaillet, and I. Molina, 2017, Structure of the plumbing system at Tungurahua volcano, Ecuador: insights from phase equilibrium experiments on July–August 2006 eruption products, *J. Petrol.*, 58(7), 1249-1278, doi: https://doi.org/10.1093/petrology/egx054.

751 Bellotti, F., L. Capra, D. Sarocchi, and M. D'Antonio, 2010, Geostatistics and multivariate analysis as a tool to characterize volcaniclastic deposits: Application to Nevado de Toluca volcano, Mexico, *J. Volcanol. Geoth. Res.*, 191(1-2), 117-128, doi: https://doi.org/10.1016/j.jvolgeores.2010.01.005.

755 Bernard, B., J. Battaglia, A. Proaño, S. Hidalgo, F. Vásconez, S. Hernandez, and M. C. Ruiz, 2016a, Relationship between volcanic ash fallouts and seismic tremor: quantitative assessment of the 2015 eruptive period at Cotopaxi volcano, Ecuador, *Bull. Volc.*, 78(11), doi: 10.1007/s00445-016-1077-5.

759 Bernard, J., K. Kelfoun, J. L. Le Pennec, and S. V. Vargas, 2014, Pyroclastic flow erosion and bulking processes: comparing field-based vs. modeling results at Tungurahua volcano, Ecuador, *Bull. Volc.*, 76(9), 858, doi: 10.1007/s00445-014-0858-y.

762 Bernard, J., J. Eychenne, J. L. Le Pennec, and D. Narváez, 2016b, Mass budget partitioning during explosive eruptions: insights from the 2006 paroxysm of Tungurahua volcano, Ecuador, *Geochem. Geophy. Geosy.*, 17(8), 3224-3240, doi: 10.1002/2016GC006431.

765 Berrisford, P., D. Dee, P. Poli, R. Brugge, K. Fielding, M. Fuentes, P. Kallberg, S. Kobayashi, S. Uppala, and A. Simmons, 2011, The ERA-Interim archive, version 2.0.

767 Biass, S., and C. Bonadonna, 2014, TOTGS: total grainsize distribution of tephra fallout: https://vhub.org/resources/3297.

769 Biondi, R., A. K. Steiner, G. Kirchengast, H. Brenot, and T. Rieckh, 2017, Supporting the detection and monitoring of volcanic clouds: A promising new application of Global Navigation Satellite System radio occultation, *Adv. Sp. Res.*, 60(12), 2707-2722, doi: https://doi.org/10.1016/j.asr.2017.06.039.

773 Bonadonna, C., and J. C. Phillips, 2003, Sedimentation from strong volcanic plumes, *J. Geophys. Res.-Sol. Ea.*, 108(B7), 2340, doi: 10.1029/2002JB002034.

775 Bonadonna, C., R. Cioni, M. Pistolesi, M. Elissondo, and V. Baumann, 2015a, Sedimentation of long-lasting wind-affected volcanic plumes: the example of the 2011 rhyolitic Cordón Caulle eruption, Chile, *Bull. Volc.*, 77(2), doi: https://doi.org/10.1007/s00445-015-0900-8.

779 Bonadonna, C., M. Pistolesi, R. Cioni, W. Degruyter, M. Elissondo, and V. Baumann, 2015b, Dynamics of wind-affected volcanic plumes: The example of the 2011 Cordón Caulle eruption, Chile, *J. Geophys. Res.-Sol. Ea.*, 120(4), 2242-2261, doi: 10.1002/2014JB011478.

783 Bonasia, R., G. Macedonio, A. Costa, D. Mele, and R. Sulpizio, 2010, Numerical inversion and analysis of tephra fallout deposits from the 472 AD sub-Plinian eruption at Vesuvius (Italy) through a new best-fit procedure, *J. Volcanol. Geoth. Res.*, 189(3-4), 238-246, doi: https://doi.org/10.1016/j.jvolgeores.2009.11.009.

787 Bonasia, R., A. Costa, A. Folch, G. Macedonio, and L. Capra, 2012, Numerical simulation of tephra transport and deposition of the 1982 El Chichón eruption and implications for
789
790

hazard assessment, *J. Volcanol. Geoth. Res.*, 231(2012), 39-49, doi:
https://doi.org/10.1016/j.jvolgeores.2012.04.006.

791
792

Bursik, M. I., 2001, Effect of wind on the rise height of volcanic plumes, *Geophys. Res. Lett.*, 28(18), 3621-3624, doi: https://doi.org/10.1029/2001GL013393.

793
794
795

Caballero, L., D. Sarocchi, E. Soto, and L. Borselli, 2014, Rheological changes induced by clast fragmentation in debris flows, *J. Geophys. Res.-Earth*, 119(9), 1800-1817, doi: 10.1002/2013JF002942.

796
797

Clemens, J. D., 1984, Water contents of silicic to intermediate magmas, *Lithos*, 17(273-287, doi: https://doi.org/10.1016/0024-4937(84)90025-2.

798
799
800

Collini, E., M. S. Osores, A. Folch, J. G. Viramonte, G. Villarosa, and G. Salmuni, 2013, Volcanic ash forecast during the June 2011 Cordón Caulle eruption, *Nat. Hazards*, 66(2), 389-412, doi: https://doi.org/10.1007/s11069-012-0492-y.

801
802
803
804
805

Connor, L. J., and C. B. Connor, 2006, Inversion is the key to dispersion: understanding eruption dynamics by inverting tephra fallout, in Mader, H. M., Coles, S.G., Connor, C.B., Connor, L.J, Statistics in Volcanology, IAVCEI Special Publication, The Geological Society of London, p. 231-241, doi: https://doi.org/10.1144/IAVCEI001.18.

806
807
808

Costa, A., F. Dell'Erba, M. A. Di Vito, R. Isaia, G. Macedonio, G. Orsi, and T. Pfeiffer, 2009, Tephra fallout hazard assessment at the Campi Flegrei caldera (Italy), *Bull. Volc.*, 71(3), 259, doi: https://doi.org/10.1007/s00445-008-0220-3.

809
810
811
812
813
814
815

Costa, A., Y. J. Suzuki, M. Cerminara, B. J. Devenish, T. Esposti Ongaro, M. Herzog, A. R. Van Eaton, L. C. Denby, M. I. Bursik, M. de'Michieli Vitturi, S. Engwell, A. Neri, S. Barsotti, A. Folch, G. Macedonio, F. Girault, G. Carazzo, S. Tait, E. Kaminski, L. G. Mastin, M. J. Woodhouse, J. C. Phillips, A. J. Hogg, W. Degruyter, and C. Bonadonna, 2016, Results of the eruptive column model inter-comparison study, *J. Volcanol. Geoth. Res.*, 326(2016), 2-25, doi: https://doi.org/10.1016/j.jvolgeores.2016.01.017.

816
817
818

de'Michieli Vitturi, M., A. Neri, and S. Barsotti, 2015, PLUME-MoM 1.0: A new integral model of volcanic plumes based on the method of moments, *Geosci. Model Dev.*, 8(8), 2447, doi: https://doi.org/10.5194/gmd-8-2447-2015.

819
820
821
822

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

823
824
825
826

Dee, D. P., S. M. Uppala, A. J. Simmons, P. Berrisford, P. Poli, S. Kobayashi, U. Andrae, M. A. Balmaseda, G. Balsamo, and D. P. Bauer, 2011, The ERA-Interim reanalysis: Configuration and performance of the data assimilation system, *Q. J. Roy. Meteor. Soc.*, 137(656), 553-597, doi: https://doi.org/10.1002/qj.828.

827
828
829
830

Del Bello, E., J. Taddeucci, M. de' Michieli Vitturi, P. Scarlato, D. Andronico, S. Scollo, U. Kueppers, and T. Ricci, 2017, Effect of particle volume fraction on the settling velocity of volcanic ash particles: insights from joint experimental and numerical simulations, *Sci. Rep.-UK*, 7(39620, doi.

831
832
833

Devenish, B. J., P. N. Francis, B. T. Johnson, R. S. J. Sparks, and D. J. Thomson, 2012, Sensitivity analysis of dispersion modeling of volcanic ash from Eyjafjallajökull in May 2010, *J. Geophys. Res.-Atmos.*, 117(D00U21), doi: 10.1029/2011JD016782.

834
835
836
837
838

Douillet, G. A., È. Tsang-Hin-Sun, U. Kueppers, J. Letort, D. A. Pacheco, F. Goldstein, F. Von Aulock, Y. Lavallée, J. B. Hanson, and J. Bustillos, 2013, Sedimentology and geomorphology of the deposits from the August 2006 pyroclastic density currents at Tungurahua volcano, Ecuador, *Bull. Volc.*, 75(11), 765, doi: https://doi.org/10.1007/s00445-013-0765-7.
839

840

841

Engwell, S., and J. Eychenne, 2016, Contribution of fine ash to the atmosphere from plumes associated with pyroclastic density currents, Volcanic Ash, Elsevier, p. 67-85, doi: https://doi.org/10.1016/B978-0-08-100405-0.00007-0.

842

843

844

Eychenne, J., and J. L. Le Pennec, 2012, Sigmoidal particle density distribution in a subplinian scoria fall deposit, *Bull. Volc.*, 74(10), 2243-2249, doi: https://doi.org/10.1007/s00445-012-0671-4.

845

846

847

848

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

849

850

851

852

Eychenne, J., J. L. Le Pennec, P. Ramon, and H. Yepes, 2013, Dynamics of explosive paroxysms at open-vent andesitic systems: high-resolution mass distribution analyses of the 2006 Tungurahua fall deposit (Ecuador), *Ear. Pl. Sci. Lett.*, 361(343-355, doi: https://doi.org/10.1016/j.epsl.2012.11.002.

853

854

Folch, A., A. Costa, and G. Macedonio, 2009, FALL3D: A computational model for transport and deposition of volcanic ash, *Comp. & Geos.*, 35(6), 1334-1342, doi.

855

856

857

Folch, A., 2012, A review of tephra transport and dispersal models: evolution, current status, and future perspectives, *J. Volcanol. Geoth. Res.*, 235(96-115, doi: https://doi.org/10.1016/j.jvolgeores.2012.05.020.

858

859

860

Folch, A., A. Costa, and G. Macedonio, 2016, FPLUME-1.0: An integral volcanic plume model accounting for ash aggregation, *Geosci. Model Dev.*, 9(431-450, doi: 10.5194/gmd-9-431-2016.

861

862

863

Folk, R. L., and W. C. Ward, 1957, Brazos River bar [Texas]; a study in the significance of grain size parameters, *J. Sediment. Res.*, 27(1), 3-26, doi: https://doi.org/10.1306/74D70646-2B21-11D7-8648000102C1865D.

864

865

866

867

García Moreno, J. D., 2016, Mapeo y determinación de parámetros físicos de las corrientes de densidad piroclásticas producidas por el volcán Tungurahua producidas desde el 2006, Tesis de Ingeniería: Escuela Politécnica Nacional, 152 p, doi: http://bibdigital.epn.edu.ec/handle/15000/15075.

868

869

870

871

872

Gaunt, H. E., B. Bernard, S. Hidalgo, A. Proaño, H. Wright, P. A. Mothes, E. Criollo, and U. Kueppers, 2016, Juvenile magma recognition and eruptive dynamics inferred from the analysis of ash time series: the 2015 reawakening of Cotopaxi volcano, *J. Volcanol. Geoth. Res.*, 328(2016), 134-146, doi: https://doi.org/10.1016/j.jvolgeores.2016.10.013.

873

874

875

876

Gouhier, M., J. Eychenne, N. Azzaoui, A. Guillin, M. Deslandes, M. Poret, A. Costa, and P. Husson, 2019, Low efficiency of large volcanic eruptions in transporting very fine ash into the atmosphere, *Sci. Rep.-UK*, 9(1), 1449, doi: https://doi.org/10.1038/s41598-019-38595-7.

877

878

879

880

Hall, M. L., A. L. Steele, P. A. Mothes, and M. C. Ruiz, 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(2013), 78-93, doi: https://doi.org/10.1016/j.jvolgeores.2013.08.011.

881

882

883

884

Hidalgo, S., J. Battaglia, S. Arellano, A. L. Steele, B. Bernard, J. Bourquin, B. Galle, S. Arrais, and F. Vásconez, 2015, SO₂ degassing at Tungurahua volcano (Ecuador) between 2007 and 2013: Transition from continuous to episodic activity, *J. Volcanol. Geoth. Res.*, 298(1-14, doi: https://doi.org/10.1016/j.jvolgeores.2015.03.022.

885

886

887

888

Hidalgo, S., J. Battaglia, S. Arellano, D. Sierra, B. Bernard, R. Parra, P. Kelly, F. Dinger, C. Barrington, and P. Samaniego, 2018, Evolution of the 2015 Cotopaxi eruption revealed by combined geochemical and seismic observations, *Geochem. Geophy. Geosy.*, 19(7), 2087-2108, doi: https://doi.org/10.1029/2018GC007514.
889

Hurst, T., and C. Davis, 2017, Forecasting volcanic ash deposition using HYSPLIT, *J. Appl. Volcan.* 6(1), 5, doi: https://doi.org/10.1186/s13617-017-0056-7.

890

891

Jay, J., F. Costa, M. Pritchard, L. Lara, B. Singer, and J. Herrin, 2014, Locating magma reservoirs using InSAR and petrology before and during the 2011–2012 Cordón Caulle silicic eruption, *Ear. Pl. Sci. Lett.*, 395(2014), 254-266, doi: https://doi.org/10.1016/j.epsl.2014.03.046.

892

893

894

895

Kalnay, E., M. Kanamitsu, R. Kistler, W. Collins, D. Deaven, L. Gandin, M. Iredell, S. Saha, G. White, J. Woollen, Y. Zhu, M. Chelliah, W. Ebisuzaki, W. Higgins, J. Janowiak, K. C. Mo, C. Ropelewski, J. Wang, A. Leetmaa, R. Reynolds, R. Jenne, and D. Joseph, 1996, The NCEP/NCAR 40-Year Reanalysis Project, *Bull. Am. Met. Soc.*, 77(3), 437-472, doi: 10.1175/1520-0477(1996)077<0437:tnyrp>2.0.co;2.

896

897

898

899

900

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

901

902

903

904

Martel, C., J. Andújar, P. Mothes, B. Scaillet, M. Pichavant, and I. Molina, 2018, Storage conditions of the mafic and silicic magmas at Cotopaxi, Ecuador, *J. Volcanol. Geoth. Res.*, 354(74-86, doi: https://doi.org/10.1016/j.jvolgeores.2018.02.006.

905

906

907

Marti, A., A. Folch, O. Jorba, and Z. Janjic, 2017, Volcanic ash modeling with the online NMMB-MONARCH-ASH v1.0 model: model description, case simulation, and evaluation, *Atmos. Chem. Phys.*, 17(2017), 4005-4030, doi: 10.5194/acp-17-4005-2017.

908

909

910

911

Michalakes, J., S. Chen, J. Dudhia, L. Hart, J. Klemp, J. Middlecoff, and W. Skamarock, 2001, Development of a next-generation regional weather research and forecast model, Developments in Teracomputing, World Scientific, p. 269-276, doi: https://doi.org/10.1142/9789812799685_0024.

912

913

914

915

NOAA, 2004, Global Data Assimilation System (GDAS1) Archive Information.

916

Parra, R., B. Bernard, D. Narváez, J. L. Le Pennec, N. Hasselle, and A. Folch, 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(2016), 1-13, doi: https://doi.org/10.1016/j.jvolgeores.2015.11.001.

917

918

919

920

921

Pistolesi, M., R. Cioni, C. Bonadonna, M. Elissondo, V. Baumann, A. Bertagnini, L. Chiari, R. Gonzales, M. Rosi, and L. Francalanci, 2015, Complex dynamics of small-moderate volcanic events: the example of the 2011 rhyolitic Cordón Caulle eruption, Chile, *Bull. Volc.*, 77(1), 3, doi: https://doi.org/10.1007/s00445-014-0898-3.

922

923

924

925

Plank, T., K. A. Kelley, M. M. Zimmer, E. H. Hauri, and P. J. Wallace, 2013, Why do mafic arc magmas contain 4 wt% water on average?, *Ear. Pl. Sci. Lett.*, 364(168-179, doi: https://doi.org/10.1016/j.epsl.2012.11.044.

926

927

928

Poret, M., A. Costa, A. Folch, and A. Martí, 2017, Modelling tephra dispersal and ash aggregation: The 26th April 1979 eruption, La Soufrière St. Vincent, *J. Volcanol. Geoth. Res.*, 347(2017), 207-220, doi: https://doi.org/10.1016/j.jvolgeores.2017.09.012.

929

930

931

932

Riley, C. M., W. I. Rose, and G. J. S. Bluth, 2003, Quantitative shape measurements of distal volcanic ash, *J. Geophys. Res.-Sol. Ea.*, 108(B10), doi: https://doi.org/10.1029/2001JB000818.

933

934

935

Samaniego, P., J. L. Le Pennec, C. Robin, and S. Hidalgo, 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, doi: https://doi.org/10.1016/j.jvolgeores.2010.10.010.

936

937

938
939

940

941

Scollo, S., A. Folch, and A. Costa, 2008, A parametric and comparative study of different tephra fallout models, *J. Volcanol. Geoth. Res.*, 176(2), 199-211, doi: https://doi.org/10.1016/j.jvolgeores.2008.04.002.

942

943

944

Skamarock, W. C., and J. B. Klemp, 2008, A time-split nonhydrostatic atmospheric model for weather research and forecasting applications, *J. Comput. Phys.*, 227(7), 3465-3485, doi: https://doi.org/10.1016/j.jcp.2007.01.037.

945

946

947

948

Steffke, A. M., D. Fee, M. Garces, and A. Harris, 2010, Eruption chronologies, plume heights and eruption styles at Tungurahua Volcano: Integrating remote sensing techniques and infrasound, *J. Volcanol. Geoth. Res.*, 193(3-4), 143-160, doi: https://doi.org/10.1016/j.jvolgeores.2010.03.004.

949

950

951

Stein, A. F., R. R. Draxler, G. D. Rolph, B. J. B. Stunder, M. D. Cohen, and F. Ngan, 2015, NOAA's HYSPLIT atmospheric transport and dispersion modeling system, *Bull. Am. Met. Soc.*, 96(12), 2059-2077, doi: https://doi.org/10.1175/BAMS-D-14-00110.1.

952

953

954

Tuffen, H., M. R. James, J. M. Castro, and C. I. Schipper, 2013, Exceptional mobility of an advancing rhyolitic obsidian flow at Cordón Caulle volcano in Chile, *Nat. Commun.*, 4(2013), 2709, doi: https://doi.org/10.1038/ncomms3709.

955

956

957

958

959

Webley, P. W., B. J. B. Stunder, and K. G. Dean, 2009, Preliminary sensitivity study of eruption source parameters for operational volcanic ash cloud transport and dispersion models—A case study of the August 1992 eruption of the Crater Peak vent, Mount Spurr, Alaska, *J. Volcanol. Geoth. Res.*, 186(1-2), 108-119, doi: https://doi.org/10.1016/j.jvolgeores.2009.02.012.

960

961

962

963

964

965

Webster, H. N., D. J. Thomson, B. T. Johnson, I. P. C. Heard, K. Turnbull, F. Marenco, N. I. Kristiansen, J. Dorsey, A. Minikin, B. Weinzierl, U. Schumann, R. S. J. Sparks, S. C. Loughlin, M. C. Hort, S. J. Leadbetter, B. J. Devenish, A. J. Manning, C. S. Witham, J. M. Haywood, and B. W. Golding, 2012, Operational prediction of ash concentrations in the distal volcanic cloud from the 2010 Eyjafjallajökull eruption, *J. Geophys. Res.-Atmos.*, 117(D00U08), doi: 10.1029/2011JD016790.

966

967

968

Wilson, L., and T. C. Huang, 1979, The influence of shape on the atmospheric settling velocity of volcanic ash particles, *Ear. Pl. Sci. Lett.*, 44(2), 311-324, doi: https://doi.org/10.1016/0012-821X(79)90179-1.

969

970

971

972

973

974

975

976

977

978

979

980

981
982

983

984

985

986

987

988

989

990

991

992

993

994

995

FIGURE CAPTIONS

997

998 **Figure 1.** Geographical locations of a) Cotopaxi and Tungurahua volcanoes in Ecuador and b) Puyehue-Cordón Caulle volcanic complex in Chile. Coordinates are in the UTM WGS84 17S

999 (a) and UTM WGS84 19S (b) systems. Basemap copyright of ESRI®, DigitalGlobe, GeoEye,

1000 Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN, and the GIS User

1001 Community.

1002

1003

1004 **Figure 2.** Stratigraphic sections used for uncertainty quantification (considering only mass

1005 loading or both mass loading and grain size) and dispersal axes from field data/simulations

1006 with different meteo data for: a) Cotopaxi C15 eruption; b) Tungurahua T13 eruption; c) Tungurahua T06 eruption; d) Puyehue-Cordón Caulle PCC11 eruption. Digital Elevation

1007 Model (30 m- resolution) from Marc Souris, IRD (a-c) and ESRI, USGS, NOAA (d).

1008

1009 Coordinates are in UTM WGS84 17S (a-c) and UTM WGS84 19S (d).

1010

1011 **Figure 3.** Cotopaxi C15 eruption. Comparison of a) column height (observed and computed);

1012 b) Δ mass loading (left axis) and observed mass loading (right axis) for different sections with

1013 respect to their orientation from North; c) Mdφ and d) σφ (computed and observed).

1014

1015 **Figure 4.** Tungurahua T13 eruption. a) Column height (observed and computed); b) Δ mass

1016 loading (left axis) and observed mass loading (right axis) for different sections with respect to

1017 their orientation from North; c) Mdφ and d) σφ (computed and observed).

1018
**Figure 5.** Tungurahua T06 eruption. a) column height (observed and computed); b) Δ mass loading (left axis) and observed mass loading (right axis) for different sections with respect to their orientation from North; c) Mdφ and d) σφ (computed and observed).

**Figure 6.** Puyehue-Cordón Caulle PCC11 eruption. a) Column height (observed and modeled); b) Δ mass loading (left axis) and observed mass loading (right axis) for different sections with respect to their orientation from North; c) Mdφ and d) σφ (computed and observed).

TABLES

<table><caption>**Table 1.** Main input parameters used for the simulations.</caption><thead><tr><th>Eruption name</th><th>Computational grid dimension (°)</th><th>Initial water content (wt%)</th><th>Particles exit velocity (m/s)</th><th>Heat capacity (J/kgxK)</th><th>Particles Shape Factor</th><th>Φ<sub>1</sub></th><th>ρ<sub>1</sub> (kg/m<sup>3</sup>)</th><th>Φ<sub>2</sub></th><th>ρ<sub>2</sub> (kg/m<sup>3</sup>)</th></tr></thead><tbody><tr><td>Cotopaxi 2015 (C15)</td><td>5x5</td><td>5.5%</td><td>135</td><td>1600</td><td>0.75</td><td>-1</td><td>1487</td><td>2</td><td>2478</td></tr><tr><td>Tungurahua 2013 (T13)</td><td>6 x 6</td><td>5.5%</td><td>135</td><td>1600</td><td>0.75</td><td>-1</td><td>1487</td><td>2</td><td>2478</td></tr><tr><td>Tungurahua 2006 (T06)</td><td>6 x 6</td><td>5.5%</td><td>275</td><td>1600</td><td>0.75</td><td>-1</td><td>1487</td><td>2</td><td>2478</td></tr><tr><td>Puyehue-Cordón Caulle 2011 (PCC11)</td><td>10 x 10</td><td>7.0%</td><td>275</td><td>1600</td><td>0.65</td><td>-4</td><td>500</td><td>5</td><td>2670</td></tr></tbody></table>

<table><thead><tr><th rowspan="2">Parameters</th><th colspan="3">Meteo Data</th></tr><tr><th>GDAS</th><th>NCEP/NCAR</th><th>ERA-Interim</th></tr></thead><tbody><tr><td colspan="4"><em>PLUME-MoM</em></td></tr><tr><td>Mean Difference seismic (m)</td><td>-479.02</td><td>-443.75</td><td>-454.83</td></tr><tr><td>Mean Difference video (m)</td><td>-466.24</td><td>-434.99</td><td>-443.23</td></tr><tr><td>Mean Difference satellite (m)</td><td>-1449.64</td><td>-1405.96</td><td>-1366.64</td></tr><tr><td colspan="4"><em>HYSPLIT</em></td></tr><tr><td>T<sup>2</sup>/MML</td><td>1.00</td><td>0.85</td><td>0.93</td></tr><tr><td>% Section Overestimation</td><td>42.9%</td><td>25.7%</td><td>22.9%</td></tr></tbody></table>
<table><caption>**Table 2.** Values calculated for the uncertainty quantification for the C15 eruption.</caption><tr><td>% Section Underestimation</td><td>57.1%</td><td>74.3%</td><td>77.1%</td></tr><tr><td>MO (kg/m²)</td><td>0.47</td><td>0.09</td><td>0.25</td></tr><tr><td>MU (kg/m²)</td><td>-3.04</td><td>-2.66</td><td>-2.60</td></tr><tr><td>MO/MML</td><td>0.21</td><td>0.04</td><td>0.11</td></tr><tr><td>MU/MML</td><td>--1.37</td><td>-1.20</td><td>-1.18</td></tr></table>

1039

1040

<table><caption>**Table 3.** Values calculated for the uncertainty quantification for the T13 eruption.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="3">Meteo Data</th></tr><tr><th>GDAS</th><th>NCEP/NCAR</th><th>ERA-Interim</th></tr></thead><tbody><tr><td colspan="4">PLUME-MoM</td></tr><tr><td>Mean Difference (m)</td><td>-2202.05</td><td>-2113.75</td><td>-2132.03</td></tr><tr><td colspan="4">HYSPLIT</td></tr><tr><td>T²/MML</td><td>0.71</td><td>1.02</td><td>1.49</td></tr><tr><td>% Section Overestimation</td><td>29.2%</td><td>33.3%</td><td>35.4%</td></tr><tr><td>% Section Underestimation</td><td>70.8%</td><td>66.7%</td><td>64.6%</td></tr><tr><td>MO (kg/m²)</td><td>0.15</td><td>0.38</td><td>0.40</td></tr><tr><td>MU (kg/m²)</td><td>-0.47</td><td>-0.34</td><td>-0.37</td></tr><tr><td>MO/MML</td><td>0.28</td><td>0.75</td><td>0.77</td></tr><tr><td>MU/MML</td><td>-0.92</td><td>-0.67</td><td>-0.73</td></tr></tbody></table>

1041

1042

<table><caption>**Table 4.** Values calculated for the uncertainty quantification for the T06 eruption.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="4">Meteo Data</th></tr><tr><th>GDAS</th><th>NCEP/NCAR</th><th>ERA-Interim</th><th>ERA-Interim/WRF</th></tr></thead><tbody><tr><td colspan="5">PLUME-MoM</td></tr><tr><td>Mean Difference (m)</td><td>-718.67</td><td>-3752.59</td><td>-1225.72</td><td>404.74</td></tr><tr><td colspan="5">HYSPLIT</td></tr><tr><td>T²/MML</td><td>5.41</td><td>19.67</td><td>1.78</td><td>4.59</td></tr><tr><td>% Section Overestimation</td><td>39.5%</td><td>44.2%</td><td>27.9%</td><td>55.8%</td></tr><tr><td>% Section Underestimation</td><td>60.5%</td><td>55.8%</td><td>72.1%</td><td>44.2%</td></tr><tr><td>MO (kg/m²)</td><td>23.96</td><td>62.57</td><td>15.64</td><td>19.73</td></tr><tr><td>MU (kg/m²)</td><td>-3.32</td><td>-2.74</td><td>-4.40</td><td>-2.88</td></tr><tr><td>MO/MML</td><td>2.94</td><td>7.68</td><td>1.92</td><td>2.42</td></tr><tr><td>MU/MML</td><td>-0.41</td><td>-0.34</td><td>-0.54</td><td>-0.35</td></tr></tbody></table>

1043

1044

<table><caption>**Table 5.** Values calculated for the uncertainty quantification for the PCC11 eruption.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="6">Meteo Data</th></tr><tr><th>GDAS</th><th>NCEP/NCAR</th><th>ERA-Interim</th><th>GDAS (inversion)</th><th>NCEP/NCAR (inversion)</th><th>ERA-Interim (inversion)</th></tr></thead><tbody><tr><td colspan="7">PLUME-MoM</td></tr><tr><td>Mean Difference (m)</td><td>296.11</td><td>71.79</td><td>182.06</td><td>195.51</td><td>-84.86</td><td>33.27</td></tr><tr><td colspan="7">HYSPLIT</td></tr><tr><td>T²/MML</td><td>17.05</td><td>9.69</td><td>22.12</td><td>11.73</td><td>8.08</td><td>7.08</td></tr><tr><td>% Section Overestimation</td><td>50.7%</td><td>30.7%</td><td>34.7%</td><td>48.0%</td><td>32.0%</td><td>38.7%</td></tr><tr><td>% Section Underestimation</td><td>49.3%</td><td>69.3%</td><td>65.3%</td><td>52.0%</td><td>68.0%</td><td>61.3%</td></tr><tr><td>MO (kg/m²)</td><td>165.27</td><td>227.60</td><td>309.99</td><td>133.10</td><td>184.06</td><td>156.82</td></tr><tr><td>MU (kg/m²)</td><td>-53.93</td><td>-15.57</td><td>-31.94</td><td>-50.57</td><td>-18.72</td><td>-36.37</td></tr><tr><td>MO/MML</td><td>3.15</td><td>4.34</td><td>5.91</td><td>2.54</td><td>3.51</td><td>2.99</td></tr><tr><td>MU/MML</td><td>-1.03</td><td>-0.30</td><td>-0.61</td><td>-0.96</td><td>-0.36</td><td>-0.69</td></tr></tbody></table>

1045
<caption>Figure 1.</caption>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>
<caption>Figure 2.</caption>
<figure><img src="image_5.png" /></figure>

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

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

<figure><img src="image_8.png" /></figure>
<caption>Figure 3.</caption>
<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>
<caption>Figure 4.</caption>
a Plume Height

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

b Mass Loading

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

c Mdφ

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

d σφ

<figure><img src="image_16.png" /></figure>
<caption>Figure 5.</caption>
a

Plume Height

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

b

Mass Loading

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

c

Mdφ

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

d

σφ

<figure><img src="image_20.png" /></figure>
Figure 6.
a

Plume Height

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

b

Dispersal axis

Δ mass loading (kg/m²)

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

Mass Loading

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

c

Mdφ

Mdφ computed

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

d

σφ computed

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