The Impact of Eruption Source Parameter Uncertainties on Ash Dispersion Forecasts During Explosive Volcanic Eruptions

Volcanic ash in the atmosphere is a hazard to aviation. To predict which areas of airspace are most likely to be affected by the presence of ash, Volcanic Ash Advisory Centers (VAACs) use observations and atmospheric dispersion models. These models are initialized with, among other parameters, a mass eruption rate (MER), which quantifies the emission rate into the atmosphere at the source. This influences the predicted spatial–temporal evolution and concentration of the ash cloud. Different models are available to estimate MER from the volcanic plume height and some models also include the weather conditions (e.g., wind speed). The REFIR software tool uses time‐series of plume height estimated from observations and weather data to provide estimates of MER through time. Here we present an updated version of REFIR that can now be used also to calculate MER for past eruptions and produce output parameters in a format suitable for use with the NAME dispersion model (UK Met Office—London VAAC). We also investigate how uncertainty in input parameters is propagated through to dispersion model output. Our results show that a +/−1 km uncertainty on a 6 km high plume can result in the affected area ranging by a factor of three between the minimum and maximum estimates. Additionally, we show that using wind‐affected plume models results in affected areas that are five times larger than using no‐wind‐affected models. This demonstrates the sensitivity of MER to the type of plume model chosen (no‐wind‐ vs. wind‐affected).


