M Alessandro Tadini

Olivier Roche

P Samaniego

A Guillin

N Azzaoui

M Gouhier

M De'michieli

F Pardini

Julia Eychenne

B Bernard

Alessandro Tadini
email: alessandro.tadini@uca.fr

M De'

Michieli Vitturi

S Hidalgo

J L Le Pennec

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

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

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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 ;[START_REF] Folch | A review of tephra transport and dispersal models: evolution, current status, and future perspectives[END_REF]. 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 [START_REF] Folch | FPLUME-1.0: An integral volcanic plume model accounting for ash aggregation[END_REF] or [START_REF] Costa | Results of the eruptive column model inter-comparison study[END_REF]; ii) for TTDMs, by comparing modelled and observed ground deposit measurements (mass loadings in kg/m 2 ) and/or ash cloud measurements (concentrations in the atmosphere in kg/m 3 ) [e.g., [START_REF] Scollo | A parametric and comparative study of different tephra fallout models[END_REF][START_REF] Costa | Tephra fallout hazard assessment at the Campi Flegrei caldera (Italy)[END_REF][START_REF] Bonasia | Numerical inversion and analysis of tephra fallout deposits from the 472 AD sub-Plinian eruption at Vesuvius (Italy) through a new best-fit procedure[END_REF][START_REF] Folch | A review of tephra transport and dispersal models: evolution, current status, and future perspectives[END_REF].

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 [START_REF] De'michieli | PLUME-MoM 1.0: A new integral model of volcanic plumes based on the method of moments[END_REF], and ii) the HYSPLIT model [START_REF] Stein | NOAA's HYSPLIT atmospheric transport and dispersion modeling system[END_REF], 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 [2015eruption, Bernard et al., 2016a] and Tungurahua [2006eruption, Eychenne et al., 2012;2013eruption, 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.

Background

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.

Cotopaxi 2015

The 2015 eruption of Cotopaxi (C15 -Fig. 1a) started with hydromagmatic explosions on August 14 th 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;[START_REF] Gaunt | Juvenile magma recognition and eruptive dynamics inferred from the analysis of ash time series: the 2015 reawakening of Cotopaxi volcano[END_REF].

The magmatic character of the eruption increased through time as was shown by microtextural analysis [START_REF] Gaunt | Juvenile magma recognition and eruptive dynamics inferred from the analysis of ash time series: the 2015 reawakening of Cotopaxi volcano[END_REF] and ash/gas geochemistry [START_REF] Hidalgo | Evolution of the 2015 Cotopaxi eruption revealed by combined geochemical and seismic observations[END_REF].

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 [START_REF] Gaunt | Juvenile magma recognition and eruptive dynamics inferred from the analysis of ash time series: the 2015 reawakening of Cotopaxi volcano[END_REF]. The hydrothermal components were dominant at the onset of the eruption but rapidly faded and were replaced by juvenile material [START_REF] Gaunt | Juvenile magma recognition and eruptive dynamics inferred from the analysis of ash time series: the 2015 reawakening of Cotopaxi volcano[END_REF]. In total, this eruption emitted ~1.2x10 9 kg of ash and was characterized as a VEI 1-2 [Bernard et al., 2016a].

Tungurahua 2013

According to [START_REF] Hidalgo | SO 2 degassing at Tungurahua volcano (Ecuador) between 2007 and 2013: Transition from continuous to episodic activity[END_REF], the eruptive phase XI (T13) at Tungurahua volcano (Fig. 1a) started on July 14 th 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 kmhigh eruptive column [START_REF] Parra | Eruption Source Parameters for forecasting ash dispersion and deposition from vulcanian eruptions at Tungurahua volcano: Insights from field data from the July 2013 eruption[END_REF]. 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 [START_REF] Parra | Eruption Source Parameters for forecasting ash dispersion and deposition from vulcanian eruptions at Tungurahua volcano: Insights from field data from the July 2013 eruption[END_REF]. The eruption intensity dropped after this paroxysm but ash emission continued with a secondary increase between July 20 th and 24 th . Finally the eruption stopped at the beginning of August.

In total, this eruption emitted ~6.7x10 8 kg of fallout deposits (~2.9x10 8 kg for the first day) and ~5x10 9 kg of pyroclastic flow deposits (mostly during the first day) [START_REF] García Moreno | 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[END_REF][START_REF] Parra | Eruption Source Parameters for forecasting ash dispersion and deposition from vulcanian eruptions at Tungurahua volcano: Insights from field data from the July 2013 eruption[END_REF]. [START_REF] Parra | Eruption Source Parameters for forecasting ash dispersion and deposition from vulcanian eruptions at Tungurahua volcano: Insights from field data from the July 2013 eruption[END_REF] performed numerical simulations of the vulcanian onset of this eruption, which occurred on July 14 th 2013, using the coupled WRF-FALL3D models [START_REF] Michalakes | Development of a next-generation regional weather research and forecast model[END_REF][START_REF] Folch | FALL3D: A computational model for transport and deposition of volcanic ash[END_REF]. 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.

Tungurahua 2006

At Tungurahua volcano (Fig. 1a), a paroxysmal eruption (T06) occurred on August 16 th 2006, which was accompanied by regional tephra fallout and many scoria flows and surges that devastated the western half of the edifice [START_REF] Douillet | Sedimentology and geomorphology of the deposits from the August 2006 pyroclastic density currents at Tungurahua volcano, Ecuador[END_REF][START_REF] Hall | Pyroclastic density currents (PDC) of the 16-17 August 2006 eruptions of Tungurahua volcano, Ecuador: Geophysical registry and characteristics[END_REF]. This eruption was characterized by vigorous lava jetting and fountaining, a vent-derived eruption column reaching 16-18 km above the vent [START_REF] Steffke | Eruption chronologies, plume heights and eruption styles at Tungurahua Volcano: Integrating remote sensing techniques and infrasound[END_REF]Eychenne et al., 2012], numerous Pyroclastic Density Currents (PDCs) descending the southern, western and northern flanks of the volcano [START_REF] Kelfoun | Testing the suitability of frictional behaviour for pyroclastic flow simulation by comparison with a wellconstrained eruption at Tungurahua volcano (Ecuador)[END_REF][START_REF] Bernard | Pyroclastic flow erosion and bulking processes: comparing field-based vs. modeling results at Tungurahua volcano, Ecuador[END_REF], and a massive blocky lava flow emplacing on the western flank while the explosive activity waned [START_REF] Samaniego | Petrological analysis of the pre-eruptive magmatic process prior to the 2006 explosive eruptions at Tungurahua volcano (Ecuador)[END_REF]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 [START_REF] Hall | Pyroclastic density currents (PDC) of the 16-17 August 2006 eruptions of Tungurahua volcano, Ecuador: Geophysical registry and characteristics[END_REF]. 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].

