A multi-scale risk assessment for tephra fallout and airborne concentration from multiple Icelandic volcanoes – Part 1 : Hazard assessment

Introduction Conclusions References


Introduction
Evaluation of the tephra hazard is necessary to carry out comprehensive risk assessments of explosive volcanoes.The process is commonly divided into a succession of logical steps, including the identification and the characterization of eruptive deposits in the field, the development of comprehensive eruptive scenarios based on field observations and the use of models to quantify the hazard related to each eruptive scenario (e.g.Biass and Bonadonna, 2013;Biass et al., 2013;Bonadonna, 2006;Bonasia et al., 2011;Cioni et al., 2003;Connor et al., 2001;Costa et al., 2009Costa et al., , 2012;;Jenkins et al., 2012a;Leadbetter and Hort, 2011;Macedonio et al., 2008;Scaini et al., 2012;Volentik et al., 2009).
Although the hazard related to tephra dispersal and fallout rarely constitutes a direct threat to human lives, tephra can deposit on the ground up to hundreds of km away from the source and be dispersed thousands of km in the atmosphere.As a result, the impact from tephra varies with distance from the vent, resulting in complex vulnerability patterns of exposed elements (Blong, 1984;Connor et al., 2001).On the ground, tephra fallout can affect aspects such as human health, buildings, lifelines, economy or the environment.In the atmosphere, the residence time of ash can be as long as weeks, thus able to paralyse air traffic far away from the source, as demonstrated by the 2010 Eyjafjallajökull and the 2011 Puyehue-Cordón Caulle eruptions (Budd et al., 2013;Davies et al., 2010;Swindles et al., 2011;Wilkinson et al., 2012).Nonetheless, probabilistic studies of tephra dispersal tend to focus either on local ground deposition (e.g.Biass and Bonadonna, 2013) or on far-ranging atmospheric concentrations (e.g.Sulpizio et al., 2012), and only a few recent studies account for comprehensive multi-scale assessments (e.g.Scaini et al., 2012).
One crucial parameter for the description of the dispersal of tephra is the total grain-size distribution (TGSD).Centimetric to millimetric particles are controlled by gravitational settling and sediment in proximal to medial distances from the eruptive vent, whereas micrometric to submicrometric particles, when not aggregated, are controlled by larger-scale atmospheric processes and transported at continental scales (Folch, 2012).Depending on the assumptions used to model these two end-members, several modelling approaches have been developed to solve the advectiondiffusion equation with analytical, semi-analytical or numerical strategies (Bonadonna et al., 2012;Folch, 2012).
Eruption scenarios are usually developed for a single volcano and are constrained by the availability of past data and the completeness of the eruptive record (Marzocchi et al., 2004).Through time, the definition of eruption scenarios has evolved from a "worst-case scenario" approach towards an evaluation of the entire possible range of activity at a given volcano.For example, early hazard maps for Cotopaxi volcano (Hall and Hillebrandt, 1988;Vink, 1984) are based upon isopach maps of two major eruptions with opposite wind directions in agreement with the regional wind patterns and the most important exposed human settlements.More recent studies considered a probabilistic approach and developed a set of eruptive scenarios of various intensities based on accurate stratigraphic studies (Biass andBonadonna, 2011, 2013).Probabilistic techniques such as Monte Carlo simulations (e.g.Hurst and Smith, 2004) are nowadays an integral part of any hazard assessment for tephra dispersal and are used to investigate both the missing or inaccessible parts of the geological record and the impact of eruptions in a representative set of atmospheric conditions (Bonadonna, 2006).
For probabilistic modelling, the identification of eruption scenarios typically requires the definition of a probability density function (PDF) for each input parameter needed by a given model in order to account for the variability of eruptive processes (i.e.aleatoric uncertainty).For tephra fallout, several approaches have been used to define eruption scenarios, based on individual eruptions (Bonasia et al., 2011;Capra et al., 2008), eruptive styles (Macedonio et al., 2008), intensities/magnitudes (Scaini et al., 2012;Yu et al., 2013) or VEI classes (Biass and Bonadonna, 2013), mostly applied to a single source.However, some regions in the world are under the threat of more than one volcano, sometimes presenting a wide range of known eruptive styles and characteristics, and the development of comparable eruption scenarios for a set of volcanoes becomes an obvious necessity (e.g.Ewert, 2007;Jenkins et al., 2012a, b;Lirer et al., 2010;Simkin and Siebert, 1994).
Here, we present a medium-to long-term multi-scale scenario-based hazard assessment for ground tephra accumulation and far-range atmospheric ash dispersal from four Icelandic volcanoes -Hekla, Katla, Askja and Eyjafjallajökull -selected either for their high probabilities of eruption in the near future or their high potential impact (Fig. 1).Due to the different eruptive styles and the varying degree of knowledge of the eruptive history at these volcanoes, we developed consistent probabilistic eruption scenarios based on field data, literature studies and historical reports.The tephra-related hazard was assessed for each eruption scenario at a local scale (i.e.ground tephra accumulation) with the analytical model TEPHRA2 (Bonadonna et al., 2005) and at a regional scale (i.e.atmospheric concentration) with the numerical model FALL3D (Costa et al., 2006;Folch et al., 2009).A population of 10 years of wind obtained from reanalysis data sets was used to assess statistical atmospheric conditions.Outputs include (i) probabilistic maps of ground tephra accumulation and atmospheric concentration for relevant thresholds, (ii) mean atmospheric arrival time and persistence time, (iii) probability maps of atmospheric arrival time and persistence time for relevant thresholds, and (iv) ground hazard curves for critical facilities in Iceland.In a companion paper, Scaini et al. (2014) present a vulnerability assessment for both Iceland and the European air traffic system and use the outcomes of this study to perform a multiscale impact analysis.This comprehensive assessment is intended to act as a first step towards the elaboration of proactive measures for the management of explosive volcanic crisis.

Geological setting
Iceland is the result of the combined effects of a spreading plate boundary and a mantle plume (Allen et al., 1999;Vink, 1984;White et al., 1995;Wolfe et al., 1997).Current active volcanic zones (i.e. the neovolcanic zones) are the superficial expression of the mid-oceanic ridge.Arranged as discrete 15-50 km wide belts of active faulting and volcanism, they collectively cover a total area of 30 000 km 2 (Gudmundsson, 2000;Thordarson and Höskuldsson, 2008;Thordarson and Larsen, 2007).Volcanic zones host volcanic systems which, in their simplest forms, contain either a fissure swarm, a central volcano or both (Gudmundsson, 1995a, b).When present, the central volcano is the focal point of eruptive activity and the largest edifice in each system under which crustal magma chambers can develop, giving rise to either silicic or basaltic magmatism.In contrast, only basaltic magmas are erupted within fissure swarms (Larsen and Eirìksson, 2008;Thordarson and Larsen, 2007).