Introduction
Explosive volcanic eruptions can release huge volumes of ash into the atmosphere. As volcanic particles can damage aircraft, the resulting ash clouds can pose a serious threat to civil aviation (Durant et al., 2010;Giehl et al., 2017). When particles settle on the ground, they can also impact the environment, buildings, and infrastructure and if resuspended can affect human health (Blake et al., 2017;Blong et al., 2017;Jenkins et al., 2015). Current mitigation measures for aviation safety rely on exclusion zones for air traffic, which can lead to airport closures and flight diversions or cancellations; the consequent travel disruption can potentially affect many countries for prolonged periods (Guffanti & Tupper, 2015;Guffanti et al., 2010;Lechner et al., 2017). The 2010 eruption of the Eyjafjallajökull volcano demonstrated that even a small-to medium-sized eruption can have a big economic impact on aviation when the event is long-lived and where the weather conditions transport the ash into busy airspace (Harris et al., 2012).
When considering whether it is safe for aircraft to fly, local civil aviation authorities, airlines, and airports use forecasts of the movement of volcanic ash clouds issued by Volcanic Ash Advisory Centers (VAACs), referred to as a Volcanic Ash Graphic (VAG) and Advisory (VAA). These forecasts are generated using observations of the ash cloud (e.g., satellite imagery and lidar retrievals) and atmospheric dispersion models to predict the expected onwards transport of the ash cloud.
Our current capability to predict the spatial-temporal evolution of volcanic ash clouds in the atmosphere depends on (1) the characterization of the eruptive source via measurable and/or inferable physical quantities, the Eruption Source Parameters (ESPs); (2) observation and modeling of weather conditions; and (3) the application of dispersion models. ESPs include particles' properties (size, density, shape), the emission rate usually referred to as the Mass Eruption Rate (MER) in kg s −1 (Wilson & Walker, 1987) and some geometric properties of the emission (e.g., the volcanic plume height) (e.g., Dellino et al., 2014;Sparks, 1997;Woods, 1988). Weather conditions play a fundamental role in controlling the evolution of the eruptive plume and the subsequent advection of the emitted volcanic ash: for example, vertical profiles of temperature and humidity above the volcano can enhance or inhibit convection in the volcanic plume, while wind can bend the eruptive column and control the direction of ash advection (e.g., Aubry et al., 2017;Bursik, 2001;Costa et al., 2016;Degruyter & Bonadonna, 2012, 2013Devenish, 2013;Hewett et al., 1971;Macedonio et al., 2016;Tupper et al., 2009;Woodhouse et al., 2013Woodhouse et al., , 2015. The subsequent transport, dispersion, and removal (wet and dry deposition including sedimentation) of ash in the atmosphere can be predicted by applying atmospheric dispersion models (e.g., Barsotti et al., 2008;Beckett et al., 2015;Costa et al., 2006;Folch, 2012;Folch et al., 2009;Jones et al., 2007;Mastin et al., 2013;Stein et al., 2015). These models approximate the source of the emission to a simple shape (e.g., point, line, top-hat). In an operational environment, a vertical line source is often used, where model particles, each representing a proportion of the mass erupted, are released with a uniform distribution between the volcano summit and the plume top height (e.g., Costa et al., 2006;Witham et al., 2019;WMO, 2012). Sensitivity studies have shown that forecasts are particularly sensitive to both the plume height and MER used to initialize dispersion models (Beckett et al., 2015;Harvey et al., 2018;Scollo et al., 2019), although particle properties like size and shape play also a major role (Beckett et al., 2015;Saxby et al., 2020).
The location and evolution of an ash cloud and the plume top height can be monitored by means of different remote sensing techniques, which can be classified as ground-based and satellite-based. Ground-based methodologies include visible (e.g., Corradini et al., 2016;Scollo et al., 2014) and infrared (e.g., Valade et al., 2014) cameras and different types of radars (e.g., Arason et al., 2011;Donnadieu et al., 2016;Marzano et al., 2013). Satellite methodologies include sensors measuring electromagnetic radiation using different spectral channels (visible and infrared), with specific bands being absorbed by volcanic ash of specific size ranges and/or volcanic gas species like SO 2 (e.g., Corradini et al., 2018;Prata, 1989;Prata & Kerkmann, 2007) thus identifying its presence in the atmosphere. Since different sensors can acquire data related to specific gas species and volcanic particle sizes and are affected by different degrees of uncertainty and limitations, research has recently focused effort on combining data obtained from different techniques so as to better approximate the whole range of particle sizes emitted and to better assess, for example, the plume top height (e.g., Corradini et al., 2016;Marzano et al., 2018;Poret et al., 2018).
To date, MER cannot be directly estimated by remote sensing techniques. Some progress has been made recently, for example, by using infrasound waves (e.g., Ripepe et al., 2013), combining thermal infrared signatures with inverse modeling (Cerminara et al., 2015) or electrostatic field signals (Calvari et al., 2012) and analyzing the spreading of the volcanic cloud from satellite imagery in certain conditions (Pouget et al., 2013(Pouget et al., , 2016. The most widely used approach for estimating MER consists of using the plume height and weather conditions to infer relevant flow parameters at the volcanic source (e.g., MER) by means of physical models of volcanic plumes. Several models are available in the volcanological literature; a comprehensive review of these can be found in the intercomparison study published by Costa et al. (2016). Volcanic plume models range from simple 0-D to complex 3-D unsteady multiphase models. 0-D models can be either empirical equations linking MER to the observed plume height, which are based on data from past eruptions (e.g., Gudmundsson et al., 2012;Sparks, 1997) or theoretical models developed from the original analysis of Morton et al. (1956) on buoyant plumes (e.g., Wilson & Walker, 1987). The analysis of Morton et al. (1956) has been subsequently revised to consider the effect of a crossflow (Hewett et al., 1971); the latter description has been used to derive wind-affected 0-D volcanic plume models (e.g., Degruyter & Bonadonna, 2012;Woodhouse et al., 2013). The 1-D models, instead, solve for the conservation equations of mass, momentum, and energy along the plume propagation direction (centerline). They can be either steady (e.g., Bursik, 2001;Degruyter & Bonadonna, 2013;de Michieli Vitturi et al., 2015;Devenish, 2013;Folch et al., 2016;Mastin, 2014;Woodhouse et al., 2013) or time-dependent (e.g., Scase & Hewitt, 2012;Woodhouse et al., 2016). In steady 1-D models, important transient processes like turbulent mixing are modeled via an entrainment parametrization. Turbulent entrainment of external air into the volcanic plume is a fundamental process controlling the evolution of the volcanic plume; it significantly contributes to buoyancy, since the entrained air reduces the bulk density of the volcanic jet near the sources and, once in contact with the still hot volcanic gases and ash, expands thus further reducing the flow bulk density. The 1-D plume models that take into account the wind distinguish between two contributions: a radial entrainment, which acts on the plume margin from all directions when the plume is vertically directed, and a wind entrainment, which accounts for the enhancing effect of wind in this process where the plume is bent over. The final category of models, 3-D unsteady models, can capture the fluid dynamics of the volcanic plume, in some cases solving for flow fields of the different phases (gas and particles) composing the eruptive mixture (e.g., Cerminara et al., 2016;Esposti Ongaro et al., 2007).
Due to their high computational demand, to date 3-D models are not suitable for operational applications, with the 0-D and 1-D models being preferred instead. Ideally, the 0-D and 1-D models should be used in real-time to quantify MER and to inform ash dispersion forecasters during an explosive volcanic eruption. This is a complex task, which requires the combination of ash plume top height data from different monitoring sources and the related uncertainty. The uncertainty of ash plume top height propagates into the estimation of MER through 0-D and 1-D models and, finally, into the dispersion model outputs. In order to simplify this task, a quasi-autonomous real-time multi-parameter system called REFIR (Real-time Eruption source parameters FutureVolc Information and Reconnaissance system) was developed within the EU-Funded Horizon 2020 Supersite project FUTUREVOLC (futurvolc.hi.is; Dürig et al., 2018). The system was designed to combine real-time data streams of plume top height from different sensors, including C-Band and X-Band radars, auto-tracking web cameras, imaging ultra-violet and infrared cameras, and electric field sensors and consider the best estimate of ash plume top height and the uncertainty on this data. The best estimate of the ash plume top height was then used to calculate MER in near real-time by using a set of plume models.
Over the last year, REFIR has been substantially improved to enable more flexibility and its application to study past eruptions in order to obtain time series of ESPs. In this manuscript, we present the new capabilities of REFIR and results of their application to the 2010 Eyjafjallajökull eruption Dioguardi et al., 2016;Dürig, Gudmundsson, Karmann, et al., 2015;Gudmundsson et al., 2012). In particular,  estimated MER in the period 8-10 May by applying video analysis of footages of the erupted ash cloud from the crater. In this research, we analyzed the period 5-8 May of the 2010 Eyjafjallajökull eruption and used a different methodology to quantify MER; our results are compared to those of  for the evening of 8 May.
We use plume height data and MER time-series generated by REFIR to initialize the atmospheric dispersion model NAME (Numerical Atmospheric-dispersion Modeling Environment; Jones et al., 2007), which is used by the London VAAC as part of its operational forecasting process. We consider how the uncertainty assessment of plume height and the choice of different combinations of 0-D plume models affect the final estimation of MER by REFIR and its uncertainty and how this, in turn, impacts the dispersion model forecasts. All the above shows the potential and flexibility of REFIR in assessing the ESPs in real-time (as an operational tool) and for past eruption (as a research tool).

The REFIR System and the Newly Introduced Capabilities in v19.0
In this section, we briefly introduce REFIR, referring to the original open-access paper of Dürig et al. (2018) and the User Manual included in the REFIRv19.0 package, which is made available with this manuscript as supporting information. Here we focus on the new capabilities recently introduced, which allowed us to conduct this study.
The current version of REFIR (v.19.0) is a collection of Python scripts (v.3 and above) of which the two main scripts are FIX.py and FOXI.py. The former is the main control unit; the user can specify sensors data to be used, model parameters, and the type of outputs via this script. The latter is the processing unit, that is, the script that performs the ESP calculations and writes the results. Input parameters and specification of the outputs are set through a graphical user interface (GUI).
REFIR was originally designed to calculate MER using five different 0-D plume models (Degruyter & Bonadonna, 2012;Gudmundsson et al., 2012;Sparks, 1997;Wilson & Walker, 1987), with the option to remotely activate the Plumerise 1-D model of Woodhouse et al. (2013). The process of plume height data acquisition is repeated every 5 min; the user can then choose a time window (known in REFIR as "time base") over which plume height data are processed for assessing the best estimate and range of variation in the plume top height (15, 30, 60, or 180 min). A small time base considers only the most up to date data, but reducing the amount of processed data increases the uncertainty. Whereas, a longer time base considers a bigger dataset and hence not only reduces the uncertainty but also reduces the time resolution.
Recently, REFIR v18.0, the last version made available in Dürig et al. (2018), has gone through major modifications in order to improve its range of possible applications and flexibility. These are explained in detail in the following.

REFIR Operation Mode
REFIR was originally developed to work in real-time during volcanic eruptions. However, it has the potential to be used for research applications as well, that is, retrieving time series of ESPs of past eruptions. We have implemented a new mode of operation, the reanalysis mode, which allows the user to obtain time series of ESPs of past or simulated explosive eruptions. Further details on how to select the desired operation mode can be found in the User Manual available in the software package.

Automatic Weather Data Retrieval
In v18.0 of REFIR the user had to specify weather conditions above the volcanic vent, specifically the temperature and pressure at the volcanic source, the temperature gradients above the vent at specific height intervals and wind speed at the tropopause.
These data were needed to apply the model of Degruyter and Bonadonna (2012), which requires the atmospheric density at the volcanic vent ρ a0 , the average buoyancy frequency N and wind speed V from the source to the plume top height (see Table 1 for notation).
The average buoyancy frequency N was calculated by the temperature profile computed by FIX.py based on the user inputs, while average wind speed was approximated to be equal to the wind speed at the tropopause. In order to ensure that more realistic weather data are used both in reanalysis and real-time mode, in v19.0, an automatic weather data retrieval has been implemented and is activated by default; the user can still deactivate this new functionality and provide weather data manually (see User Manual). If this option is active, FIX.py calculates the necessary weather parameters from data retrieved from forecasts or reanalysis archives. The data retrieval and processing is carried out by new scripts stored in the folder "weather" in the REFIR package. Further details can be found in the User Manual.
It must be noted that ECMWF has recently announced that ERA-Interim will not be updated from the 1 September 2019. For this reason, future versions of REFIR will make use of the ERA5 dataset instead (https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5), which will be kept up to date by ECMWF and has a higher spatial (~30 km grid) and temporal (1 h) resolution. In this study, we used ERA-Interim data.
Atmospheric temperature at volcanic vent K T m0 Eruptive mixture temperature at volcanic vent Once data have been downloaded, the weather package extracts the relevant data over the vertical coordinate above the volcanic vent and saves these data in a profile file, an example of which can be found in the User Manual. Finally, it calculates the weather parameters needed by FOXI.py to calculate MER with the wind-affected plume models (atmospheric pressure and temperature at the vent location; plume height-averaged buoyancy frequency and wind speed; the g W s parameter of the Woodhouse et al., 2013, model [see below]), which are stored in separate files.

Implementation of Woodhouse et al. (2013) 0-D Model
The 0-D model of Woodhouse et al. (2013) (hereafter referred to as "Woodhouse0D") has been implemented in REFIR v19.0. This is presented in Equation 1 : where g W s is a parameter quantifying the strength of the wind shear from the ground to a reference height H 1 (e.g., the plume top height H 1 = H): where V 1 is the wind speed at the reference height H 1 and N is the average buoyancy frequency, defined as: where g is the gravitational acceleration, C is the heat capacity, T is the temperature, z is the vertical coordinate above the source and the subscripts a and 0 refer to the atmosphere and the volcanic source vent height, respectively (Table 1). Note that plume heights H are always above the ground (volcano summit) in this manuscript unless otherwise stated.
With the implementation of the Woodhouse 0-D model, the number of internal MER models included in REFIR v19.0 increased to six. The other five models are listed below: a. "Wilson & Walker," a theoretical model by Wilson and Walker (1987): where c is a constant which is calibrated to be 236 m(s kg −1 ) 1/4 . b. "Sparks," an empirical model by Sparks (1997): where ρ DRE is the Dense Rock Equivalent density of the erupted material and c is a constant calibrated to be 1,670 m(s m −3 ) 1/3.86 .
c. "Mastin," an empirical model by : with ρ DRE being fixed to 2,500 kg m −3 and c is a constant calibrated to be 2,000 m(s m −3 ) 1/4.15 .

10.1029/2020JD032717
Journal of Geophysical Research: Atmospheres d. "Gudmundsson," an empirical model by Gudmundsson et al., 2012, which has been obtained by calibrating the Mastin model (Equation 6) to the 2010 Eyjafjallajökull eruption data: where H avg and H max are the average and maximum estimate of top plume height, hence taking into account the uncertainty of H determination (e.g., with ground based radar), c equals the Mastin's constant, a is a dimensionless constant calibrated to be 0.0564, k 1 is a scaling factor which was found to be 2.15 for the first stage (14-16 April 2010) of the Eyjafjallajökull eruption, then dropping to 1.58 for the rest of the eruption. The default value used by REFIR is 1.6, which is the value used in the application described in this manuscript.
e. "Degruyter & Bonadonna," a theoretical model based on an algebraic relationship calibrated applying a 1-D numerical model, based on the combination of Morton et al. (1956) and Hewett et al. (1971). It is designed for wind-affected bent-over plumes: where V is the plume height-averaged wind speed, α and β are the radial and wind entrainment coefficients, respectively, z 1 = 2.8 is the maximum nondimensional height determined by numerical integration of the nondimensional governing equations of Morton et al. (1956) and g′ is the reduced gravity at the source defined as: where the subscript m refers to the eruptive mixture.
The six models contribute to the calculation of the so-called REFIR Conventional models MER (CMER, Dürig et al., 2018), and the user can select which model to use by means of weight factors w ranging from 0 to 1. CMER, combined with other sources of MER estimates (see User Manual for further explanations), results in the Final best estimate of MER (FMER as in Dürig et al., 2018). It should be noted that Dürig et al. (2018) suggested some combinations of weight factors suitable for different eruption styles, though these recommendations did not include Woodhouse0D model yet. In this work, we exploited this capability of selecting-deselecting specific models as to investigate the impact of this choice on the REFIR outputs and, subsequently, on the ash dispersion simulations.
By default, REFIR combines MER obtained with the six internal 0-D models (Equations 1 and 4-8) to obtain the best estimate of MER and related uncertainty (Dürig et al., 2018) (Equations 10a-10c). However, the user can manually switch off one or more models by acting on the weight factors (see User Manual). The user can specify any value ranging from 0 to 1: A weight factor of 0 for a specific model would ignore MER obtained with that model; vice versa, a weight factor >0 would include that model in the final estimate of MER.
where i is the index identifying the MER model, ranging from 1 to 6, and min, max, and avg are the subscripts referring to the minimum, maximum and average estimate.

Database of Volcanoes
FoxSet.py is a script located in the refir_config folder of the REFIR package that can be used to, for example, include a new volcano in the REFIR simulation. In v18.0, the user had to manually specify the coordinates and the elevation of the volcano of the volcanic source. In the new version, when adding a new volcano the user should just type the Smithsonian Institute ID (Global Volcanism Program, 2013); FoxSet.py then interrogates the dataset now included in refir_config (SI_volcanoes_list.xlsx) to obtain the coordinates and summit elevation of the selected volcano automatically.

Time-Averaged Outputs
Since REFIR was originally designed to analyze the radar data streamed by the Icelandic Met Office at a time resolution of 5 min, by default REFIR v18.0 produced output files which were updated every 5 min. We have introduced new functionality in v19.0 which allows the user to write time-averaged results over a time interval that can be selected from five options: 15, 30, 60, 180, and 360 min. This new feature has been designed to make REFIR outputs easier to use in ash dispersion modeling applications, in which the volcanic source is generally modeled at a time resolution significantly longer than 5 min. Additionally, the user can also activate the option to save the ESP time series in a format that is directly readable by the NAME model (Jones et al., 2007).

The Eyjafjallajökull 2010 Eruption
We have focused the application of REFIR v19.0 on the April and May 2010 Eyjafjallajökull eruption, since this eruption showed the crucial importance of timely and reliable estimates of ESPs for accurate ash dispersion forecasts. Additionally, a time series of plume height data retrieved from the C-band radar based in Keflavík airport has been made available by Arason et al. (2011).
Eyjafjallajökull is an ice-capped central volcano in southern Iceland, rising to 1,666 m above sea level (Jonsson, 1988;Saemundsson, 1979). The eruption was characterized by different phases, described in detail by Dellino et al. (2012) and Gudmundsson et al. (2012). A flank eruption (20 March-12 April) was characterized by lava fountaining and flows from the ice-free area of Fimmvörðuháls between the Eyjafjallajökull and Mýrdalsjökull ice caps; the first explosive phase at the summit (14-18 April) was characterized by a plume rising to 5-10 km high alternating between dark gray (loaded with tephra) and some steam-rich plumes. The second phase (18 April-4 May) was characterized by mainly lava effusion and much weaker explosive activity. The third phase (5-17 May) was characterized by increased explosive intensity before a final phase of decline (18)(19)(20)(21)(22). During most of the eruption, the volcanic plume was generally bent-over by the wind, a situation known in volcanology as a "weak plume." The ash plume reached a maximum height of 10 km a.s.l. in the first explosive stage. In the three-day period analyzed in this research (5-8 May 2010), the plume height increased up to~10 km a.s.l. on the 5 May, then decreased to an almost stable value of 6 km a.s.l.

Ground-Based Radar Retrieval of Plume Top Height
A time series of ash plume top height data was recorded by the C-band ground based weather radar installed at Keflavik Airport (Arason et al., 2011). The time series has a time resolution of 5 min, covers the whole eruption duration and has some data gaps when the plume was not detected by the radar (e.g., below the minimum visible elevation, which in this specific case is 2.9 km a.s.l., Arason et al., 2011). For our application, we selected the period 5-8 May 2018, since this period is affected by few data gaps, which have been filled with plume height data obtained from linear interpolation. Furthermore, video data of the vent taken on 8 May allowed inference of the MER with an alternative approach, which is independent of plume height models Dürig, Gudmundsson, Karmann, et al., 2015). Figure 1a shows the time series used for the period analyzed.
Black dots represent the plume heights; the error bars have been calculated by REFIR. By default, REFIR applies the following equation of Arason et al. (2011) to obtain the absolute uncertainty of plume height ΔH linked to ground based radars:

10.1029/2020JD032717
Journal of Geophysical Research: Atmospheres where d is the distance between the volcano and the radar and δ is the radar beam width. REFIR uses information about the radar type (C-or X-band, which, together with the geometry of the antenna, controls the beam width) and location available from the FUTUREVOLC project and constantly updated on this website: http://brunnur.vedur.is/radar/vespa/ In this case, d ≈ 130 km and δ = 0.9°resulting in ΔH ≈ 1.2 km.
The plume height time series produced by REFIR with associated uncertainties calculated using Equation 11 and with a time resolution of 5 min are shown in Figure 1b. Note that the original graphical output of REFIR has been slightly modified for a better representation in this manuscript. In this specific case, in which only one plume height data series is used, the uncertainty coincides with that calculated in Equation 10 for the radar based in Keflavik. Consequently, the curves of H max and H min are symmetrically positioned around the average plume height H avg . Theoretically, REFIR can combine both plume height data and related uncertainties when more than one data series are available (Dürig et al., 2018).

Weather Conditions
For the time interval under analysis, REFIR automatically retrieved ERA-Interim data. Since ERA-Interim analyses are available every 6 h and the REFIR internal time-step is 5 min, REFIR uses the latest available

Journal of Geophysical Research: Atmospheres
weather data at each time step. For example, the first available ERA-Interim analysis for the selected time interval is at 00:00 UTC on 5 May 2010. This analysis has been used by REFIR to compute the necessary weather parameters in the time interval 00:00-05:55 UTC. The weather data have been interrogated to obtain vertical profiles above the eruptive vent by interpolating to the volcanic vent location and, from these, to obtain the weather parameters needed to run Woodhouse0D and Degruyter and Bonadonna models: N (Equation 3), V, the wind speed at plume top height V(H) and g W s (Equation 2). All these parameters depend on H; hence, Figure 2 shows the time series of N, V, V(H), and g W s over the analyzed interval for H avg , H max , and H min . Note that all parameters change significantly with the weather conditions, as can be inferred by the abrupt changes every 6 h, which is the time resolution of the available weather data. However, the average buoyancy frequency N looks less sensitive. The smaller scale variations reflect changing of plume height, hence showing that N, unlike the other parameters, is more sensitive to plume height changes than weather conditions themselves.
From 5 to 9 May 2010, the weather conditions were characterized by average wind speeds >10 m s −1 , with peak velocities >20 m s −1 . Moderately windy conditions, together with a weak to moderate explosive activity, made the plume almost continuously weak and bent-over.