In total, the whole August 2006 eruption produced 39.3±5.1x10 6 m 3 of fallout deposit (both vent-derived and co-PDC derived) of which 24.9±3.3x10 9 kg were related to the ventderived fall [Bernard et al., 2016b].

Puyehue-Cordón Caulle 2011

According to [START_REF] Collini | Volcanic ash forecast during the June 2011 Cordón Caulle eruption[END_REF], the Puyehue-Cordón Caulle 2011 eruption (PCC11 -Fig. 1b) started on June 4 th at 14:45 LT (18:45 UTC) with the opening of a new vent 7 km 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 [START_REF] Jay | Locating magma reservoirs using InSAR and petrology before and during the 2011-2012 Cordón Caulle silicic eruption[END_REF] and comprised both explosive and effusive activity [START_REF] Tuffen | Exceptional mobility of an advancing rhyolitic obsidian flow at Cordón Caulle volcano in Chile[END_REF]. The main explosive phase, which dispersed most of the tephra toward E and SE, lasted approximately 17-27 hours [START_REF] Jay | Locating magma reservoirs using InSAR and petrology before and during the 2011-2012 Cordón Caulle silicic eruption[END_REF]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 14 th [Bonadonna et al., 2015a;[START_REF] Biondi | Supporting the detection and monitoring of volcanic clouds: A promising new application of Global Navigation Satellite System radio occultation[END_REF].

During the eruption, the mass eruption rate (MER) fluctuated between 2.8x10 7 (during the first days) and less than 5x10 5 kg/s after June 7 th [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 4 th -5 th . Unit I had a total erupted mass of 4.5±1.0×10 11 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]. [START_REF] Collini | Volcanic ash forecast during the June 2011 Cordón Caulle eruption[END_REF] performed numerical modellings of this eruption between June 4 th to June 20 th using the above-mentioned WRF-FALL3D code. The authors compared both the column mass load (in ton/km 2 ) 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 [START_REF] Marti | Volcanic ash modeling with the online NMMB-MONARCH-ASH v1. 0 model: model description, case simulation, and evaluation[END_REF], who simulated the eruption from June 4 th up to Jun 21 st using the NMMB-MONARCH-ASH model and compared the same parameters as in [START_REF] Collini | Volcanic ash forecast during the June 2011 Cordón Caulle eruption[END_REF]. 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 [START_REF] Pistolesi | Complex dynamics of smallmoderate volcanic events: the example of the 2011 rhyolitic Cordón Caulle eruption, Chile[END_REF], finding a good agreement between modeled and observed data.

Numerical modeling

Models used and coupling of the codes

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.

Following the approach adopted in [START_REF] Bursik | Effect of wind on the rise height of volcanic plumes[END_REF], 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 [START_REF] De'michieli | Uncertainty quantification and sensitivity analysis of volcanic columns models: Results from the integral model PLUME-MoM[END_REF] 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 height determination (i.e. the plume height decreased by a factor of ~1.54 when increasing water content from 1 to 5 wt%).

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. [START_REF] Webley | 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[END_REF] 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.

In the present study we coupled the PLUME-MoM and HYSPLIT models with an adhoc 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 bestfitting inverse version of this coupling, which was based on the approach first described by [START_REF] Connor | Inversion is the key to dispersion: understanding eruption dynamics by inverting tephra fallout[END_REF] and applied, among others, by [START_REF] Bonasia | Numerical inversion and analysis of tephra fallout deposits from the 472 AD sub-Plinian eruption at Vesuvius (Italy) through a new best-fit procedure[END_REF] 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 [START_REF] Wilson | The influence of shape on the atmospheric settling velocity of volcanic ash particles[END_REF][START_REF] Riley | Quantitative shape measurements of distal volcanic ash[END_REF].

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 2 function

𝑇 2 = ∑ 𝑤 𝑖 [𝑀𝐿 𝑜,𝑖 -𝑀𝐿 𝑚,𝑖 ] 2 𝑁 𝑖=1
where the sum is extended over N stratigraphic sections used in the inversion, w i are weighting factors (in our case all are equal to 1), ML o,i denotes the observed mass load (in kg/m 2 ) and ML m,i are the values predicted by the model (in kg/m 2 ). The values of T 2 will be then compared to the standard Chi-2 distribution of N-p degrees of freedom, with p=3 the number of free parameters.

Modelling features and input parameters

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.

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 bestfitting 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).

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 [START_REF] Andújar | Structure of the plumbing system at Tungurahua volcano, Ecuador: insights from phase equilibrium experiments on July-August 2006 eruption products[END_REF] 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% [START_REF] Martel | Storage conditions of the mafic and silicic magmas at Cotopaxi, Ecuador[END_REF]; 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/kgxK following [START_REF] Folch | FPLUME-1.0: An integral volcanic plume model accounting for ash aggregation[END_REF]; 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 [START_REF] Riley | Quantitative shape measurements of distal volcanic ash[END_REF]. For the inverse simulations of the PCC11 eruption, the particle shape factor values at each iteration were sampled between 0.6 and 0.8 [START_REF] Riley | Quantitative shape measurements of distal volcanic ash[END_REF]; h) the particle density values were assumed to vary linearly between two values (ρ 1 and ρ 2 ) specific of two grain sizes (φ 1 and φ 2 ) according to [START_REF] Bonadonna | Sedimentation from strong volcanic plumes[END_REF].