Target volcanoes
In this paper, we focus on eruptions occurring at the central volcanoes of four different volcanic systems located within two different volcanic zones: the Eastern Volcanic Zone (EVZ) and the Northern Volcanic Zone (NVZ; Fig. 1), Table 1.Historical eruptions at the central volcano of Hekla (see Thordarson and Larsen, 2007, and references therein).Tephra volumes are recorded as "freshly fallen" (i.e.40 % larger than volumes of old eruptions inferred from field mapping; Thorarinsson, 1967).Following the typical pattern of mixed eruptions at Hekla, plume heights correspond to the maximum altitude reached a few minutes after the onset of the eruption.NA: not available.ranked first and third in terms of volcanic activity in Iceland throughout the Holocene (Thordarson and Höskuldsson, 2008).No volcanic system was considered within the second most active volcanic zone, the Western Volcanic Zone (WVZ), because the EVZ is an active axial rift propagating southwards, thus taking over the activity of the WVZ (Mattsson and Höskuldsson, 2003;Thordarson and Larsen, 2007).The EVZ is divided into two sectors.In the north, the axial rift zone is characterized by a thick crust, high heat flow, well-developed tensional features and the production of tholeitic basalts.The southern propagating tip of the EVZ is often referred to as a flank zone, which lies on an older and thinner crust presenting a lower heat flow.Here, tensional features are poorly developed and the magma production consists mainly of transitional alkali to alkali basalts (Loughlin, 2002;Mattsson and Höskuldsson, 2003).We consider three volcanoes from the EVZ (Hekla, Katla, and Eyjafjallajökull; Fig. 1) and one from the NVZ (Askja).Amongst the selected volcanoes and in terms of activity since the settlement during the last 1100 years, Hekla is the most active with 18 eruptions from the central vent followed by Katla (18 eruptions), Eyjafjallajökull (3 post-17th century eruptions) and Askja (1 central vent eruption).Tables 1 and 2 summarize the recent eruptions of Hekla and Katla, respectively.The Supplement comprises a detailed description of  (Thordarson and Larsen, 2007).Uncompacted volumes are presented either as moderate (> 0.1-0.5 km 3 ) or large (> 0.5 km 3 ).the typical activity of these four volcanoes.The most striking features will be reviewed in Sect.4.1.together with the eruption scenarios.

Method
We aim to assess the hazard caused by the ground deposition and the atmospheric dispersion of tephra.Probabilistic approaches are adopted in order to account for the variability (i.e.aleatoric uncertainty) related to both eruptive and atmospheric conditions.We quantify the probability of hazardous mass load on the ground: where M(x, y) is the tephra mass load (kg m −2 ) accumulated at a given location and M T a mass load threshold (e.g.Bonadonna, 2006).For a given eruption scenario, the probability P M at a pixel (x, y) is quantified by counting the number of times a given threshold of load is reached over the total number of runs N: where For atmospheric concentrations, we start by quantifying where C(x, y, z, t) is the tephra mass concentration in the atmosphere (mg m −3 ) at a given point and time instant and C T a mass concentration threshold.For a given eruption scenario, the probability of disruption P C at a point (x, y, z) is quantified by counting the number of times a given mass concentration threshold is exceeded over the total number of runs N : where Disruption can be calculated at a given height or flight level (FL) or be comprehensive of all FLs, that is, considering that disruption occurs at a point (x, y) if the critical condition is achieved at any height z above the point.Note that for a given run, disruption occurs regardless of the number of model time steps during which Eq. ( 6) is verified.In order to provide a comprehensive picture of interest for air traffic management (ATM), we also quantify persistence and the first arrival time of a concentration threshold C T .The persistence is calculated by counting, for each run, the time interval in which the critical concentration threshold C T is exceeded at a pixel (x, y, z).The first arrival time (hereafter referred to as arrival time) computes the time from the beginning of the eruption to the first detection of the concentration C T at a pixel (x, y, z).Both persistence and arrival time were quantified as mean (i.e.average value at a pixel (x, y, z) over the entire number of runs), as standard deviation and as the probabilities of exceeding relevant thresholds of persistence and arrival times following the same approach as Eq. ( 5) at both specific FL and considering all vertical levels simultaneously.
In this section, we review (i) the models used at different scales, (ii) the probabilistic strategies adopted in this study, and (iii) the strategy used to account for particle aggregation processes.

Ground accumulation
The hazard related to ground deposition of tephra was assessed at the scale of Iceland using the steady semi-analytical advection-diffusion model TEPHRA2 (Bonadonna et al., 2005) following the approach detailed in Biass and Bonadonna (2013).TEPHRA2 requires five main input parameters: plume height, eruption duration, erupted mass, TGSD and particle density.It also requires a vertical wind profile, a calculation grid and three empirical parameters: a fall-time threshold acting as a threshold for the modelling of the diffusion of small and large particles (i.e.power law vs.linear diffusions), a diffusion coefficient used for the linear diffusion law and an apparent eddy diffusivity fixed at 0.04 m 2 s −1 for the power law diffusion (Bonadonna et al., 2005;Suzuki, 1983;Volentik et al., 2009).Empirical parameters can be estimated either using TEPHRA2 in inversion mode when enough field data are available (Connor and Connor, 2006;Volentik et al., 2009) or using analogue eruptions.Wind conditions for the 2001-2010 period were extracted from the ECMWF Era-Interim Reanalysis data set at a 1.5 • resolution.The calculation grid covers the small computational domain (Fig. 1) at a resolution of 1 km.When needed, smaller calculation grids were used at a resolution of 500 m.
All simulations compute 168 h of dispersal on a 120 × 60 grid with a horizontal resolution of approximately 0.5 degrees.Vertical layers are defined from 0 to 36 000 m with a variable vertical interval.The mass distribution follows Suzuki (1983) and Pfeiffer et al. (2005), with A = 4 and λ = 5.Such values were chosen since (i) all eruptions are of subplinian/Plinian types and (ii) only the fine fraction is modelled with FALL3D to assess the far-range dispersal.The terminal velocity is given by Ganser (1993), with a sphericity fixed to 0.9.The vertical diffusion coefficient was fixed to 10 m 2 s −1 and the horizontal diffusion coefficient was calculated using the CMAQ model parameterization (Folch et al., 2009).Meteorological fields for the considered period were extracted from the ECMWF Era-Interim Reanalysis at 1.5 • every 3 h, covering northern and central Europe (Fig. 1).Outputs are produced at a regular 5000 ft vertical interval from FLs 050-400.

Probabilistic strategies
Several approaches exist to assess the probability distribution of reaching a hazardous accumulation of tephra given an eruption (Bonadonna, 2006).In order to account for variable parameters (i.e.eruptive and atmospheric conditions), a large number of model runs are performed varying input parameters, including eruption date (i.e.wind profile for TEPHRA or 4-D variables for FALL3D).Each run consists either of a single occurrence of the model (i.e.short-lasting eruptions) or a set of simulations in time (i.e.long-lasting eruptions).When an approach with variable eruptive parameters is adopted, a PDF must be defined to constrain the stochastic sampling.The definition of the PDF, which reflects the knowledge of the system, is relevant to the definition of eruptive scenarios and will be tackled later.