Dispersion Modeling
We use NAME v7.2, the Met Office's Lagrangian atmospheric dispersion model (Jones et al., 2007), to simulate the transport and dispersion of ash emitted from the Eyjafjallajökull volcano between 00:00 UTC on 5 May and 23:30 UTC on 8 May 2010. Model particles are advected by 3-D wind fields; for this study, we use analysis data from the Global configuration of the Met Office's Unified Model (UM; Davies et al., 2005) which, for the time period considered, has a horizontal resolution of~25 km (at mid-latitudes). Small-scale three-dimensional atmospheric turbulence and unresolved mesoscale motions are parameterized within NAME using random-walk techniques . NAME also includes parameterizations of sedimentation, dry deposition, and wet deposition that simulate the removal of volcanic ash from the atmosphere (Webster & Thomson, 2011).
To initialize NAME, we apply the operational set-up used by the London VAAC to represent the eruption source. Model particles, each representing a mass of volcanic ash, are released as a vertical line source with a uniform distribution, between the vent and the plume top height. Particles are assumed to be spherical with a density of 2,300 kg m −3 , and the particle size distribution (PSD) represents particles with diameter between 0.1 and 100 μm (Witham et al., 2019). The plume height and source strength (MER) are taken from the time-averaged REFIR outputs over a time interval of 30 min exploiting the new capabilities of REFIR (see section 2.5). As we are considering only the long-range transport of the smallest particles (with diameters ≤100 μm), the total MER calculated by REFIR (which consider all the emitted material) is scaled to account for the near-source fallout. We use just 5% of the total MER (Dacre et al., 2011;Devenish et al., 2012), We present 6 h averaged peak ash concentration charts, as generated by the Met Office during operational response to supplement the London VAAC graphics (VAG) and advisories (VAA) (ICAO, 2016). As the vertical structure of the distal ash plume is difficult to resolve, due to uncertainty in the vertical distribution of ash at the source and limitations in our ability to model atmospheric turbulence (Devenish et al., 2012;Webster et al., 2012), a "thin-to-thick" conversion is applied to output air concentrations. The peak concentration within thin layers having a thickness of 25FL (where FL is Flight Level in hundreds of feet) across prescribed thick layer depths (FL000-FL200, FL200-350, FL350-550) is taken and applied to the whole thick layer (Witham et al., 2019). This does not imply that the maximum concentration applies to the whole layer, just that within the layer we expect the peak concentration.