Values of ρ 1 , ρ 2, φ 1 , and φ 2 were taken from Eychenne and Le Pennec [2012] (C15/T13/T06) and [START_REF] Pistolesi | Complex dynamics of smallmoderate volcanic events: the example of the 2011 rhyolitic Cordón Caulle eruption, Chile[END_REF] (PCC11). For each eruption, all the other most relevant features of input parameters are described below.

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

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

For the Tungurahua T06 eruption, the simulations covered 4 hours corresponding to the climatic phases I and II described in [START_REF] Hall | Pyroclastic density currents (PDC) of the 16-17 August 2006 eruptions of Tungurahua volcano, Ecuador: Geophysical registry and characteristics[END_REF]. Plume heights were derived from [START_REF] Steffke | Eruption chronologies, plume heights and eruption styles at Tungurahua Volcano: Integrating remote sensing techniques and infrasound[END_REF]. 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).

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 [START_REF] Pistolesi | Complex dynamics of smallmoderate volcanic events: the example of the 2011 rhyolitic Cordón Caulle eruption, Chile[END_REF] 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).

Uncertainty quantification procedure

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.

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.

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φ and σφ [START_REF] Folk | Brazos River bar [Texas]; a study in the significance of grain size parameters[END_REF]. For mass loading we use hereafter the notation "Δ 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) Δ 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),