One eruption scenario (OES)
The OES is an approach used to compile the probability of reaching a given threshold of tephra accumulation in variable wind conditions, with ESP chosen deterministically.Figure 2a summarizes the algorithm applied to short-lasting eruptions (Bonadonna, 2006).First, the plume height, the eruption duration, the total mass and the TGSD are fixed deterministically.Then, for each single run of the model, an eruption date is sampled from which the corresponding wind conditions are extracted from the meteorological database.

Eruption range scenario (ERS)
The ERS approach allows for a stochastic sampling of ESP at each run as well as wind conditions, where each variable parameter is characterized by a range and a PDF (Bonadonna, 2006).Figure 2b shows the algorithm used for the ERS.First, a plume height, an eruption duration and an eruption date are sampled from their respective PDF, and a maximum total mass for the eruption scenario is set.From the eruption date, the respective meteorological conditions are loaded, which, combined with the plume height, allows for the calculation of the MER using the method of Degruyter and Bonadonna (2012).A test is then performed to assess whether the resulting mass, calculated by combining the MER and the eruption duration, fits into the initial assumptions of mass range.If the test is negative, all parameters are resampled, otherwise the selected input parameters are sent to the model.

Long-lasting eruptions
New eruption scenarios were developed to assess the hazard related to long-lasting eruptions (Fig. 2c and d).For longlasting eruptions, ESPs are expressed as time series at constant time intervals, t, defined depending upon the availability of data (i.e.measurements of plume height, wind profiles).The application of algorithms shown in Fig. 2c and d varies depending on the scale of the hazard assessment, and thus the model used.When used with steady models (e.g.TEPHRA2), they consist of discrete model runs at constant time intervals, and the final hazard maps are the sum of all individual runs (e.g.Scollo et al., 2013).When used with nonsteady models (e.g.FALL3D), ESPs are updated at a constant time interval.For clarity, we will refer to any single run or update of the model as "occurrence", i.e. for long-lasting eruptions, a run (i loop in Fig. 2) consists of several occurrences (j loop in Fig. 2).

Long-lasting one eruption scenario (LLOES)
The LLOES relies on the same concept as the OES, i.e. eruptive parameters chosen deterministically with varying wind conditions, only applied to long-lasting eruptions.Here, the total eruption duration, the time series of plume heights and the TGSD are set deterministically, and the time interval ( t) is set based on the availability of data.At each run of the model, an eruption date is sampled and the corresponding wind conditions are extracted based on the eruption duration and t.At each occurrence, knowing t, the plume height and the wind conditions at the given time, the MER and the mass are calculated averaged over the interval length, and a new occurrence of the model is performed.Each run is the sum of all occurrences performed.

Long-lasting eruption range scenario (LLERS)
The LLERS applies the ERS strategy to long-lasting eruptions.Eruptive parameters are stochastically sampled as time series at a constant t.Following the algorithm in Fig. 2d, a mass boundary is first set.At each run, an eruption date and an eruption duration are sampled.At each occurrence, a plume height is sampled from a PDF, and the MER is calculated based on the wind conditions for that specific date using the method of Degruyter and Bonadonna (2012).If the sum of the mass of all occurrences of a single run falls out of the initial mass assumptions, the sampling process is restarted.
If not, input values for all occurrences are sent to the model.Eventually, each run is the sum of all occurrences performed.

Ash aggregation
Aggregation processes are known to modify deposition trends along the dispersal axis by aggregating fine particles (typically < 100 µm) into larger clusters (Brown et al., 2012;Rose and Durant, 2011).The aggregation of fine particles into larger aggregates results in premature sedimentation in the proximal area and in a relative depletion of fines away from the vent (Carey and Sigurdsson, 1982;Hildreth and Drake, 1992), possibly leading to an underestimation of the hazard in proximal areas and an overestimation in distal sectors.The fallout of aggregates has been characterized during numerous eruptions, including the 1980 eruption of Mount St. Helens (Carey and Sigurdsson, 1982;Durant et al., 2009), the ongoing eruption of Montserrat (Bonadonna et al., 2002) and the 2010 eruption of Eyjafjallajökull (Bonadonna et al., 2011;Taddeucci et al., 2011).Although aggregation is a topic of intense research, no satisfactory parameterization of this process has been achieved yet (e.g. Brown et al., 2012;Costa et al., 2010;Folch et al., 2010;Van Eaton and Wilson, 2013;Van Eaton et al., 2012).Several models attempt to describe aggregation processes using either empirical (Bonadonna and Phillips, 2003;Carey and Sigurdsson, 1982;Cornell et al., 1983) or numerical approaches (Costa et al., 2010;Folch et al., 2010).Here, we use the empirical approach of Bonadonna et al. (2002) and observations of Bonadonna et al. (2011) to modify the TGSD before running the models.Following this strategy, we remove an equal proportion of masses of fine particles from phi classes ≥ 4 , which are equally redistributed between classes −1 and 3 .The amount of fine particles removed, i.e. the aggregation coefficient, is stochastically sampled between 20 and 80 % on a uniform distribution at every loop increment on the algorithms shown in Fig. 2. The choice of such a range is based on field observations presented in Bonadonna et al. (2002Bonadonna et al. ( , 2011)).  1) from the ERA-Interim Reanalysis database.These wind roses show the probability that the wind will blow in given directions and speed intervals.

Identification of eruption scenarios
The identification of eruption scenarios is based upon the eruption history presented in the Supplement for each volcano.Here, only the historical catalogue of eruptions was used for three reasons.Firstly, with a mean eruption of > 20 events per century (Thordarson and Höskuldsson, 2008), reports and eye-witness accounts since human settlement constitute a valuable source of information when developing eruption scenarios.Trying to merge such a data set with an eruption catalogue based on stratigraphic studies only results in discrepancies in the completeness of the eruptive record, making any comparison difficult.Secondly, the climate of Iceland is prone to fast erosion of freshly fallen deposits, inducing a bias towards large eruptions when trying to assess the prehistoric eruptive activity.Thirdly, since it is recognized that mean eruption frequencies are strongly correlated to glacier load (Albino et al., 2010), any attempt to assess eruptive patterns during older time periods might not be representative of the actual climate and load.As a result, the hazard assessment presented here implies that a future eruption at any of the four volcanoes will follow behaviours similar to the eruptive style displayed in historical times.
Additionally, we only focus on the activity occurring at the central vent of the selected volcanic systems as no explosive eruptions occurred from fissure swarms.Both ground accumulation and atmospheric dispersal were tackled probabilistically, using the models TEPHRA2 and FALL3D, respectively.Since statistical analysis shows neither monthly nor seasonal trends, we used a population of 10 years of wind data (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010).Figure 3 shows wind roses at three altitude levels for points specified in Fig. 1.