Results
In this section, we present the time series of MER obtained with REFIR and results of their usage, together with H, as source conditions of ash dispersion modeling with NAME.
We exploited the capability of REFIR to activate/deactivate selected models by acting on their weight factors (setting these to 1 when the model is active and 0 when not active) to quantify how the choice of a specific

10.1029/2020JD032717
Journal of Geophysical Research: Atmospheres model or a type of models affects the final estimate of MER. In the suite of 0-D models available in REFIR, it is possible to distinguish between "not wind-affected" models (Wilson & Walker, Sparks, Mastin, Gudmundsson) and "wind-affected" models (Degruyter & Bonadonna,Woodhouse0D), hereafter referred to as NWA and WA, respectively. We therefore distinguished between four cases of MER estimation: 1. with all models (hereafter referred to as ALL); 2. with the Mastin model only (hereafter referred to as MASTIN). This case has been selected since the Mastin model is the default model used operationally by London VAAC (Witham et al., 2019); 3. with WA models only (Degruyter & Bonadonna,Woodhouse0D); 4. with NWA models only (Wilson & Walker, Sparks, Mastin). Note that we excluded Gudmundsson

REFIR Assessed Time Series of Eruption Source Parameters
The time series of MER min , MER avg , and MER max obtained in the four cases listed above are shown in Figure 3: For a few minutes in the time series, H min = H − ΔH fell below the volcano's summit elevation. This is due to a very low plume height detected at that time; in this case, the value of H min above the vent, which is the one used in the plume models, would be negative, causing REFIR to declare an error and stop the iteration. To deal with these situations and complete the simulation, the current version of REFIR sets H min = 1 when it is found to be negative. This very low value, even if it allows REFIR to continue its computations, produces values of MER ≪ 1 kg s −1 for the empirical not wind-affected models, which is the reason for the very low values observed in the time series.
A visual inspection of Figure 3 reveals some important aspects: 1. WA models produce higher values of MER compared to NWA. 2. The case MASTIN produces the lowest values of MER among the four tested cases. MASTIN is part of the NWA set of models (Wilson & Walker, Sparks, Mastin), which on average produce higher estimates of MER than MASTIN. In fact, the other NWA models (Wilson & Walker, Sparks) predict higher MER taken individually for H < 22 km. For H > 22 km, Sparks' model predicts lower MER than Mastin's one, but this is not the scenario under analysis. This can easily be inferred by looking at Equations 4-6; even if the Mastin relationship is characterized by the highest exponent of the power law, the constant dividing H has the highest value among all the NWA models. 3. The time series of MER closely follow the plume height data, without evidence of jumps when the weather conditions changed (compare Figure 3 with Figures 1 and 2). This suggests that, while realistic weather data are preferable when using wind affected models, plume height is the most important parameter for the plume models considered in this study. 4. The difference between MER max and MER avg is greater than that between MER min and MER avg , with the former being lower in the WA case than for NWA (see Table 2).
These differences have been quantified by means of the formula in Equations 12a-12e and the results are listed in Table 2: where the index j refers to a time step of the time series and N is the total number of time steps.
When comparing the single solutions of MER between the WA and NWA cases (Table 2), the differences are much higher, ranging from 320% for MER max; WA − MER max; NWA (the maximum solution) to 331% for MER avg; WA − MER avg; NWA up to 478% for MER min; WA − MER min; NWA (the minimum solution). These differences are highly significant and can be explained by the different formulations employed in WA and NWA models. NWA models are generally simple power laws of H, with an exponent~4; for this reason, very low plume heights result in very low values of MER, as previously discussed, and the difference between MER max and MER avg is larger than the difference between MER min and MER avg (at constant ΔH) ( Table 2). WA models are more complex functions of H and other meteorological parameters, but both show the tendency to reduce the difference between MER max and MER avg . The higher values of MER in the WA case are due to the fact that the weather conditions were characterized by a moderate wind field throughout the analyzed period. In this case, WA models compute higher values of MER compared to NWA models, and this tendency would be the opposite in the case of no wind (V ¼ 0).
To prove this, in Figure 4 we display curves of MER vs. H for the five considered models (Equations 1, 4, 5, 6, and 8) in the case of no wind ( Figure 4a) and windy conditions (Figure 4b). For simplicity, in the latter case a uniform average wind speed for the whole plume