{ 𝑀𝑂 = ∑ ∆ 𝑖 𝑁 𝑜 𝑖=1 𝑁 𝑜 𝑓𝑜𝑟 ∆ 𝑖 > 0 𝑀𝑈 = ∑ ∆ 𝑖 𝑁 𝑢 𝑖=1 𝑁 𝑢 𝑓𝑜𝑟 ∆ 𝑖 < 0
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 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.

Results

For all the eruptions, Fig. 2 S4-S7) and the simulation outputs in PDF (Figures S1-S16).

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 [START_REF] Gaunt | Juvenile magma recognition and eruptive dynamics inferred from the analysis of ash time series: the 2015 reawakening of Cotopaxi volcano[END_REF] 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 seismicderived data, exceptions are the days around the 23 rd 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 and3b). Notice that the deposits simulated, despite in extremely low quantities (i.e. 10 -10 -10 -11 kg/m 2 ) 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 2 ; 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 2 for all the simulations). The T 2 /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.

Tungurahua 2013

For the T13 eruption, a total of 48 mass loading measurements [unpublished and from [START_REF] Parra | Eruption Source Parameters for forecasting ash dispersion and deposition from vulcanian eruptions at Tungurahua volcano: Insights from field data from the July 2013 eruption[END_REF] and 29 grain size analyses [unpublished and from [START_REF] Parra | Eruption Source Parameters for forecasting ash dispersion and deposition from vulcanian eruptions at Tungurahua volcano: Insights from field data from the July 2013 eruption[END_REF] 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 and4b).

The observed values of mass loading (Fig. 4b and Table S3 from the Supporting Information) are all <3 kg/m 2 , 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 2 /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).

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

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 and5b).

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 Mdφ 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).

Puyehue-Cordón Caulle 2011

For the PCC11 eruption, a total of 75 mass loading measurements and 24 grain size analyses [Bonadonna et al., 2015a;[START_REF] Pistolesi | Complex dynamics of smallmoderate volcanic events: the example of the 2011 rhyolitic Cordón Caulle eruption, Chile[END_REF]unpublished] were used for the comparison (Fig. 2d). For the mass loadings, the thickness data of [START_REF] Pistolesi | Complex dynamics of smallmoderate volcanic events: the example of the 2011 rhyolitic Cordón Caulle eruption, Chile[END_REF] were multiplied by the bulk deposit density value of 560 kg/m 3 reported in Bonadonna et al.