Hekla
Out of all volcanoes considered in this study, Hekla presents the most accurate historic record with 18 identified and reasonably well-described eruptions, a good characterization of 5 of them and a TGSD of the 2000 eruption (Gronvold et al., 1983;Gudmundsson et al., 1992;Höskuldsson et al., 2007;Thorarinsson, 1967;Thorarinsson and Sigvaldason, 1972).From the 18 eruptions presented in Table 1, we discarded the eruption of 1104 from the eruptive record as it belongs to the larger magnitude and lower frequency H1-5 series and occurred after a repose interval of > 230 years, more than twice the maximum repose time in historical times (Larsen and Eirìksson, 2008;Thordarson and Höskuldsson, 2008;Thordarson and Larsen, 2007).Since a large majority of the tephra is produced during the opening phase of "mixed" eruptions (Supplement), we only model the first part of the eruption.
Based on the historical period summarized in Table 1, two eruption scenarios were identified for Hekla and named after the best-described eruption of each cluster.First, a 2000-type eruption was identified to describe eruptions that have occurred at a ∼10-year frequency.This eruption scenario considers all available data for the eruptions of 1222, 1970, 1981, 1991 and 2000.Second, a 1947-type eruption is considered representative of large Plinian eruptions occurring on an average of one per century starting from the eruption of 1158 (Thorarinsson, 1967).

2000-type scenario
The good knowledge of eruptions considered in the 2000type scenario allows for the identification of well-constrained ranges of ESP for plume heights and erupted masses, and, subsequently, the use of the ERS strategies.Table 3 and Fig. 4 summarize ranges of ESPs and their corresponding PDFs.In agreement with historical witnesses, the duration of the intense opening phase was fixed between 0.5 and 1 h and sampled on a uniform PDF.Based on Table 1 and following the algorithm presented in Fig. 2, a mass constraint was fixed between 6.9 × 10 9 and 6.9 × 10 10 kg considering a deposit density of 691 kg m −3 .Based on observations, plume heights were defined between 6 and 16 km a.s.l. and sampled on a logarithmic scale in order to slightly favour small events over large ones.
The TGSD for the 2000 eruption shows bimodality with peaks at −3 and 3 (Fig. 4c) and a minimum at 1 .This TGSD was used as representative for the 2000-type eruption scenario and was varied at each run.First, considering 1 as a boundary between the two populations, we varied the relative weight percentage of the two populations between 30 and 70 %.Second, we applied the aggregation model described in Sect.3.3 on the resulting modified TGSD by removing between 20 and 80 % of fines.Figure 4c shows the three states of the TGSD for one of the 1000 runs.

1947-type scenario
A similar approach was applied to the 1947-type scenario.
Since only sparse estimates of plume heights are reported in the literature for this eruption type, the sampling was constrained between 16 (i.e. the highest boundary of the 2000-type eruption) and 30 km a.s.l. on a logarithmic PDF (Fig. 4d).Based on Table 1, the mass constraint was defined between 6.9 × 10 10 and 3.5 × 10 11 kg (Fig. 4e).The resulting distribution of plume heights is shown in Fig. 4d and displays a maximum of sampling around 18-20 km a.s.l.No solution compatible with the mass constraint was found for plumes higher than 27 km.Figure 4e shows a rather uniform mass distribution, though slightly biased towards lower values.
Since no sufficient measurements exist to infer the TGSD for any of these eruptions, we consider here a Gaussian distribution ranging from −5 to 8 .We allowed a variability of Md and σ on uniform PDF between −1 to 1 and 1 to 2 , respectively.The aggregation model was applied, with aggregation coefficients varying between 20 and 80 % for all bins ≥ 4 .

Katla
Numerous eruptions from Katla have been well described and documented, but only a few quantitative constraints exist.Based on Table 2 and Larsen (2010), about 10 historical eruptions produced tephra volumes > 0.1 km 3 , with only the 934-940 Eldgjá eruption responsible for a volume > 1 km 3 .Since the Eldgjá eruption originated from the surrounding fissure swarm rather than the central volcano, we discard it from the eruption record used in this study, resulting in relevant tephra volumes ranging from 0.1 to 1 km 3 .Historical eruptions at Katla are known to have lasted from 2 weeks to 5 months, with most of the tephra produced during the first days.No silicic eruptions were witnessed during historical times, with the last SILK layer erupted 1675 years BP (Larsen et al., 2001).
As a result, a LLERS strategy was applied to Katla volcano in order to assess the hazard related to a future moderate to large basaltic eruption (Table 2).According to the existing literature, plume heights were sampled between 10 and 25 km a.s.l. on a logarithmic PDF (Einarson et al., 1980;Larsen, 2000Larsen, , 2002;;Óladóttir et al., 2008, 2006, 2011;Thordarson and Höskuldsson, 2008).Only the paroxysmal phase was modelled and assumed to last between 1 and 4 days stochastically sampled on a uniform PDF.A volume constraint was set between 0.1 and 1 km 3 , converted into a mass constraint between 0.7 and 7 × 10 12 kg using a bulk density of 700 kg m −3 .Since no TGSD is available, we used a reconstructed TGSD from the 10 points available for the 1357 eruption in the study of Einarson et al. (1980) using the method of Bonadonna and Houghton (2005), which results in a Md of −1 and a σ of 2. However, due to the southward dispersal of the eruption and the narrow area of land between the volcano and the sea, these points present a proximal cross section of the deposit.As a result, a Md of −1 is considered as a maximum value, and the TGSD adopted here is a Gaussian distribution between −7 and 8 with Md sampled between −1 and 1 , and σ sampled between 1 and 2 , both on uniform PDF.The resulting distribution is aggregated with an aggregation coefficient sampled between 20 and 80 %.The same TGSD is used for all occurrences of a given run.The time interval between two occurrences was set to 6 h based on the availability of wind data from the ECMWF database.
Figure 5 and Table 3 summarize the ESP for Katla.Figure 5a shows the resulting PDF for plume heights displaying a slight logarithmic trend, and Fig. 5b shows eruption durations.Figure 5c  individual occurrences and results in a strongly logarithmic shape with masses comprised between 10 9 and 6 × 10 11 kg.
Figure 5d shows the resulting PDF for the mass per run, and Fig. 5e the TGSD at one of the 1000 runs.

Eyjafjallajökull
The limited knowledge of eruptions at Eyjafjallajökull constrains the identification of eruption scenarios.Two prehistoric eruptions are recognized in the field but poorly constrained and three post-17th century eruptions were witnessed, amongst which the eruptions of 1612 and 1821-1823 lack any constraint.However, since detailed observations and measurements of eruptive parameters exist for the 2010 eruption, we applied here a LLOES strategy in order to assess the entire range of possible hazard related to the occurrence of a similar eruption.
We model here the 40 days of explosive phase that occurred from 14 April to 20 May 2010.The algorithm used is shown in Fig. 2. The total duration, the time series of plume heights and the TGSD are deterministically set a priori.Figure 6a shows measurements of plume heights every 6 h for the 40 days of eruption (Arason et al., 2011), converted into MER using the method of Degruyter and Bonadonna (2012) and wind conditions extracted from the ECMWF database (Fig. 6b).As a result, the time interval between occurrences within a run was set to 6 h.The TGSD used here is as described by Bonadonna et al. (2011), who reconstructed disaggregated and aggregated TGSD by combining ground-based and satellite-based measurements.Here, the same TGSD is used for all runs and does not vary through occurrences.
Following the algorithm in Fig. 2, an eruption date is sampled at each run, after which wind conditions are extracted for the 40 days of the eruption every 6 h.At each occurrence, the MER is calculated accounting for the wind velocity and converted to 6 h averaged mass.The occurrence is sent to the model, and each run is the sum of the 240 occurrences.
For computational reasons, it was not possible to run the Eyjafjallajökull probabilistically with FALL3D.