Journal of Geophysical Research: Atmospheres
height of 20 m s −1 is assumed. A uniform value of the average buoyancy frequency N ¼ 0:01 is also assumed for both scenarios.
In the case of no-wind (Figure 4a), NWA models produce higher MER values than WA models, with the latter computing very similar MER values. Interestingly, the curves of the two empirical models (Mastin and Sparks) fall between the Wilson & Walker curve and the WA models' curves, although closer to the Wilson & Walker's curve. This could be attributed to the fact that the empirical models come from data sets of past eruptions that include cases in windy conditions. Hence, the models of Sparks and Mastin actually include an influence from the wind, although not directly quantifiable. In the case of windy conditions (Figure 4b), the situation is the opposite, with WA models estimating higher MER values than NWA. The model of Wilson & Walker, though, approaches the prediction of the Woodhouse0D model at high plume heights (~18 km), which are much higher than the plume heights under analysis in this study. It is therefore very likely that, in windy situations like the ones under analysis, NWA models underestimate MER and vice versa. Finally, it should be noted that WA models depend on the plume centerline height rather than plume top height. For this study, though, there was no information on the plume radius and hence we used the plume top height measured by the ground-based radar to feed WA models. This leads to an overestimation of MER, since it would imply using the top plume height (which is always higher than the centerline height) as an estimate of the centerline height in the WA models. For example, for a plume height of approximately 8 km in a cross-flow of 20 m s −1 , the overestimation can be quantified as ranging from 3% to 15% for a radius ranging from 100 to 500 m for both WA models. These findings agree with those of Degruyter and Bonadonna (2012), where they compared their output MERs to those calculated using the Mastin model and time-averaged MER estimates obtained for the period 4-8 May 2010 by Bonadonna et al. (2011) and Bonadonna and Costa (2012). In particular, the REFIR-assessed MER avg (~1.6 × 10 5 kg s −1 ) compares well with the time-averaged MER obtained by estimating the erupted volume with the Weibull methodology (~2 × 10 5 kg s −1 , Bonadonna & Costa, 2012). We note that on the evening of 8 May the activity of Eyjafjallajökull temporarily dropped resulting in decreased plume heights and REFIR-assessed MER values (Figure 3). During that eruptive stage, the independent measurements obtained by Dürig, Gudmundsson, Karmann, et al. (2015) found values of MER ranging from 2.2 to 3.5 × 10 4 kg s −1 , based on video-analysis of the expansion of the discrete jets observed in video footages taken between 19:35 and 22:05 UTC. Given the considerable range of MER in the period in question, the MER values agree well with the REFIR output for MASTIN, NWA, and ALL ( Figure 4) and are on the lower end of range for estimates based on wind affected models only (WA).
Finally, it is worth noting that, having used the same weight factor for all the models in the case ALL, and as NWA models are more numerous than WA in REFIR (4 vs. 2), the resulting MER in the case ALL is more influenced by the estimate obtained with NWA models. This is also due to the fact that NWA and WA models tend to predict similar values of MER within the respective class of models, which can also be inferred from Figure 4. This effect could be avoided by first averaging the results of each class separately (NWA and WA) and subsequently merging these to obtain the final estimate of MER in the ALL case. This new