[2015a] for Unit I, in order to obtain kg/m 2 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 and6b).

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 2 and from -18 kg/m 2 up to -54 kg/m 2 . 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 [START_REF] Pistolesi | Complex dynamics of smallmoderate volcanic events: the example of the 2011 rhyolitic Cordón Caulle eruption, Chile[END_REF] 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 2 close to the vent (section n° 57, Table S3 from the Supporting Information) to 160 kg/m 2 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 2 , 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 2 ).

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.

Discussion

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. [START_REF] Devenish | Sensitivity analysis of dispersion modeling of volcanic ash from Eyjafjallajökull in May 2010[END_REF][START_REF] Webster | Operational prediction of ash concentrations in the distal volcanic cloud from the 2010 Eyjafjallajökull eruption[END_REF], 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., [START_REF] Webley | 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[END_REF][START_REF] Bonasia | Numerical simulation of tephra transport and deposition of the 1982 El Chichón eruption and implications for Confidential manuscript submitted to[END_REF][START_REF] Folch | A review of tephra transport and dispersal models: evolution, current status, and future perspectives[END_REF]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, [START_REF] Poret | Modelling tephra dispersal and ash aggregation: The 26th April 1979 eruption[END_REF]. 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 2 and the MO-MU values) for the lower magnitude C15 and T13 eruptions, 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 2 /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 [START_REF] Parra | Eruption Source Parameters for forecasting ash dispersion and deposition from vulcanian eruptions at Tungurahua volcano: Insights from field data from the July 2013 eruption[END_REF]. 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, [START_REF] Steffke | Eruption chronologies, plume heights and eruption styles at Tungurahua Volcano: Integrating remote sensing techniques and infrasound[END_REF] 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., [START_REF] Bonadonna | Sedimentation from strong volcanic plumes[END_REF]Eychenne and Le Pennec, 2012;[START_REF] Pistolesi | Complex dynamics of smallmoderate volcanic events: the example of the 2011 rhyolitic Cordón Caulle eruption, Chile[END_REF]. This is particularly important for the PCC11 eruption that has the highest T 2 /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 [START_REF] Clemens | Water contents of silicic to intermediate magmas[END_REF], among which the direct measurement from melt inclusions in crystals are the most used [see for example [START_REF] Plank | Why do mafic arc magmas contain 4 wt% water on average?[END_REF]. As a comparison, for this study we relied on estimates made both using direct measurements from melt inclusions and experimental approaches [START_REF] Martel | Storage conditions of the mafic and silicic magmas at Cotopaxi, Ecuador[END_REF] or considering only experimental approaches [START_REF] Andújar | Structure of the plumbing system at Tungurahua volcano, Ecuador: insights from phase equilibrium experiments on July-August 2006 eruption products[END_REF]: results gave H 2 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 [START_REF] De'michieli | Uncertainty quantification and sensitivity analysis of volcanic columns models: Results from the integral model PLUME-MoM[END_REF], its careful estimation is therefore of primary importance.

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 2 /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 [START_REF] Folch | A review of tephra transport and dispersal models: evolution, current status, and future perspectives[END_REF]. 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 [START_REF] Hurst | Forecasting volcanic ash deposition using HYSPLIT[END_REF]. 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 [START_REF] Bello | Effect of particle volume fraction on the settling velocity of volcanic ash particles: insights from joint experimental and numerical simulations[END_REF][START_REF] Gouhier | Low efficiency of large volcanic eruptions in transporting very fine ash into the atmosphere[END_REF]. However, the problem of the effect of diffusion on volcanic plumes dispersal and therefore on particle sedimentation is complex [see for example [START_REF] Devenish | Sensitivity analysis of dispersion modeling of volcanic ash from Eyjafjallajökull in May 2010[END_REF]: a more rigorous study is therefore needed for HYSPLIT to investigate the influence of different available diffusion equations on final results.

The failure to take into account such mechanisms implies that the simulated finestgrained 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 [START_REF] Collini | Volcanic ash forecast during the June 2011 Cordón Caulle eruption[END_REF], 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.

Regarding the simulated grain size data, the Mdφ values are systematically coarsergrained 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.

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.

Conclusions

This paper presents the coupling of the PLUME-MoM model with a renewed version of the HYSPLIT tephra dispersal model. These two coupled models have been tested against four eruptions of different magnitudes and styles from three Andean volcanoes. A procedure of uncertainty quantification has been applied by computing the differences between modeled and observed data of plume height, mass loading and grain size (in terms of Mdφ and σφ).  For the above-mentioned reasons, high amounts of fine particles might reduce the accuracy of the model when simulating mass loading and grain size data.

 If the PLUME-MoM/HYSPLIT model is meant to be employed for hazard assessment purposes, we recommend considering high magnitude eruptions (i.e.

sub-Plinian or Plinian) as target cases, and mass loading as primary parameter.

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

  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

Four

  different meteorological datasets (GDAS, NCAR/NCEP, ERA-Interim, ERA-Interim+WRF) have been tested as well. The main conclusions and future perspectives of this uncertainty quantification are: None of the meteorological datasets tested produced systematically the best results for all the eruptions. This implies that if a specific dataset is employed for numerical modelling, its uncertainty (as quantified here) should be considered. The PLUME-MoM model tends to underestimate measured plume heights, except for the eruption with the highest magnitude tested (i.e. PCC11). Though underestimation might be related to uncertainties in plume height measurements and input data (e.g. mass eruption rate and initial water mass fraction), we note that for most of the cases we investigated the mean underestimations of the model for each eruption were lower with respect to the uncertainty in observed data. The PLUME-MoM/HYSPLIT model tends generally to have more points underestimating the mass loading data. If absolute mean differences are considered, however, then overestimation is almost always higher than underestimation. The distribution of sections with overestimation and underestimation does not highlight systematically homogeneous areas of either overestimation or underestimation. The advective part of the HYSPLIT model appears to work more efficiently than the diffusive part. Moreover, the failure to take into account any collective settling mechanisms of volcanic ash in HYSPLIT might cause important discrepancies between observed and modeled data of mass loading and, above all, grain size distributions.

Figure 1 .

 1 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 (a) and UTM WGS84 19S (b) systems. Basemap copyright of ESRI®, DigitalGlobe, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN, and the GIS User Community.

Figure 2 .

 2 Figure 2. Stratigraphic sections used for uncertainty quantification (considering only mass loading or both mass loading and grain size) and dispersal axes from field data/simulations 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 Model (30 m-resolution) from Marc Souris, IRD (a-c) and ESRI, USGS, NOAA (d).Coordinates are in UTM WGS84 17S (a-c) and UTM WGS84 19S (d).

Figure 3 .

 3 Figure 3. Cotopaxi C15 eruption. Comparison of 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 4 .

 4 Figure 4. Tungurahua T13 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 5 .

 5 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 .

 6 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).

Figure 1 .

 1 Figure 1.

Figure 2

 2 Figure 2.

Table 1 .

 1 Main input parameters used for the simulations.

	TABLES
	Eruption name	Computational grid dimension (°)	Initial water content (wt%)	Particles exit velocity (m/s)	Heat capacity (J/kgxK)	Particles Shape Factor	φ 1	ρ 1 (kg/m 3 )	φ 2	ρ 2 (kg/m 3 )
	Cotopaxi
	2015	5x5	5.5%	135	1600	0.75	-1	1487	2	2478
	(C15)
	Tungurahua
	2013	6 x 6	5.5%	135	1600	0.75	-1	1487	2	2478
	(T13)
	Tungurahua
	2006	6 x 6	5.5%	275	1600	0.75	-1	1487	2	2478
	(T06)
	Puyehue-
	Cordón Caulle 2011	10 x 10	7.0%	275	1600	0.65	-4	500	5	2670
	(PCC11)
	Parameters	GDAS	Meteo Data NCEP/NCAR		ERA-Interim
				PLUME-MoM
	Mean Difference seismic (m)	-479.02	-443.75		-454.83
	Mean Difference video (m)	-466.24	-434.99		-443.23
	Mean Difference satellite (m)	-1449.64	-1405.96		-1366.64
				HYSPLIT
	T 2 /MML	1.00		0.85			0.93
	% Section Overestimation	42.9%	25.7%			22.9%

Table 2 .

 2 Values calculated for the uncertainty quantification for the C15 eruption.

	Parameters	GDAS	Meteo Data NCEP/NCAR	ERA-Interim
		PLUME-MoM
	Mean Difference (m)	-2202.05	-2113.75	-2132.03
		HYSPLIT
	T 2 /MML	0.71	1.02	1.49
	% Section Overestimation	29.2%	33.3%	35.4%
	% Section Underestimation	70.8%	66.7%	64.6%
	MO (kg/m 2 )	0.15	0.38	0.40
	MU (kg/m 2 )	-0.47	-0.34	-0.37
	MO/MML	0.28	0.75	0.77
	MU/MML	-0.92	-0.67	-0.73

Table 3 .

 3 Values calculated for the uncertainty quantification for the T13 eruption.

	Parameters	GDAS	Meteo Data NCEP/NCAR ERA-Interim ERA-Interim/WRF
			PLUME-MoM
	Mean Difference (m)	-718.67	-3752.59	-1225.72	404.74
			HYSPLIT
	T 2 /MML	5.41	19.67	1.78	4.59
	% Section Overestimation	39.5%	44.2%	27.9%	55.8%
	% Section Underestimation	60.5%	55.8%	72.1%	44.2%
	MO (kg/m 2 )	23.96	62.57	15.64	19.73
	MU (kg/m 2 )	-3.32	-2.74	-4.40	-2.88
	MO/MML	2.94	7.68	1.92	2.42
	MU/MML	-0.41	-0.34	-0.54	-0.35

Table 4 .

 4 Values calculated for the uncertainty quantification for the T06 eruption.

				Meteo Data
	Parameters	GDAS	NCEP/NCAR	ERA-Interim	GDAS (inversion)	NCEP/NCAR (inversion)	ERA-Interim (inversion)
				PLUME-MoM
	Mean Difference (m)	296.11	71.79	182.06	195.51	-84.86	33.27
				HYSPLIT
	T 2 /MML	17.05	9.69	22.12	11.73	8.08	7.08
	% Section Overestimation	50.7%	30.7%	34.7%	48.0%	32.0%	38.7%
	% Section Underestimation	49.3%	69.3%	65.3%	52.0%	68.0%	61.3%
	MO (kg/m 2 )	165.27	227.60	309.99	133.10	184.06	156.82
	MU (kg/m 2 )	-53.93	-15.57	-31.94	-50.57	-18.72	-36.37
	MO/MML	3.15	4.34	5.91	2.54	3.51	2.99
	MU/MML	-1.03	-0.30	-0.61	-0.96	-0.36	-0.69

Table 5 .

 5 Values calculated for the uncertainty quantification for the PCC11 eruption.

Acknowledgements

This research was financed by the French government IDEX-ISITE initiative 16-IDEX-0001 (CAP 20-25), the French Research Institute for Sustainable Development (IRD) in the context of the Laboratoire Mixte International "Séismes et Volcans dans les Andes du Nord" and the CNRS Tellus programme. Marco Pistolesi is warmly acknowledged for providing additional unpublished data for the Puyehue-Cordón Caulle 2011 eruption. The authors are grateful to Sandra Banson and Cèline Planche for the help in setting up the simulations with WRF. We thank two anonymous reviewers for useful comments that improved the quality of the manuscript, and we acknowledge as well the editorial handling provided by Nico Fournier and Martha Savage.