Askja
At least two large tephra deposits associated with strong Plinian eruptions of VEI 5 are recognized at Askja including the 10 ka BP and the 1875 eruptions.Since no accurate mapping or constraints of eruptive parameters are available for the 10 ka BP eruption, we use the 1875 as a reference eruption.Previous studies of Sparks et al. (1981) and Carey et al. (2010) provide an accurate chronology of the different phases of the eruption as well as constraints of the associated eruptive parameters.Two phases are responsible for most of the production of tephra, namely the 1 h long phreato-Plinian phase Askja C followed 6 h later by the 1.5 h long Plinian phase Askja D. As a result, we apply here OES strategies modelling two consecutive eruptions separated by a 6 h break.
Figure 2 shows the algorithm developed for a single OES modelling.Here, the hazard related to a 1875-type eruption consists of the sum of one OES for Askja C and one OES for Askja D. All ESPs (i.e.plume height, erupted mass, eruption duration and TGSD) are fixed deterministically and are summarized in Table 3 and Fig. 7 for both phases.At each run, an eruption date is sampled and wind data for the consecutive phases are extracted.Both TGSD are aggregated with an aggregation coefficient sampled between 20 and 80 %.

Hazard assessment
This section presents the results of the different model runs for all volcanoes.These results should be viewed as part of a scenario-based hazard assessment, showing the geographical distribution of probability of tephra fallout knowing that an eruption of a given magnitude is occurring.As described in Sect.3.2, the compilation of probability maps requires a threshold -i.e.either a mass load (kg m −2 ), a concentration (mg m −3 ), a persistence or an arrival time (h) -in order to calculate the exceeding probability for each eruption scenario.Based on the available literature, we use three relevant thresholds for ground accumulation, one for atmospheric concentration, one for persistence and one for arrival time (Table 4).The Supplement contains the entire collection of maps produced including probability maps for atmospheric concentration as well as persistence and arrival times in terms of mean, standard deviation and probability, computed for all FL and for all critical thresholds.The resulting impact is presented in the companion paper of Scaini et al. (2014).

Ground deposition
Figures 8-11 show the probability maps of exceeding a given threshold of tephra accumulation on the ground for Hekla, Katla, Eyjafjallajökull, and Askja, respectively.In agreement with Fig. 3, there are preferential eastwards dispersals, leaving the Reykjavík area with a negligible probability of being affected by tephra fallout for the volcanoes considered here.As a result, eruptions from volcanoes in the EVZ are likely to affect the area between and east of Gullfoss and Vík í Mýrdal (hereafter referred to as Vík; Fig. 1) .
Figure 8 shows the probability maps for Hekla.Figure 8a  and b show the probability of reaching an accumulation of 1 kg m −2 for the 2000-type and the 1947-type scenarios, respectively.In the case of a 2000-type eruption, there is a > 10 % probability of reaching such an accumulation up to 50 km east of the volcano and a negligible probability to affect Vík.However, Vík and the southernmost coast have a ∼ 15 % probability of reaching an equal accumulation for a 1947-type eruption, with the > 10 % probability line extending 150-200 km eastwards and 50 km westward from the volcano (Fig. 8b).This scenario has a ∼ 10 % probability of producing tephra accumulation of 10 kg m −3 in the vicinity of Gullfoss (Fig. 8c) and has a 10 % probability of affecting an area with an accumulation of 100 kg m −3 25 km east of the volcano (Fig. 8d).
Figure 9a-c show the spatial distribution of probabilities of reaching tephra accumulations of 1, 10 and 100 kg m −2 , respectively, associated with an eruption at Katla.At Vík, such an eruption results in probabilities of 40 %, ∼ 30 %  (Arason et al., 2011); (b) corresponding MER for the same period based on wind conditions extracted from the ERA-Interim database and the method of Degruyter and Bonadonna (2012); (c) disaggregated and aggregated TGSD as inferred by Bonadonna et al. (2011).See caption of Fig. 4 for explanation of symbols.and 10 % of reaching tephra accumulations of 1, 10 and 100 kg m −2 , respectively.The 10 % probability line of reaching a tephra accumulation of 1 kg m −2 extends about 200 km northwards on the mainland and eastwards along the coast.
Figure 10 displays the probability distribution for an Eyjafjallajökull 2010-type eruption, resulting in 80 and 20 % probabilities of reaching tephra accumulations of 1 and 10 kg m −2 in Vík, respectively (Fig. 10a and b).Due to the low probability level, the map for an accumulation of 100 kg m −2 is not shown.
Located in the NVZ, Askja is most likely to impact the eastern part of the country, with half of the territory having a 5-10 % probability of reaching a tephra accumulation of 1 kg m −2 should an 1875-type eruption occur (Fig. 11a).The main town under the threat of an eruption of Askja, Egilsstaðir, has ∼ 35 and ∼ 15 % probabilities of reaching tephra accumulations of 1 and 10 kg m −2 , respectively (Fig. 11a and b).The towns of Akureyri and Húsavik, which both have airports used for internal flights, have 15 and 20 % probabilities of reaching tephra accumulations of 1 kg m −2 , respectively.A 1875-type eruption also has a 10 % probability of depositing 100 kg m −2 of tephra 50 km east of the vent (Fig. 11c).
Along with probability maps, the hazard related to tephra accumulation can also be expressed as hazard curves, for which the probability of exceeding any tephra accumulation is quantified for a given location.Figure 12 shows hazard curves for relevant eruptions for the locations of Vík, Gullfoss, Akureyri and Egilsstaðir.Although Gullfoss is only a tourist facility, this location was used to assess the probability of tephra accumulation inland.Figure 12a shows that Vík has 15, 40 and 80 % probability of exceeding a tephra accumulation of 1 kg m −2 following the considered eruptions of Katla, Hekla 1947-type and Eyjafjallajökull, respectively.