10.1029/2020JD032717
Journal of Geophysical Research: Atmospheres methodology of combining results of different type of models will be implemented in future versions of REFIR. However, the user can correct this effect by acting on the weights and, for example, assign a higher weight to WA models especially in situation of bent-over plumes. The User Manual gives some first indications on weight factor combinations depending on different scenarios, and this aspect is the matter of ongoing research.

NAME Outputs
We next consider the impact on the modeled extent of the ash cloud and the predicted peak concentrations, when using the MER time-series generated using the three plume height solutions: (H min , MER min ), (H avg , MER avg ), (H max , MER max ) and four different REFIR model configurations (ALL, WA, NWA, and MASTIN). Figure 5 shows the modeled ash cloud for 6 May 2010 over FL000-FL200 using H avg and MERs from the ALL, WA, NWA, and MASTIN configurations. Outputs for 5, 7, and 8 May 2010 and for FL200-350 are provided in the supporting information. Note that there is no ash in FL350-550 because the plume heights for this event are too low (Figure 1).
Not accounting for the wind when assessing the source strength for a weak bent-over plume has a significant impact on the forecast air concentrations and the extent of the predicted plume. Using the WA models in REFIR gives the highest MER, and as such when used to initialize NAME the greatest predicted ash concentration in the atmosphere. For the period 00:00-06:00 UTC the peak ash concentration is 11,605 μg m −3 , whereas using the MER from the ALL model configuration the peak is much lower, 6,831 μg m −3 . Figure 6 shows the difference in the modeled plume extent and concentration of ash for 06:00-12:00 UTC on 6 May over FL000-FL200 ( Figure 6a) and FL200-350 ( Figure 6b) when using (H min , MER min ), (H avg , MER avg ), and (H max , MER max ), and the ALL, WA, NWA, and MASTIN configurations. The biggest discrepancy between predicted peak concentrations in FL000-200 is observed when using the MASTIN MER min (459 μg m −3 ) versus the WA MER max (16,085 μg m −3 ). Using (H max , MER max ) ash is released over a greater vertical extent at the source and predicted peak concentrations in FL200-350 are higher. Using (H min , MER min ) results in the lowest forecast concentrations and very little ash is predicted above FL200 ( Figure 6b); modeled peak concentrations do not exceed 2,000 μg m −3 .
The results of the NAME modeling are summarized in Figure 7. It shows that the area where airborne ash concentration exceeds 2 mg m −3 and the maximum concentration of airborne ash are five and three times larger respectively for wind-affected models when compared to nonwind-affected models. The Mastin model, which is commonly used in volcanic ash advisory centers, is typical of nonwind-affected models. In both cases, the lowest estimates for the wind-affected models are similar to the highest values of the

Journal of Geophysical Research: Atmospheres
nonwind-affected models. The results for ALL models lies between those of wind-affected and nonwindaffected.