Atmospheric concentration
The hazard assessment for atmospheric concentration is shown in Figs. 13 and 14, Table 5 and the Supplement.For practical reasons, we present here only probability maps accounting for the presence of ash above a threshold of 2 mg m −3 at any FL (Sect.3, Table 4).Arrival time and persistence are compiled here in the form of probability maps of exceeding arrival and persistence times of 24 and 12 h, respectively.The choice of 24 h for the threshold of arrival time is based on Guffanti et al. (2010), who showed that 89 % of the aircrafts encounters with volcanic ash in the period 1953-2009 occurred within the first 24 h after the onset of the eruption.Since no threshold of persistence time has been officially outlined, we adopted a threshold of 12 h based on qualitative observations found in the literature (e.g.Ulfarsson and Unger, 2011).The Supplement comprises probability maps for other thresholds (i.e. 2 × 10 −6 and 0.2 mg m −3 ) and for separate FL as well as maps of mean and standard deviation of persistence and arrival time.Due to the high computing demand required to run FALL3D in a probabilistic mode for a 40 day long eruption, the Eyjafjallajökull LLOES 2010type eruption was omitted from the hazard assessment for atmospheric dispersal.
Figure 13a-c show the results for Hekla.Due to the local dispersal following a 2000-type eruption, only maps for a 1947-type eruption are presented here.Such an eruption would result in a 5-10 % probability of reaching a concentration of 2 mg m −3 above the northern Atlantic Ocean and probabilities of reaching London and Oslo of 0.8 and 0.5 % after 13 ± 3 and 17 ± 5 h, respectively (Fig. 13a and b, Table 5).Persistence times for these locations are negligible and are of 3 ± 2 and 5 ± 2 h, respectively (Fig. 13c, Table 5).As shown in the Supplement, similar observations can be made at all separate FLs.These results are in agreement with the findings of Leadbetter and Hort (2011), although differences in scenario identification, modelling strategies and definition of critical thresholds make detailed comparison difficult.
Following the scenario used here, an eruption at Katla has a 5-20 % probability of affecting the UK and Scandinavia with concentrations of 2 mg m −3 .Such concentrations are expected to arrive above London and Oslo after ∼ 45 ± 22 h in both cases but persist in the atmosphere for less than 10 h (Fig. 13d-f, Table 5).Due to the long-lasting nature of the Katla scenario, note (i) the longer arrival time compared to other eruption scenarios and (ii) the high uncertainty on the mean persistence and arrival time.The Supplement shows that at separate FL, the impacts of a Katla eruption tend to decrease with altitude.
The atmospheric dispersal of tephra following a 1875-type eruption at Askja is presented in Fig. 13g-i, which results in 5-20 % probabilities to reach a concentration of 2 mg m −3 over Scandinavia and Western Europe (UK, Northern France, Netherlands, Belgium, and Western Germany).Such a concentration has a 5-10 % probability of reaching the UK and Scandinavia within 24 h, with mean arrival times above London and Oslo of 22±13 h and 25±12 h, respectively (Fig. 13g  and h, Table 5).The airports of Paris and Frankfurt can potentially be impacted after ∼ 50 ± 20 h in both cases.In all cases, persistence in the atmosphere would be in the range of ∼ 6-8±4 h, with a 5-10 % probability of western Norway to be locally impacted by concentrations of 2 mg m −3 persisting for more than 12 h.Similar probability distributions, arrival and persistence times are to be expected at all separate FL (see the Supplement).

Short vs. long-lasting eruptions
In order to compare the potential impact resulting from different types of activity at the selected volcanoes (i.e.shortvs.long-lasting eruptions), this section presents deterministic scenarios based on historical and well-constrained eruptions.Note that these simulations do not aim at presenting "worstcase" scenarios, which would require the combined identification of worst-case eruption scenarios and wind condition, but can be viewed as a comparison of key historical eruptions happening under similar meteorological conditions.
The eruptions of Hekla 1947, Katla 1918, Eyjafjallajökull 2010, and Askja 1875 were selected as case-study scenarios for which sufficient data were available to produce a realistic forecast of potential impact.In order to scale and compare the effect of these eruptions, simulations were run using the wind conditions of Eyjafjallajökull 2010, starting from 14 April and lasting for 10 days.ESPs for Eyjafjallajökull and Askja are summarized in Table 3.The 1918 eruption of Katla is the most recent eruption to break through the Mýrdalsjökull ice cap.The available literature suggests that the eruption lasted for 3 weeks, with the most intense tephra production during the first days, plume heights up to 14 km a.s.l. and a total volume varying between 0.7 and 1.6 km 3 (Larsen, 2000;Sturkell et al., 2010).In the absence of any detailed variations of plume heights, radar observations of Arason et al. (2011) were used and scaled to fit observed minimum and maximum plume heights of the Katla 1918 eruption.Using wind conditions specified above and the method of Degruyter and Bonadonna (2012) to estimate the MER, we obtained a total mass of 1.24×10 12 kg, which is consistent with published volume estimates.Given the similarities between the two systems (Sturkell et al., 2010), the TGSD of Eyjafjallajökull defined by Bonadonna et al. (2011) was also used for this run.ESP for the Hekla 1947 eruption were set using the literature, with a plume height of 27 km a.s.l., a total erupted tephra volume of 0.18 km 3 and a duration of 1.5 h.
Figure 14 summarizes the expected concentrations at FL150 over the main European airport hubs of London Figure 14 shows that a Hekla 1947-type eruption bears the largest impact in terms of airborne concentration.A 1947type eruption would result in concentrations above the threshold of 0.2 mg m −3 over London, Paris, Amsterdam, Frankfurt and Copenhagen (i.e.0.35-0.68mg m −3 ) arriving after 2-4 days and potentially disrupting the air traffic between 1 and 3 days.Oslo would be the most impacted area with concentration peaks above 10 mg m −3 .A Katla 1918type eruption, characterized by a pulsatory regime with repeated emissions of ash, would potentially be most problematic for the European airspace in terms of duration of disruption.Although most likely reaching concentrations comprising only between 0.1 and 0.25 mg m −3 , Fig. 14 reflects the repeated arrival of clouds over time and its potential implications for the management of a volcanic crisis.The eruptions of Askja 1875 and Eyjafjallajökull 2010 only seem to be problematic for Oslo, reaching concentrations above 2 and 0.2 mg m −3 respectively.

Eruption scenarios and probabilistic strategies
Three steps were taken to develop probabilistic hazard scenarios including (i) the identification of the most probable and potentially problematic eruptive styles at given volcanoes, (ii) the development of adapted algorithms to model each eruption type (Fig. 2), and (iii) the definition of eruption scenarios with constrained ESPs (Figs. 4-7).At the selected volcanoes, this led to the identification of both shortand long-lasting eruptions scenarios with both fixed and variable ESP.

Fixed vs. variable ESPs
Amongst chosen volcanoes, a large discrepancy exists in terms of the knowledge of the eruptive history, which is directly related to the frequency of activity during historical times: 17, 10, 3, and 1 eruptions of interest (i.e.explosive eruptions at the central vent) occurred during historical times at Hekla, Katla, Eyjafjallajökull and Askja, respectively.On the other hand, a discrepancy also exists in the degree of detail to which eruptions have been mapped and characterized.For example, the single eruption of Askja is thoroughly characterized in terms of chronology of eruptive phases, plume height, erupted volume and TGSD whereas eruptions of Katla are mainly bounded by rough estimates of volume.As a result, the choice of expressing eruption scenarios as either a single set of ESPs deterministically defined (OES) or as a stochastic sampling on a PDF (ERS) is made upon the combined knowledge of the eruptive history and the degree of characterization of eruptions.The end-user should account for the limitations considered with each approach.3 and Fig. 7.

Sampling of ESPs
Figure 2 shows the algorithms used to produce both ERS and LLERS, where the erupted mass is indirectly derived from the plume height, the eruption duration and wind conditions (Degruyter and Bonadonna, 2012) and tested against a mass range.Although the sampling of plume height was constrained on a logarithmic distribution as a prior knowledge, the resulting PDF only including values validated by our algorithm shows a wide variety of shapes.For example, Fig. 4a shows that the initial assumptions of erupted mass (i.e.6.9 × 10 9 -6.9 × 10 10 kg) for a 2000-type eruption at Hekla for a 0.5-1 h long eruption cannot be realistic with plume heights under 10 km a.s.l., resulting in (i) a PDF for plume heights biased towards the largest end-members and (ii) a PDF for erupted mass in agreement with an initial assumption of a logarithmic distribution of ESP.Similarly, the 1947-type scenario results in a PDF with a maximum at plume heights of 18-20 km a.s.l.but without solution for plume heights above 27 km a.s.l.satisfying the initial mass (6.9 × 10 10 -3.5 × 10 11 kg) and duration (0.5-1 h) conditions (Fig. 4d).As a result, this method accounts for a prior knowledge of the system (i.e. initial choice of a PDF for the sampling of ESP) and helps to correct the sampling of dependent ESPs (i.e.plume height, eruption duration, MER and erupted mass) in order to produce realistic events within an eruption scenario.

Short-vs. long-lasting eruption scenarios
Assessing the hazard related to tephra dispersal from longlasting eruption is commonly done using non-steady models but rarely using steady models.Scollo et al. (2013) already used the model TEPHRA2 to evaluate tephra hazard associated with long-lasting violent Strombolian activity at Mt Etna, Italy (e.g. the 21-24 July 2001 eruption).They defined it as weak long-lived plume scenario (OES-WLL and ERS-WLL) with an eruption duration of 4-100 days, in contrast to short-lived plume scenarios (OES-SSL and ERS-SSL) associated with the paroxysmal phase of subplinian eruptions (e.g. the 22 July 1998 eruption).
Here, we developed new algorithms for the sampling of ESPs to assess the ground deposition from long-lasting eruptions (Fig. 2).Conceptually, the total ground accumulation calculated with TEPHRA2 consisted of consecutive occurrences of the model run at a given time interval t, after which all outputs are summed.With FALL3D, the continuous computation allowed simply updating of the ESPs every t without interruption.Here, the typical 6 h time resolution of reanalysis data sets conditioned the duration of t, implying constant eruption conditions between either different runs or updates.

Ground accumulation
Figures 8-11 show that although computed accumulations are not sufficient to cause structural damage to buildings (i.e.> 100 kg m −2 ), deposition of 1-10 kg m −2 is likely to occur, which is consistent with historical chronicles primarily reporting impacts on agricultural activities (e.g.crops destruction, poisoning of animals; Thorarinsson, 1967;Thorarinsson and Sigvaldason, 1972).In addition, eruptions from ice-capped volcanoes such as Katla and Eyjafjallajökull are typically associated with jökullhaups, which can cause significant structural damage to buildings, roads and bridges.If these observations are valid for the selected volcanoes and potentially for most of central vent eruptions with VEI up to 5 -excluding maybe the volcanoes located in the vicinity of Reykjavík and Keflavik -"fires"-type eruptions would result in larger magnitude impact.A review of the environmental changes produced by the Eldgjá fires can be found in Larsen (2000).
Hazard maps produced here show preferential dispersals towards the E-ENE, consistent with wind observations (Fig. 3).However, compilations of dispersal axes for historical eruptions are available for Hekla (Thorarinsson and Sigvaldason, 1972) and Katla (Larsen, 2000) and show the existence of deposition in all directions around the vents.For example, only 7 out of the 14 historical eruptions of Hekla were dispersed in directions between 0 and 180 • , and 8 eruptions dispersed tephra in a 340-20 • sector.Similarly, half of the historical eruptions of Katla were dispersed and deposited with a bearing comprised between 0 and 180 • .By comparing our probability maps for Eyjafjallajökull to the isomass maps of Gudmundsson et al. (2012) compiled for land deposition in the period 14 April to 22 May, we observe a ground deposition slightly more directed towards the south than predicted by our model (Fig. 10).However, deposition observed on the ground for both 1 and 10 kg m −2 (converted from the isopach maps of Gudmundsson et al. (2012) with a density of 1400 kg m −3 ) fall between our 10 and 30 % probability lines.Similarly, the isopach maps and ground measurements for the C and D units of the Askja 1875 produced by Carey et al. (2010) are in agreement with our 10 and 20 % probability lines for ground tephra accumulations of 1 and 10 kg m −2 (Fig. 11).

Atmospheric concentration
Figure 13 summarizes the most likely dispersal trends, in agreement with the wind transect of Fig. 2, and shows that the areas most probably affected by far-range dispersal of ash are Scandinavia and the northern UK.Such results are in agreement with the compilation of the tephrochronological studies of Swindles et al. (2011), who show that the past 7000 years of volcanic activity in Iceland resulted in the identification of 38, 33, and 11 tephra layers in Scandinavia, Ireland and Great Britain, respectively.As suggested by Lacasse (2001), Scandinavia is subject to zonal airflow, whereas Ireland is more likely to be affected than the rest of Europe as it is most probably in the path of anticyclonic airflows from Iceland.As a result, minimum estimates provided by Swindles et al. (2011) show that, based on the record of the past 1000 years, northern Europe is affected by volcanic ash with a mean return interval of 56 ± 9 years and that there is a 16 % probability of tephra fallout every decade based on a Poisson model.
When using a deterministic approach, Fig. 14 shows that amongst the selected eruptions, an eruption of Hekla 1947 is the most likely to produce critical concentrations above the main European airports.Interestingly, scaling Hekla 1947 with the two main phases of Askja 1875 reveals that an erupted volume falling within the boundaries of a low VEI 4 (Hekla 1947) can produce concentrations more than one order of magnitude larger than an eruption of low-medium VEI 5 (Askja 1875).For example, London Heathrow would suffer ash concentrations of 0.68 and 0.025 mg m −3 following eruptions of Hekla 1947 and Askja 1875, respectively (Fig. 14).Similarly, Davies et al. (2010) report that five eruptions of Hekla with volumes varying between 0.18 and 0.33 km 3 produced tephra beds in Norway, Scotland and Finland, more than 1500 km beyond the source.Such observations support the growing idea that the tephra volume of an eruption is not the primary factor controlling the distal dispersal of fine ash and that the TGSD, the nature of the fragmentation process (i.e.dry vs. phreatomagmatic) and the  weather patterns play important roles (Davies et al., 2010;Swindles et al., 2011).The concept of defining critical thresholds, typically 0, 0.2 or 2 mg m −3 depending on the approach adopted, implies that the hazard level is at its maximum once the concentration threshold has been reached.If, for instance, a value of 0.2 mg m −3 is adopted as critical, the shape of the curve above the threshold for the eruption of Hekla 1947 displayed on Fig. 14 does not provide any relevant information as the level of maximum impact is reached.However, for crisis management purposes, the duration during which concentrations are above the threshold becomes critical.In this perspective, the eruption type has a major control on the hazard and the potential associated consequences.For example, concentration plots using a deterministic approach for all flight levels shown in the Supplement illustrate how short-lasting and intense Plinian eruptions result in single peaks of critical concentration typically lasting for a couple of days, reaching up to 12 mg m −3 above the Oslo airport following a Hekla 1947 eruption.In contrast, an eruption of Katla 1918, although not reaching such high levels of critical thresholds, results in more diffuse signals spanning over longer periods of time over single geographic points.When considering a 3-D volume above the European territory representing the airspace and observing the potential disruptions from a management perspective, continuous emission of tephra with a pulsatory regime, though not producing the high concentrations of Plinian eruptions, are potentially able to become more problematic than short-lasting powerful eruptions.