Discussion and Conclusions
The new capabilities of REFIR v19.0 have enabled the collation of time series of plume heights and MER, with their associated uncertainties, for a period of the 2010 explosive eruption at Eyjafjallajökull volcano. Summarizing the results, the new version of REFIR allows: 1. The estimation of ESPs (top plume height and MER) in reanalysis mode, hence enabling the analysis of time series of plume height data of past eruptions. 2. The usage of the Woodhouse et al. (2013) 0-D model, complementing the model of Degruyter and Bonadonna (2012) in the suite of WA models now available in REFIR; 3. The possibility of retrieving weather data from GFS or ECMWF repositories, allowing for the use of realistic weather conditions that are required for running WA models. 4. The output of time-series of input parameters for the NAME dispersion model, hence simplifying the process of preparing model runs.
The REFIR-assessed MER time series for the period 5-8 May 2010, obtained by using the plume height data from the C-band ground-based radar of Keflavík airport (Arason et al., 2011) have shown that the final estimate of the MER is strongly dependant on the plume model used and, more importantly, the type of model (WA and NWA). In this investigation, WA models produced MER estimates that are significantly higher than those outputs from NWA models (from 320% for the maximum solution to approximately 480% for the minimum solution). Similar trends were obtained for the weak plume scenario discussed in Costa et al. (2016). REFIR simulations considering both wind affected and no-wind-affected models (ALL) tend to smooth these differences out and the associated uncertainties. The time-averaged MER avg obtained in the ALL case compared well with published ground-truth data obtained from estimates of total erupted volumes (Bonadonna et al., 2011;Bonadonna & Costa, 2012), but are larger than the estimates obtained by Dürig, Gudmundsson, Karmann, et al. (2015), which are more comparable with MER avg obtained using NWA. Taken individually, NWA models would underestimate the MER compared to the ground-truth data, while WA models would overestimate them. The underestimation of MER from NWA models in the case of bent-over weak plumes is a well-known effect linked to the nature of the equation in these models . Using plume models to obtain H from MER at the source, NWA models would predict much larger H than WA models in case of windy conditions. Consequently, when using these models in the reverse way (i.e., from H to infer MER), the same plume height H would lead to a smaller value of MER using NWA than WA models. It is important to stress here that these conclusions apply to the eruptive scenario here presented and discussed; eruptions of the same size with no wind and larger eruptions less affected by winds would lead to opposite results, as also shown in the example of Figure 4.
Another interesting aspect is that the deviation of MER max from MER avg is smaller for the WA models than NWA; generally, curves of MER max and MER min are not evenly distributed around the average estimate MER avg , unlike the plume height estimates from which MER values have been obtained. Therefore, all models produce asymmetric values of MER around the best estimate due to the form of the equation, with the deviation increasing at increasing plume height. Currently, REFIR does not account for the intrinsic uncertainty of the plume models taken individually but mainly propagates the plume height determination's uncertainty into the MER computation of each single model and then combines the MER ranges of variation of all the model. We will explore the possibility to consider the single plume model's uncertainty in future versions of REFIR.
We have shown how uncertainty on the plume height propagates into MER estimation and, eventually, into the ash dispersion forecast predictions carried out with NAME, leading to uncertainty in both the predicted extent of the ash cloud and forecast peak concentrations, which needs to be accounted for when providing advice to local civil aviation authorities. With the new version of REFIR, it is straightforward to explore different modeling strategies for estimating MER in real-time and for past eruptions. In fact, REFIR can be used to consider differences in output from individual plume models or a combination of the models (ALL, NWA, WA) and to capture the uncertainty associated with both plume height monitoring and modeling approaches.
We initialized NAME with a scaled MER by taking just 5% of MER to represent the mass in the distal ash cloud, following the setup of the London VAAC. This assumption does not affect the results of this study, since we have not compared the results of dispersion simulations with observations, while instead focusing on the sensitivity of the model output to the input MER. However, it should be noted that this approach also represents a significant source of uncertainty in the setup of the dispersion model, which will also be propagated onto the predictions of the ash cloud extent and the peak concentrations. Further work is needed to validate modeled ash concentrations with observations, in which the uncertainty on both the modeled and observed mass loadings must be well constrained. Real-time assessment of the uncertainty on all the ESPs used to initialize operational dispersion models and observations of mass loadings in the distal ash cloud (e.g., from satellite and lidar retrievals) would allow operational meteorologists to better assign confidence in predictions of both the forecast extent and concentration of ash in the atmosphere in the future.
This study has identified and highlighted several current limitations of REFIR, which are being addressed by ongoing research. In particular, we are further investigating the best combination of the plume models' weight factors, based on the eruptive scenario (e.g., strong or weak plume). This will be achieved by comparing results obtained with REFIR against the ground truth (e.g., erupted mass, MER obtained with other methods) of well-studied past eruptions. Once a set of scenario-dependent weight factors has been defined, we will consider how to implement REFIR to automatically recognize the scenario, by combining weather and plume observations, and eventually assign the weight factors combinations automatically, but still leaving the operator with the possibility to manually control this process. Finally, a new method for combining results of the two classes of plume models (WA and NWA), consisting in combining the results of the models of each class taken separately before merging all the results, will be implemented in future versions as to avoid MER estimates being biased by the number of plume models in each class.

Data Availability Statement
Supporting data (including input and output data of REFIR and NAME simulations) can be retrieved at https://doi.org/10.5281/zenodo.3697186.