Conclusions
The present work highlights the challenges of achieving a multi-scale hazard assessment from multiple and heterogeneous sources in order to compare and combine outcomes of the most likely range of possible eruptions.For the selected volcanoes, we could define both semi-probabilistic (i.e.stochastic sampling of wind conditions and ESP deterministically fixed) and fully probabilistic (i.e.stochastic sampling of both wind profiles and ESP) eruption scenarios based on the available data.In each case, we developed new algorithms to assist the identification of eruption parameters for both short-and long-lasting eruptions, which help to achieve the sampling of realistic ESP and account for particle aggregation processes.For the atmospheric dispersal of fine ash, a sound deterministic approach demonstrated the different hazards posed by short-and long-lasting eruptions and showed the importance of the potential disruption time over high concentrations.As a result, the outcomes of this work constitute a first step towards an improved management of future volcanic crises, accounting for most critical aspects of both the geological and atmospheric science sides of the problem.The second step toward a sound impact and risk assessment typically involves the identification of the exposed elements and their vulnerability to the stress constituted by ground tephra accumulation and distal atmospheric ash.Such an approach is tackled in a companion paper by Scaini et al. (2014).

S. Biass et al.: Part 1: Hazard assessment
In terms of hazard assessment, we can conclude that: -Eruption scenarios and ESP must be defined using probabilistic strategies based on strong field observations.
-Based on our probabilistic scenarios (e.g.OES, ERS), Askja represents the most hazardous volcano.
-Based on our deterministic scenarios, Hekla is likely to produce the largest atmospheric concentrations of ash but Katla will result in longer disruptions of air traffic.
-At the Icelandic scale, expected accumulations will mainly be a concern for electrical power-lines and agricultural activities (i.e.accumulations of 10 kg m −2 ).
-Our empirical approach to describe aggregation, although simplistic, is a suitable field-based method for computationally heavy probabilistic analysis, which allows the investigation of a wide range of aggregation conditions.
-Results suggest that the erupted tephra volume is not the primary control on the dispersal, whereas the eruptive style (i.e.long-lasting vs. short-lasting) and the TGSD (i.e.fine vs. coarse distributions) might play primary roles.

Figure 3 .
Figure3.Wind analysis at three atmospheric levels along a N-S section from Reykjavík to the UK (Fig.1) from the ERA-Interim Reanalysis database.These wind roses show the probability that the wind will blow in given directions and speed intervals.

Figure 4 .
Figure 4. ESP used for (a-c) the Hekla 2000-type and (d-f) the Hekla 1947-type scenarios (see Table 3for details).(a) and (d): plume height (m a.s.l.);(b) and (d): erupted mass; (c) and (f): total grain-size distribution.Histograms in (c) and (f) show both the fractions of individual particles (light grey) and aggregates (dark grey) generated based on the algorithm explained in the text; original indicates the original grain-size distribution obtained from sieving (i.e.not containing any aggregates).

Figure 5 .
Figure 5. ESP for the Katla LLERS; (a) plume heights considering all occurrences (i.e.all model runs); (b) eruption duration; (c) erupted masses considering all occurrences; (d) erupted masses considering single runs; (e) example of a TGSD of a single run.See caption of Fig. 4 for explanation of symbols.

Figure 6 .
Figure 6.ESP for the Eyjafjallajökull LLOES; (a) 6 h interval measurements of plume height for the period 14 April-21 May 2010(Arason et al., 2011); (b) corresponding MER for the same period based on wind conditions extracted from the ERA-Interim database and the method ofDegruyter and Bonadonna (2012); (c) disaggregated and aggregated TGSD as inferred byBonadonna et al. (2011).See caption of Fig.4for explanation of symbols.

Figure 7 .
Figure 7. Original and aggregated TGSD of the 1875 Askja eruption based on Sparks et al. (1981) for (a) the phreato-Plinian phase Askja C and (b) the dry Plinian phase Askja D. See caption of Fig. 4 for explanation of symbols.

Figure 8 .
Figure 8. Probability maps (%) for ground accumulation for an eruption at Hekla.(a) ERS for a 2000-type eruption, threshold of 1 kg m −2 ; (b) ERS for a 1947-type eruption, threshold of 1 kg m −2 ; (c) ERS for a 1947-type eruption, threshold of 10 kg m −2 ; (d) ERS for a 1947-type eruption, threshold of 100 kg m −2 .Eruption parameters are summarized in Table 3 and Fig. 4.

Figure 13 .
Figure 13.Atmospheric dispersion of tephra for a threshold of 2 mg m −3 for all FL for the eruption scenarios of Hekla ERS 1947-type (a, b, c), Katla LLERS (d, e, f) and Askja OES 1875-type (g, h, i).Maps show (a, d, g) probability maps of exceeding a concentration of 2 mg m −3 , (b, e, h) probability maps of exceeding an arrival time of 24 h for a concentration of 2 mg m −3 and(c, f, i) probability maps of exceeding a persistence time of 12 h for a concentration of 2 mg m −3 .Probability maps for other thresholds and separate FL are available in the Supplement.

Figure 14 .
Figure14.Airborne concentration of ash at FL150 above the airports of London Heathrow (EGLL), Paris Charles de Gaulle (LFPG), Amsterdam Schiphol (EHAM), Frankfurt (EDDF), Oslo Gardermoen (ENGM) and Copenhagen Kastrup (EKCH) resulting from the eruptions of Hekla 1947, Katla 1918, Eyjafjallajökull 2010 and Askja 1875.Each eruption was run for 10 days starting from 14 April 2010 using similar wind conditions that occurred during the Eyjafjallajökull 2010 eruption.Concentrations for other FL can be found in the Supplement.

Table 4 .
Critical thresholds for probability calculation.

Table 5 .
Probabilities, mean arrival time and mean persistence time for a concentration of 2 mg m −3 for all FL combined above the selected airports shown in Fig.1.NR: not reached.