Dispersion modelling of volcanic emissions

UNIVERSITATIS ACTA UPSALIENSIS

UPPSALA

Digital Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology 1433

Dispersion modelling of volcanic emissions

ADAM DINGWELL

ISSN 1651-6214

ISBN 978-91-554-9704-0

Dissertation presented at Uppsala University to be publicly examined in Axel Hambergssalen, Villavägen 16, Uppsala, Thursday, 17 November 2016 at 10:00 for the degree of Doctor of Philosophy. The examination will be conducted in English. Faculty examiner: Professor David Simpson (Chalmers University of Technology, Department of Earth and Space Sciences).

Abstract

Dingwell, A. 2016. Dispersion modelling of volcanic emissions. (Spridningsmodellering av utsläpp från vulkaner). Digital Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology 1433. 53 pp. Uppsala: Acta Universitatis Upsaliensis.

ISBN 978-91-554-9704-0.

Gases and particles released by volcanoes pose a serious hazard to humans and society.

Emissions can be transported over long distances before being reduced to harmless concentrations. Knowing which areas are, or will be, exposed to volcanic emissions is an important part inreducing the impact on human health and society. In this thesis, the dispersion of volcanic emissions is studied using a set of atmospheric models.

The work includes contribution to the development of the Lagrangian Particle Dispersion Model FLEXPART-WRF. Three case studies have been performed, one studying potential ash emissions from potential future eruptions on Iceland, a second covering SO

2

emissions from Mt. Nyiragongo in D.R. Congo, and a third studying the SO

2

emission rate of the Holuhraun eruption (Iceland) in 2014–2015.

The first study covers volcanic ash hazard for air traffic over Europe. Three years of meteorological data are used to repeatedly simulate dispersion from different eruption scenarios.

The simulations are used to study the probability of hazardous concentrations in ash in European airspace. The ash hazard shows a seasonal variation with a higher probability of efficient eastward transport in winter, while summer eruptions pose a more persistent hazard.

In the second study, regional gas exposure around Mt. Nyiragongo is modelled using flux measurements to improve the description of the emission source. Gases are generally transported to the north-west in June–August and to the south-west in December–January. A diurnal variation due to land breeze around lake Kivu contributes to high concentrations of SO

2

along the northern shore during the night. Potentially hazardous concentrations are occasionally reached in populated areas in the region, but mainly during the nights.

The third study uses inverse dispersion modelling to determine the height and emission rates based on traverse measurements of the plume at 80–240 km from the source. The calculated source term yields better agreement with satellite observations compared to commonly used column sources.

The work in this thesis presents improvements in dispersion modelling of volcanic emissions through improved models, more accurate representation of the source terms, and through incorporating new types of measurements into the modelling systems.

Keywords: dispersion modelling, atmospheric, volcano, gas emissions, volcanic ash, FLEXPART, FLEXPART-WRF

Adam Dingwell, Department of Earth Sciences, LUVAL, Villav. 16, Uppsala University, SE-75236 Uppsala, Sweden. ,

© Adam Dingwell 2016 ISSN 1651-6214 ISBN 978-91-554-9704-0

urn:nbn:se:uu:diva-303959 (http://urn.kb.se/resolve?urn=urn:nbn:se:uu:diva-303959)

Akademisk avhandling som för avläggande av teknologie doktorsexamen vid Uppsala univer- sitet kommer att offentligt försvaras i Axel Hambergssalen, Villavägen 16, Uppsala, torsdag 17 november 2016 kl. 10:00. Disputationen sker på engelska. Opponent: Professor Simpson David (Chalmers Tekniska Högskola, Institutionen för Rymd och geovetenskap).

Sammanfattning

Dingwell, A. 2016. Dispersion modelling of volcanic emissions. (Spridningsmodellering av utsläpp från vulkaner). Digital Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology 1433. 53 pp. Uppsala: Acta Universitatis Upsaliensis.

ISBN 978-91-554-9704-0.

Gas- och partikelutsläpp från vulkaner utgör en fara för människor och för vårt samhälle. Ut- släppen kan transporteras över långa avstånd innan de reduceras till oskadliga halter. Att känna till vilka områden som utsätts för, eller kommer utsättas för, utsläppen är ett viktigt verktyg föratt minska påverkan på folkhälsa och samhälle. I avhandlingen studeras spridningen av ut- släpp från vulkanutbrott med hjälp av en uppsättning numeriska atmosfärsmodeller.

Den Lagrangiska Partikelspridningsmodellen FLEXPART-WRF har förbättrats och applice- rats för spridningsmodellering av vulkanutbrott. Tre studier har utförts, en fokuserar på vulka- naska från potentiella framtida utbrott på Island, den andra studerar SO

2

-ustläpp från vulkanen Nyiragongo i Demokratiska Republiken Kongo, och den tredje studerar SO

2

-ustläpp från ut- brottet i Holuhraun (Island) 2014–2015.

Den första studien uppskattar sannolikheten för att vulkanaska från framtida vulkanutbrott på Island ska överskrida de gränsvärden som tillämpas för flygtrafik. Tre år av meteorologisk data används för att simulera spridningen från olika utbrottsscenarier. Sannolikheten för skad- liga halter aska varierar med årstid, med en högre sannolikhet för effektiv transport österut un- der vintermånaderna, sommarutbrott är istället mer benägna att orsaka långvariga problem överspecifika områden.

In den andra studien undersöks spridningen av SO

2

från Nyiragongo över en ettårsperiod.

Flödesmätningar av plymen används för att förbättra källtermen i modellen. Gaserna transpor- teras i regel mot nordväst i juni–augusti och mot sydväst i december–februari En dygnsvariat- ion, kopplad till mesoskaliga processer runt Kivusjön, bidrar till förhöjda halter av SO

2

nattetid längs Kivusjöns norra kust. Potentiellt skadliga halter av SO

2

uppnås av och till i befolkade områden men huvudsakligen nattetid.

Den tredje studien utnyttjar inversmodellering för att avgöra plymhöjd och gasutsläpp base- rat på traversmätningar av plymen runt 80–240 km från utsläppskällan. Den beräknade källter- men resulterar i bättre överensstämmelse mellan modell- och satellitdata jämfört med enklare källtermer.

Arbetet i den här avhandlingen presenterar flertalet förbättringar för spridningsmodellering av vulkanutbrott genom bättre modeller, nogrannare beskrivning av källtermer, och genom nya metoder för tillämpning av olika typer av mätdata.

Nyckelord: Spridningsmodellering, atmosfär, vulkan, gasutsläpp, vulkanaska, FLEXPART, FLEXPART-WRF

Adam Dingwell, Department of Earth Sciences, LUVAL, Villav. 16, Uppsala University, SE-75236 Uppsala, Sweden.

© Adam Dingwell 2016 ISSN 1651-6214 ISBN 978-91-554-9704-0

urn:nbn:se:uu:diva-303959 (http://urn.kb.se/resolve?urn=urn:nbn:se:uu:diva-303959)

The beginning of knowledge is the discovery of something we do not understand

Frank Herbert

List of papers

This thesis is based on the following papers, which are referred to in the text by their Roman numerals.

I Brioude, J., Arnold, D., Stohl, A., Cassiani, M., Morton, D., Seibert, P., Angevine, W., Evan, S., Dingwell, A., Fast, J. D., Easter, R. C., Pisso, I., Burkhart, J. and Wotawa, G. (2013). The Lagrangian particle dispersion model FLEXPART-WRF version 3.1. Geoscientific Model Development, 6(6):1889–1904.

II Dingwell, A. and Rutgersson, A. (2014). Estimating volcanic ash hazard in European airspace. Journal of Volcanology and Geothermal Research, 286:55–66.

III Dingwell, A., Rutgersson, A., Claremar, B., Arellano, A., Mapendano, Y. and Galle, B. (2016). Seasonal and diurnal patterns in the dispersion of SO ₂ from Mt. Nyiragongo. Atmospheric Environment, 132:19–29 IV Dingwell, A., Rutgersson, A., Arellano, S., Galle, B. Using DOAS

traverses and atmospheric modelling to determine plume height and eruption rate of the 2015 Holuhraun eruption. Manuscript.

Reprints were made with permission from the publishers.

In Paper I, I participated in development of the model, including the design and implemention of a new output format; I further improved on the model in Papers II and III.

I came up with the initial ideas for Papers II, III and IV and drafted and implemented the methods. I configured and executed all model simulations.

All co-authors have contributed through advice and suggestions for both data analysis and writing of the papers.

In Paper II, I had the main responsibility for designing and implementing the

method and data analysis and for the writing of the paper. In Paper III, I

was responsible for most of the data analysis — except for the processing of

the flux-data — and also implemented 2 new land-use options in FLEXPART-

WRF. I had the main responsibility for writing the paper. In Paper IV, I

designed the study and carried out the short- and long-range modelling. I had the main responsibility for the analysis and for writing the paper.

The following related publication is not included in the thesis:

Arellano, S.; Yalire, M.; Galle, B.; Bobrowski, N.; Dingwell, A.;

Johansson, M. and Norman, P. (2016). Long-term monitoring of

SO 2 quiescent degassing from Nyiragongo’s lava lake. Journal

of African Earth Sciences.

Contents

1 Introduction

. . . .

11

1.1 Aim of this thesis

. . .

12

2 Modelling tools

. . .

13

2.1 Meteorological data

. . .

13

2.1.1 Downscaling meteorological data

. . .

14

2.1.2 Meteorological model

. . .

15

2.2 Dispersion modelling

. . . .

16

2.2.1 Advection of particles

. . .

17

2.2.2 Deposition processes

. . .

18

2.2.3 Aggregation and Chemistry

. . .

20

2.2.4 Plume rise from volcanic eruptions

. . . .

21

2.2.5 Mass eruption rates

. . . .

22

2.2.6 Inverse modelling

. . .

23

3 Model application

. . . .

26

3.1 Ash hazard climatology

. . .

26

3.2 Regional exposure to passive degassing

. . .

28

3.3 Inverse modelling of SO 2 emissions

. . . .

29

3.4 Model enhancements

. . .

30

4 Results

. . . .

32

4.1 Ash hazard climatology

. . .

32

4.2 Regional exposure to volcanic degassing

. . .

34

4.3 Inverse modelling of SO 2 emissions

. . . .

36

5 Concluding remarks

. . .

39

5.1 Outlook

. . .

40

6 Acknowledgments

. . .

42

7 Sammanfattning p˚a svenska (Summary in Swedish)

. . .

45

References

. . . .

49

1. Introduction

Planet Earth is host to many volcanic systems, about 1,300 of these have had eruptions within the past 10,000 years. Some eruptions are short lasting, such as the frequent eruptions at Sakurajima (Japan) or Vulcano (Italy); others can last for several decades or millennia, such as Stromboli volcano in Italy. Volca- noes in general are natural emitters of aerosols and numerous hazardous gases.

These emissions can have a negative impact on human health, infrastructure and/or the environment. Emission plumes from large eruptions can maintain harmful concentrations for days or even weeks, before being reduced through removal processes and dispersion. Gas emissions can occur either as part of eruptions or through passive degassing over longer periods. The emissions can affect cloud chemistry and microphysics and cause acid rain (e.g. Delmelle et al., 2001; Rose et al., 2006). Exposure to gas plumes can increase wear on man-made structures (e.g. antennas, machinery, roofs). Acidic compounds (e.g. SO 2 , H 2 S, HCl, HF) can also have an irritating effect on airways when inhaled (Delmelle et al., 2002). Fluoride (F) — which is mainly emitted as HF — is toxic at high concentrations and can disturb dental development for children (Baxter et al., 1999; Delmelle et al., 2002); it can even prove fatal for grazing animals, which might ingest large quantities of fluoride coated plants (Cronin et al., 2000; Thorarinsson and Sigvaldason, 1972).

Volcanic emissions have been dispersed over Europe many times through- out history, one of the most notable cases was the Laki fissure eruption in Iceland in 1783–1784. The eruption had severe impacts on Iceland; most of Iceland’s livestock died within a year from the eruption by ingestion of ash coated grass. During the this period, the human population suffered from famine and exposure to the toxic emissions; the population was reduced by 19–22 %, or approximately 22,000 people (Demar´ee and Ogilvie, 2001). The eruption also affected the European continent, where a persistent haze was re- ported in mid to late June 1783, lasting for about two months and eventually extending as far as Moscow, Baghdad and northern Africa (Stothers, 1996).

The eruption received attention recently, when a fissure eruption started in Holuhraun (Iceland) on the 29th of August, 2014; the eruption lasted until late February 2015. The Holuhraun eruption was one of the largest emissions of SO 2 since the Laki eruption.

Gas emissions from volcanoes are not uncommon, the main contribution

being from passive degassing (Andres and Kasgnoc, 1998). Individual cases

of passive degassing generally have weaker emission rates than events like

the Laki eruption, however, they can still pose a regional hazard. Passively

degassing volcanoes often entails a negative impact on human health and the environment in the downwind areas (Hansell and Oppenheimer, 2004).

Explosive eruptions pose a different hazard by releasing ash at high alti- tudes; the ash can be transported long distances and cause serious damage to jet engines. The eruption of Eyjafjallaj¨okull in April 2010 resulted in almost complete closing of the European airspace through pre-cautionary decisions.

Fine ash was injected in the jet stream and transported over 1000 km, reaching continental Europe within a matter of days. Around 104,000 flights were can- celled form the 15th to the 22nd of April (EUROCONTROL, 2010). While the total ash emissions of the eruption were not exceptional, the fraction of fine ash produced was unusually high. Persistent northerly winds resulted in a large portion of the fine ash to reach continental Europe (Stohl et al., 2011).

In 2010, the European countries were not sufficiently prepared to efficiently manage such a situation. As a response, new guide-lines for volcanic ash were decided upon and training programmes were started.

The impact of extreme natural events is highly dependent on the informa- tion available to decision makers. Increased knowledge and better predictions enable our society to make preparations and better manage the situation when eruptions do occur.

1.1 Aim of this thesis

In this thesis, exposure to volcanic ash and gases are studied and techniques for estimating hazards for different types of volcanoes are further developed.

• Paper I presents the first stable release of the multi-scale Lagrangian particle dispersion model FLEXPART-WRF.

• The probability of encountering volcanic ash in European airspace as the result of potential future eruptions on Iceland is investigated in Paper II.

• Regional exposure to gases emitted from Mt. Nyiragongo (eastern Demo- cratic Republic of Congo) is covered in Paper III.

• Different methods for calculating fluxes and plume height are presented

and evaluated in Paper IV.

2. Modelling tools

The transport of pollutants in the atmosphere can be affected by a variety of processes (advection, turbulent mixing, fallout, impaction, rain-out, wash-out, radioactive decay, and/or chemical reactions); a full understanding of the com- bined effect of these processes is hard to achieve even if each process in itself is well understood. Atmospheric models in general are used in order to study the combined effect of various physical and/or chemical process.

A combination of two types of models have been used throughout this the- sis: a meteorological model, mainly used to produce wind and precipitation data, and a dispersion model, which uses the meteorological data to simulate transport and removal of the emissions. In Papers II and III the Weather Re- search and Forecasting (WRF) model was used together with the Lagrangian Particle Dispersion Model (LPDM) FLEXPART-WRF. In Paper IV simula- tions made with FLEXPART-WRF were compared to results from the similar FLEXPART model.

2.1 Meteorological data

Dispersion modelling requires certain meteorological data, mainly winds and precipitation. For short scales, single meteorological measurements can be used directly to force a dispersion model. For application over larger areas gridded data is required, which best achieved through meteorological mod- els. One option is to use a combination of analysis and forecast data, which is continuously produced by a number of meteorological centres around the world. The European Centre for Medium-range Weather Forecasts (ECMWF) produce global, hourly forecasts using their Integrated Forecast System (IFS).

Today, the global forecasts can be retrieved by member states down to a res- olution of 0.1 degrees ( ∼ 10 km), however, this is a fairly new product and older data have a coarser spatial resolution. The global forecasts are suitable for use in atmospheric dispersion model as long as the temporal resolution is sufficient.

If higher resolution is required it might be a better option to use reanaly-

sis data, e.g. when studying short range dispersion or during a period when

high resolution forecast data is not available. Reanalysis data is produced in a

similar fashion as forecast data but is to a large extent produced retroactively

(Warner, 2011). A reanalysis product is produced using a large set of observa-

tions and a meteorological model; data is typically produced covering several

decades to provide a long time series of consistent data. This means that re- analysis data seldom uses the latest versions of models and usually run at a coarser resolution compared to the most recent forecasts.

ERA-Interim (Dee et al., 2011) is a reanalysis product produced by ECMWF which has been used extensively in this thesis. The data comes on a horizontal resolution of 0.75 degrees with an output interval of 6 hours. The large gaps between output times can be problematic in dispersion modelling, e.g. due to meandering of winds or propagation of frontal systems which might appear to skip over grid cells if the temporal resolution is too coarse. Furthermore, a higher spatial resolution is sometimes desired, in which case an additional step is necessary to increase the resolution. An efficient way to solve most of these problems is dynamical downscaling.

2.1.1 Downscaling meteorological data

In order to improve the resolution of a meteorological dataset, one can either interpolate (which only smooths the data without adding much information), use a statistical downscaling (which is more suitable for climatology) or one can use a limited-area meteorological model (LAM). The coarse data is used to initiate the LAM, which itself runs at a higher resolution. Forcing is applied along the boundaries of the model domain, throughout the simulation, using data from the global model. In this way, the large-scale dynamics are used to feed smaller scale dynamics in the LAM which cannot be resolved by the global model. The LAM is able to use more detailed land-use and topography data than a global model; both topography and land-use data will influence the development of the smaller-scale dynamics.

The method of increasing resolution by the use of a LAM, is called dynam- ical downscaling. The result is a more detailed description of the meteoro- logical situation than what can be achieved by interpolation. The downside is that dynamical downscaling requires far more computational resources than other downscaling techniques. A comparison between interpolated and dy- namically downscaled meteorological data is shown in Figure 2.1; both cases produce smooth solutions but the interpolated data lacks the additional infor- mation which is added with dynamical downscaling.

The goal of dynamical downscaling is to add information about smaller-

scale processes while retaining the information from the large scale processes

in the global model. Downscaling should be made incrementally if there is

a large difference between the initial resolution of the global model and the

final resolution of the LAM; this is due to how forcing is applied during a

simulation. If there is a large resolution-gap between parent (low resolution)

and child (high resolution) domains, then processes from intermediate scales

will never be resolved. This leaves the child domain without proper forcing,

allowing it to develop features which, in reality, do not exist.

245 250 255 260 265 270 275 280 Temperature (K)

66°N 64°N 62°N 60°N

58°N

40°W 35°W 30°W 25°W 20°W 15°W 66°N 64°N 62°N 60°N

58°N

40°W 35°W 30°W 25°W 20°W 15°W

Figure 2.1. Comparison of near-surface temperature from basic interpolation (left) versus dynamical downscaling (right) from Dingwell (2012). The left image was produced by interpolating data from ERA-Interim at 0.75 degree resolution ( ≈ 35 × 85 km) down to 13 .5 × 13.5 km. The right image has the same resolution, but was produced by the WRF-model, using the interpolated data from ERA-Interim as forc- ing.

Incremental nesting is especially important if there are e.g. mountains, lakes or islands nearby which are too small to be properly resolved by the global model but which have an important influence on the dynamics in the region.

With nesting, these influences can be accounted for in an intermediate domain which might not require much computational time compared to expanding the innermost domain. Furthermore, numerical problems may arise if there is a too large difference between the input data and the computational grid (Warner, 2011). There are however cases where direct nesting can be just as accurate as incremental nesting, e.g. climatology studies (Beck et al., 2004).

We have used incremental downscaling with resolutions increasing by a factor of 3 between nests; this is the most common approach with the meteorological model we used.

2.1.2 Meteorological model

When downscaling meteorological data, we have used the Advanced Weather Research and Forecasting (WRF) model (Skamarock et al., 2008). The WRF- model is a 3-dimensional, non-hydrostatic, mesoscale meteorological model.

It is designed with a modular approach, which enables wide use-cases from operational forecasting to idealized scenarios. Multiple options are available for most physics schemes in the model, several of which were tested for Pa- pers II, III and IV.

For the purpose of this thesis, WRF has been set up to downscale ERA-

Interim reanalysis data. First, the model is initiated in each grid-cell using the

interpolated reanalysis data. The model then requires 6–12 simulation hours to

relax the dynamics and precipitation fields (spin-up) before producing reliable results. This type of initialization is called a cold start. During a simulation, the reanalysis data is used to force the model at the boundaries of the outermost model domain. No forcing is applied within the domains; this might cause the model to deviate from the forcing data during longer simulations in a similar way as when a long-term forecast deviates from the actual weather. Sequential initiations were used in order to prevent the model from deviating too far from its forcing data. A restart interval of 6 + 24 hours (spin-up + production) was used for the mid-latitude simulations in Paper II. A longer spin-up time was required in Paper III; therefore, the restart interval was extended to 12 + 48 hours. There are several possible reasons why a longer spin-up was required.

Paper III covered the inland tropical region of eastern D.R. Congo. Numer- ical weather prediction is more difficult in the tropics compared to the mid- latitudes, where predictable synoptic systems mostly dominate. In the tropics, mesoscale processes have a larger influence on the weather (Laing and Evans, 2011) but the forcing data are too coarse to fully resolve them. Furthermore, the model resolution was increased in Paper III; with a large difference be- tween the forcing data and the final product, more work could be required by WRF to “fill in the gaps”.

The WRF model simulates soil moisture and temperature down to 200 cm, divided into 4 layers. Soil moisture has a strong influence on surface fluxes in mesoscale models (Angevine et al., 2014). Soil layers need more time for spin-up than the atmosphere; the topmost soil layer can spin up in a couple of days but the deeper layers take longer. Therefore, in Paper III, the soil layers were retained between restarts allowing them to develop throughout the entire simulation period. Similarly, accumulated fields (e.g. precipitation) were retained between restarts (in Papers II and III), producing continuously increasing values; this was not necessary in Paper IV since the simulation periods were short enough to cover in a single model run.

2.2 Dispersion modelling

Atmospheric dispersion models are used to simulate the transport of emissions in the atmosphere; some models are designed to run in parallel with a meteo- rological model (online) while others run separately using prepared meteoro- logical data (offline). Online models have the benefit of allowing emissions to affect the dynamics in the meteorological model (active trace feedback); how- ever, these models are more complicated to design and tend to require more computational resources than offline models.

The two most common types of dispersion models are Eulerian grid models

and Lagrangian Particle Dispersion Models (LPDM). Eulerian models have a

fixed grid through which emissions are transported by calculating fluxes be-

tween adjacent grid cells. Calculations for emissions, removal processes and

chemistry are applied to each grid cell. A Lagrangian model works differently, emissions can be assigned anywhere within the model domain and are not lim- ited by a superimposed grid. The emissions are assigned to computational par- ticles which are released within a point-, line-, area-, or volume source. The computational particles represent air parcels containing a mixture of chemical species or particulate matter. The computational particles are moved within the model domain through advection and turbulent mixing. Computational particles simulating clouds of aerosols are also affected by gravitational set- tling. Calculations are made on a per-particle basis rather than for a set of predetermined grid cells; this allows LPDMs to efficiently simulate emissions from a single source in a large model domain without wasting calculations on empty grid cells. Another benefit from LPDMs is their ability to generate re- alistic thin plumes from point sources; Eulerian models suffer from numerical diffusion which causes plumes to “fan out” excessively unless compensated.

LPDMs are not perfect, like Eulerian models they have certain limitations.

An LPDM requires a large number of computational particles to be released;

if the distance between particles become to large, the output becomes patchy and harder to interpret reliably. The distance between particles will increase further from the source as the plume fans out due to mixing processes. For large domains, one must either release an excessive amount of particles at the source or add particles along the way as the plume fans out. It is also harder to implement chemistry and aggregation of particulates in a Lagrangian model;

chemical reactions are typically limited to reactions between species within a particle or between emitted species and background species.

In this thesis, two LPDMs were used, FLEXPART-WRF (Papers II, III and IV) and FLEXPART (Paper IV). The main difference between the mod- els is the meteorological data they require; FLEXPART-WRF uses meteoro- logical data from WRF, while FLEXPART can run on different data products from ECMWF. FLEXPART-WRF was developed based on the FLEXPART model, technical aspects of its first stable release is described in Paper I.

2.2.1 Advection of particles

The core of any LPDM is the advection of computational particles. Particle positions are updated by moving particles according to the wind field, i.e. by numerically solving a differential equation for each particle’s position. The simplest advection scheme is to use the Euler method — not to be confused with Eulerian models — as in FLEXPART and its derivatives:

¯r (t + Δt) = ¯r(t) + ¯V(¯r,t)Δt (2.1)

where ¯r (t) is the particle’s position at any given time, t; Δt is the time step of

the simulation and ¯r (t + Δt) is the particle’s position at the next step. Some

models use more elaborate methods, e.g. the HYSPLIT model which uses

the improved Euler’s method (also called Heun’s method) for calculating the advection (Draxler and Hess, 1998). Additional terms are typically added to (2.1), representing horizontal and vertical turbulent perturbations and/or sub- grid convection, which displace particles from the mean flow.

When modelling aerosols, the dispersion model needs to account for gravi- tational settling by adding a term for terminal velocity, v _g , to (2.1). For spher- ical particles larger than ∼10 μm, v g can be calculated as a balance between downward gravitational pull and the upward drag force:

v g =

 4 3

d _p g C _D

 ρ − ρ air

ρ air



≈

 4 3

d _p g C _D

 ρ ρ air



(2.2) where d _p is the particle’s diameter, g is acceleration due gravity, ρ is the par- ticle density and ρ air is the density of air. C _D is the drag coefficient which, in FLEXPART and its derivatives, is calculated using the method by N¨aslund and Thaning (1991). The drag coefficient varies with Reynolds number which in turn depends on the velocity, therefore, (2.2) and C _D must be determined iteratively.

Smaller particles ( ∼1 μm) are affected by the mean free path between mole- cules in the air. In this case, the air behaves less like a continuum and there- fore, the drag force exerted on a small particle will be weaker than what is assumed in (2.2). In order to accurately describe terminal velocities of parti- cles down to micrometre size, (2.2) can be adjusted by dividing C _D with the Cunningham correction factor as proposed by Cunningham (1910); Knudsen and Weber (1911):

C _cun = 1 + 2 λ d _p



A 0 + Q · e ^−C

^2dpλ



(2.3) where λ is the mean free path of gas molecules in air. A 0 , Q and C are di- mensionless constants, which in FLEXPART-WRF are set to 1.257, 0.400 and 1.10, respectively. While (2.2) and (2.3) enable the calculation of settling velocities of a wide size-range of particles, there are still limitations. The ex- pression (2.3) assumes spherical particles, which is rarely the case. Correcting for non-spherical particles is usually done by defining particles by their aero- dynamic diameter, i.e. the diameter which a spherical particle should have in order to have the same terminal velocity as the studied particle.

2.2.2 Deposition processes

Deposition of aerosols or gases from the atmosphere can take two paths: dry

or wet deposition. Dry deposition includes gravitational settling to the surface

(large particles), impaction with the surface (small particles), and adsorption

or absorption of gases by various surface elements. Wet deposition includes

the nucleation of rain droplets within a cloud (rain-out) and the collection of particles or gases through collision with existing rain drops (wash-out).

Dry deposition can be described in terms of a deposition velocity, V _d , de- fined as the flux, F divided by the concentration, C, of a given species. Dif- ferent processes apply to particles and gases, so the two cases need to be de- scribed separately.

Dry deposition of gases can be determined using resistances similar to an electric circuit (McMahon and Denison, 1979; Wesely and Hicks, 2000):

V _d = 1

r _A + r B + r C

(2.4)

where r _A is the aerodynamic resistance of the surface layer which can be esti- mated from Monin-Obukhov similarity theroy (Stull, 1988). r _B is the viscous resistance of the quasi-laminar sublayer, and r _C is the bulk surface resistance which collects a number o off surface-specific properties. An illustration of the processes described in (2.4) is shown in Figure 2.2. The most common ap- proach to determining r _C is by using a number of land-use categories with pre-determined resistances for a number of known chemical species. We- sely (1989) proposed a parameterization for 11 land-use types and 2 chemical species (O 3 and SO 2 ); it was suggested that dry deposition for other species could be parameterized as well by comparing their oxidizing properties with O 3 and their solubility with SO 2 .

A

B

C

net flux

Figure 2.2. Illustration of the three sub-regions for dry deposition. (A) represents the surface layer, where the resistance is determined by turbulent eddies. (B) is the laminar sublayer, a thin layer of air in contact with surfaces, where turbulence is negligible and viscous forces dominate. Transport through this layer depend on Brownian motion and inertia of particles or molecular diffusion of gases. (C) represents the bulk surface properties, which combines numerous processes for absorption of gases by the surface.

Dry deposition of particles require a different expression than what is used for gases. Surface resistances are negligible for particles, which tend to “stick”

to surfaces due to their own weight; instead, settling velocities need to be taken

into account as a parallel process to dynamic and diffusive fluxes. This can be

described by (e.g. Seinfeld and Pandis, 2016):

v _d = 1

r A + r B + r A r B v g + v g (2.5) which is analogous to (2.4). The third term in the denominator is not a physical resistance, but appears when deriving the resistances from fluxes and concen- trations. In LPDMs the dry deposition is typically implemented as a mass loss per computational particle, using a decay function:

m(t + Δt) = m(t)exp

 −V d (z)Δt 2z



(2.6) where z is the altitude above ground level and m (t) is the particle mass at a given time. (2.6) is typically applied to particles below a given altitude, e.g.

within an assumed constant flux layer.

Wet deposition can be expressed similarly, as was suggested by McMahon and Denison (1979):

m (t + Δt) = m(t)exp(−ΛΔt) (2.7)

where Λ is a scavenging coefficient which depends on the precipitation inten- sity. Empirical tests by McMahon and Denison (1979) suggested the following from for the scavenging coefficient:

Λ = AI ^B (2.8)

where A and B are coefficients which need to be determined for each modelled species (i.e. for each particle type and for each gas). Equation (2.8) is designed for below cloud scavenging, i.e. for the collection of particles/gases by falling rain drops.

Scavenging can also be the result of nucleation of cloud droplets on parti- cles or the dissolution of gases in existing cloud droplets; both of these pro- cesses contribute to in-cloud scavenging. In-cloud scavenging can occur either prior to or simultaneously as precipitation. The actual removal of species from the atmosphere is through precipitation, regardless of the scavenging process.

Hertel et al. (1995) described a parameterization which includes in-cloud scav- enging according to:

Λ i = S i I

H _i (2.9)

where S _i is the scavenging ratio and H _i is the height over which scavenging

takes place. The scavenging ratio is the ratio between concentrations in pre-

cipitation and concentrations in air; it is determined differently for particles

and gases. For particles, S _i is determined by the fraction of particles activated

as condensation nuclei. For gases, S _i is determined from the solubility of a

given gas (based on Henry’s law). For details, see Hertel et al. (1995).

2.2.3 Aggregation and Chemistry

Chemistry and particle aggregation might affect the removal rate and aero- dynamic properties of emissions; neither of the two are easy to include in Lagrangian models. The most common solutions are to superimpose an Eu- lerian grid for chemistry and aggregation calculations (as in NAME) or re- strict chemistry to reaction between primary species and climatological back- ground species (as in FLEXPART). Another example is the GRAL-C model, which combines Lagrangian and Eulerian frameworks in a way where primary species follow Lagrangian dispersion while secondary species are dispersed through an Eulerian grid (Oettl and Uhrner, 2011). Businger et al. (2015) demonstrated a fully Lagrangian model where SO ₂ is converted to sulphate particles over time; both species were represented in each computational parti- cle, therefore, no gravitational settling could be applied. The fully Lagrangian approach to chemistry works as long as there is only a single (or a few) simple reactions involved.

Volcanic plumes contain numerous reactive compounds in gas, liquid, and solid phases; they contain large amounts of sulphuric gases which, given time, will form droplets of sulphuric acid. Furthermore, particles in a volcanic plume might adhere to one-another and form aggregate particles with differ- ent aerodynamic properties. Particles released directly into the atmosphere are called primary particles. The term secondary particles is used to refer to par- ticles formed from gas-to-particle formation, aggregation of other particles, or condensation of gases on existing particles; “secondary” since they form from the primary emissions. The formation of secondary particles will alter wet and dry deposition rates over time.

One option for handling secondary particles is to adjust the primary emis- sions (at the source) to emulate secondary particle formation which, in reality, occurs further downwind. This approach is not perfect since the physical and chemical properties of the emulated secondary particles will apply directly at the source, instead of gradually building up further downwind.

Furthermore, current parameterizations for atmospheric chemistry are mainly designed for industrial emissions. Volcanic plumes contain different mixtures of emissions than emissions from industries, therefore the parameterizations used in atmospheric modelling might not be a good representation of the pro- cesses within volcanic plumes.

2.2.4 Plume rise from volcanic eruptions

Volcanic eruptions can be intense events; large eruptions can transport mate- rials far into the stratosphere, through a combination of plume rise processes.

Sparks and Wilson (1976) presented a general model for plume rise from vol-

canoes, where the eruption column was divided into two parts: the gas thrust

and convective regions. In the vent and conduit of an erupting volcano, ma-

terials can be accelerated to speeds greater than 100 ms ⁻¹ due to the rapid expansion of gases (Blackburn et al., 1976; Wilson et al., 1980). As the ma- terials are released into the atmosphere, they continue to flow upward due to inertia.

The lowermost part of the eruption column is non-buoyant so the ascent velocity decreases with altitude. This part of the eruption columns makes up the gas-thrust region. Within the gas-thrust region, heavier particles fall out, effectively decreasing the bulk density of the column. Ambient air is entrained into the eruption column which is heated by the erupted material; the heated air expands, which further decreases bulk density and increases buoyancy (Sparks and Wilson, 1976; Sparks, 1986).

If enough air is entrained into the eruption column, it becomes buoyant and continues to rise, forming the convective region. The convective region typically makes up the major contribution to the plume rise when present.

The extent of the convective region has been shown to depend on how much thermal energy is released at the vent, ambient air temperature and lapse rate (Wilson et al., 1978; Settle, 1978) as well as ambient wind conditions (Graf et al., 1999).

Eventually, the eruption column will reach an altitude where its density equals that of the ambient air, buoyancy becomes neutral but the material con- tinues to rise due to inertia. This led Sparks (1986) to introduce a third compo- nent of the eruption column: the umbrella region. In the umbrella region, the upward flux is inertia driven; the column reaches beyond the level of neutral buoyancy and begins to expand laterally due to gravity flow (i.e. by displacing the now lighter ambient air). Figure 2.3 shows an illustration of an eruption column with three components, with corresponding temperature and velocity profiles.

gas thrust convective umbrella

velocity

height

atmosphere

density plume

height

Figure 2.3. Qualitative illustration of an eruption column. The column is divided into three regions with different velocity and buoyancy characteristics. (Original image:

Sparks (1986))

Not all eruptions will produce a three-component eruption column, e.g.

Wilson et al. (1978) found that a too large eruption rate could prevent a con-

vective component from forming. For weak eruptions, or passive degassing,

the column can consist solely of a convective component.

Volcanic eruption columns are small scale processes when compared to meteorological modelling, including eruption dynamics in a meteorological model requires a high resolution which is computationally expensive. There are models designed to bridge the two areas, such as the Active Tracer High resolution Atmospheric Model, ATHAM, (Textor et al., 2006). However, in broader application of dispersion modelling it is usually not possible to include a fluid-dynamic representation of the source due to the added computational demand.

2.2.5 Mass eruption rates

While the plume rise from volcanic eruptions can often be observed directly, determining the mass eruption rate (MER) is usually more difficult. Building upon the model of a convective eruption column, there should be a relation between the heat flux at the vent and the total plume rise. Assuming an even distribution of heat, the MER should be proportional to the heat release. Sev- eral empirical relations have been proposed to describe the relation between eruption rate and plume rise. Analogous to plume-rise from smoke stacks, Settle (1978) proposed the following relation based on data from 6 historical eruptions:

ΔH = 0.117

 dM dt

 0 .22

(2.10) where dM /dt is the MER and ΔH is the plume height above the vent. The MER was calculated based on deposit volumes and pyroclast densities, the eruption heights were based on eye-witness accounts. A more recent attempt was made by Mastin et al. (2009), who proposed the following relation based on data from 34 historic eruptions:

ΔH = 2.00

 dV dt

 0 .241

(2.11) where dV /dt is the volume eruption rate of solid and liquid material. Since the mass release should be dominated by liquid and solid particles, dV /dt should be proportional to dM /dt. In reality, eruption dynamics are far more complicated, with possibly rapid changes in MER and the composition of the emissions (e.g. shapes and sizes of particles or mixture of gases) and with both temporal and spatial variations in atmospheric properties (e.g. stability, wind shear and/or humidity). Modellers rarely have the possibility to include but a few of these processes.

2.2.6 Inverse modelling

Designing eruption parameters based on empirical relations introduces large

uncertainties to the dispersion modelling. An alternative method is to use

inverse modelling to determine the MER and emission height. Eckhardt et al.

(2008) demonstrated how inversion schemes could be used to improve results from dispersion modelling when combined with satellite observations of the plume. The method involves setting up a number of potential sources. A dispersion simulation made with the potential sources yields an a-priori plume which can then be compared with measurements.

The aim of the a-priori plume is not to replicate the true plume but to present the modeller with the relation between gridded concentration data and emis- sion rates from each of the potential sources. By comparing the a-priori results with e.g. satellite observations, it is possible to optimize the source in order to create a posteriori plume which better resembles the observation. This is done using a regression algorithm which yields new emission rates for the potential sources. This is called an inversion scheme, since the model output is used to determine the source parameters.

Inverse modelling is a type of regression problem, i.e. the goal is to find a solution for:

y≈Mx (2.12)

where x is a n × 1 matrix of source elements, y is a m × 1 matrix of observa- tions (e.g. column densities of a SO 2 ), and M is an m × n matrix calculated using the dispersion model, determining the relation between x and y. In other words, (2.12) is an equation system with m equations and n unknown variables (x). If n > m, the problem is overdetermined and needs to be solved through regression. One method of solving (2.12) for x is to use a linear-least squares solution, which changes the problem to:

x = 

M ^T M  ₋₁

M ^T y (2.13)

The least-squares solution cannot be directly applied to dispersion modelling for several reasons. First, the solution allows negative emission sources which is unphysical. Second, a least-squares solution can provide a reliable result as long as M is well conditioned and correlated errors are small; this is often not the case for dispersion modelling. For these reasons, a more robust method is needed.

Eckhardt et al. (2008) developed an inversion scheme suitable for use with volcanic emissions and satellite observations. The potential emission sources (x) were set to represent different levels in a column source. A-priori emission strengths (x _a ) were assumed to have an influence on the final result. Including a-priori emission in the inversion required a modified form of (2.12) which included both a-priori and posteriori source terms:

M˜x ≈ ˜y (2.14)

where ˜x = x − x a is the difference between posteriori and a-priori source terms,

and ˜y = y − Mx a is the difference between observations and a-priori modelled

values. The least-squares solution (2.13) was updated with the new terms (J ₁ ) and expanded to include two terms, J ₂ and J ₃ , representing the deviation from a-priori source terms and the smoothness of the source terms, respectively:

J ₁ = σ _o ⁻² (M˜x − ˜y) ^T (M˜x − ˜y) J ₂ = ˜x ^T diag (σ x −2

)˜x (2.15)

J ₃ = ε (D˜x) ^T D˜x

where σ o is the observation error (assumed to be the same for all observations, diag(σ _x ⁻² ) is a diagonal matrix with the main diagonal consisting of squared inverse errors of the a-priori emissions, D is the discrete representation of the second derivative and ε is a regularization parameter which is used to adjust the weight of the smoothness term. The inversion was completed by solving (2.14) for ˜x. Negative values in x were avoided by decreasing the values of corresponding elements in σ x (nudging toward a-priori values) and iterating until only negligible negative values remain.

Throughout the rest of this thesis, this method will be referred to as the Eckhardt inversion scheme. The Eckhardt inversion scheme is specifically designed to work with dispersion modelling of volcanic emissions and is more robust than simpler methods such as ordinary least-squares or non-negative least squares. However, the method requires some tuning before each use case. It is possible to apply (2.15) to different configurations, e.g. Stohl et al.

(2009) used it to determine the contribution to global emission of greenhouse

gases from different regions around the world.

3. Model application

During the development of the stable release of FLEXPART-WRF (Paper I), several validation studies were made (e.g. Fast and Easter, 2006; Brioude et al., 2012; Angevine et al., 2013). The model has since been used in various cases. For the purpose of this thesis the model has been applied to study dispersion from volcanic emissions for three different cases.

3.1 Ash hazard climatology

In Paper II the hazard from various eruption scenarios was studied. The prob- ability of exposure to airborne ash was estimated over the north Atlantic Ocean and Europe by modelling the dispersion of emissions from different eruption types. A total number of 90,000 dispersion simulations were performed cov- ering various weather conditions over a period of 3 years. The study required realistic representation of the eruption parameters in each scenario. Basic eruption parameters — such as plume height, eruption rate, and grain-size distribution — were retrieved from records of historic eruptions (Table 3.1).

The goal was not to completely replicate any of the eruptions, but rather to create a set of realistic parameters for each scenario. FLEXPART-WRF was used for the dispersion modelling; the model can only simulate particles of a single size at a time so separate simulations were performed for the different particle sizes. The particle-specific simulations were then combined and the full plume was analyzed.

The evolution of the plume from each simulated eruption was studied for 96 hours from the onset of the eruption. Hourly averages of the modelled concentrations were produced and compared with threshold values used to determine flying conditions for commercial aircrafts. Maximum hourly con- centrations were calculated for three different altitude ranges (flight levels in units of 100 ft above sea level: 0–FL200, FL200–FL350, FL350–FL550). The probability of exceeding each threshold value was calculated for each altitude range on a 15 ×15 km grid over northern Europe (see Figure 3.1A).

An important part in setting up the study was assigning the vertical dis-

tribution of emissions. As was discussed in Chapter 2, plume height is of-

ten governed by the convective updraft in the eruption column. In dispersion

modelling, a common approach has been to assign a uniform column source

representing the whole eruption column (e.g. Witham et al., 2007). Using

a uniform column source introduces some errors into the model. First there

Table 3.1. Eruption parameters used for the scenarios simulated in Paper II. Each scenario is based on a historic eruption, some of which with multiple explosive phases;

for those cases, all phases were included in each simulation. The columns present eruption column height (ΔH), duration (dt), mass eruption rate ( ^ΔM _Δt ) as well as the fraction of the mass ( ×10 ⁻² ) allocated to the different particle classes (6Φ − 11Φ).

Table from: Paper II.

Case Phase ΔH Δt

^ΔM_Δt

6 .0Φ 7 .5Φ 9 .0Φ 11 .0Φ VEI

[km] [h] [kgs

⁻¹

] 16 μm 5.5 μm 2.0 μm 0.49 μm

Mt.Spurr 1992 1 14.0 3.5 3 × 10

⁶

13 5.3 1.1 0.05 3

1 8.0 6.5 2 × 10

⁶

20 11 13 1.1 4

Eyjafjallaj¨okull 2010 2 5.5 3.5 7 × 10

⁴

3 6.0 5.4 6 × 10

⁵

1 8.0 1.0 3 × 10

⁶

20 11 13 1.1 5

Askja 1875 2 22.8 1.0 1 × 10

⁸

3 26.0 6.0 2 × 10

⁷

Pinatubo 1991 1 35.0 3.0 4 × 10

⁸

9.7 2.6 0.39 0.26 6

2 30.0 6.0 2 × 10

⁸

Tambora 1815 1 35.0 28.0 5 × 10

⁸

30 9.5 3.5 0.90 7

is no accounting for adiabatic expansion of the column as the gas mixture is de-pressurized; the concentration of emissions should actually decrease with altitude. Second, the majority of the fine ash should be released at the upper portion of the eruption column due to it’s lower terminal velocity (Carey and Sparks, 1986; Koyaguchi and Ohno, 2001).

~

Lake Kivu 60°N

50°N

40°N

30°N

(A) (B)

Kitchanga KaheMasisi SakeGoma E Goma W Nyamulagira (3058 m) Nyiragongo (3470 m)

Figure 3.1. Geographical coverage of the two studies included in this thesis. (A) shows three nested domain, all of which were used in the dispersion simulations in Paper II; a similar setup was used in Paper IV but focusing mainly on Iceland. (B) shows the single domain (third order nest) used in the dispersion simulations in Pa- per III.

Some problems from uniform column sources can be overcome by dis- tributing mass differently within the column, making an non-uniform column source. This was done by e.g. Peterson and Dean (2008); Steensen et al.

(2013), who used top-weighted functions to determine the particle release

in the source column. While this method reduces the error of how mass is distributed vertically, it increases the error in the initial concentrations. Top- weighted column sources will result in higher concentrations at higher altitude, which is the opposite of what one would expect from a convective eruption column. This problem was addressed in Paper II, where the source volume was designed to maintain a fixed emission rate per unit-volume, while still re- leasing most of the ash at high altitudes. A Poisson distribution was used to determine the emission rate at different altitudes. The mass release rate was then used to determine the horizontal extent of the source volume at each level according to:

d = M ˙ (z i )

U(z i ) ˜CΔz (3.1)

where ˙ M(z i ) is the mass release rate in a given source segment (z i ), U (z i ) is the average wind speed within a segment, Δz is the segment thickness, and C is a constant proportional to the concentration within each segment. This ˜ method reduces the overestimation of concentration in the upper portions of the plume. The method could be further improved by including atmospheric density profiles.

3.2 Regional exposure to passive degassing

Exposure to gases emitted from Mt. Nyiragongo (eastern D.R. Congo) was studied in Paper III. Mt. Nyiragongo is an active volcano which has been emitting large amounts of SO ₂ since its most recent eruption in 2002. Ground based monitoring of SO 2 in the plume has been in place since 2007 as part of the global Network for Observation of Volcanic and Atmospheric Change (NOVAC) since 2007 (Galle et al., 2010). The measurements are made by four scanning Differential Optical Absoroption Spectroscopy (DOAS) sta- tions. The stations are located 10–15 km south-west of the summit of Mt.

Nyiragongo. By combining DOAS measurements with wind profiles above the summit, it is possible to calculate the SO ₂ -flux and altitude of the plume;

this was done by Arellano et al. (2016), who used meteorological data pro- duced by WRF-simulations from Paper III.

The main focus of Paper III was the regional distribution SO 2 . The WRF- model and FLEXPART-WRF were used to simulate the dispersion of emis- sions from Mt. Nyiragongo. A map over the studied area is shown in Fig- ure 3.1B; note the mountain-range extending from north to south (the Alber- tine rift) and lake Kivu which have an important influence on the dynamics of the region. Also note the close proximity of Nyiragongo to the city of Goma.

Meteorological data was produced by dynamically downscaling ERA-Inter-

im, using the WRF-model in four steps to a final resolution of 2 ×2 km. The

data was then used in dispersion simulations using FLEXPART-WRF. Emis-

sion rates and altitude were determined based on the DOAS-observations;

however, since the observations can only be made during daylight, there were regular nightly gaps between observations. Due to various reasons there were also occasional longer gaps in the time series. Since the dispersion model re- quired more than 12 hours to spin-up properly, a technique for filling the gaps was required in order to study the long-term exposure.

A technique was developed where fluxes and emission heights were ran- domly sampled from the rest of the time-series whenever a gap longer than 1 hour was encountered; each sample was used to cover a period of ∼30 min.

In order to estimate the error introduced when using this technique, the process was repeated 30 times, creating 30 times series of emission rates and heights.

Each time series was used in a separate dispersion simulation, which — when combined — formed an ensemble of 30 members.

Dispersion simulations were made covering one year and seasonal changes were identified. Near-surface concentrations were calculated for a number of populated areas downwind of the volcano.

3.3 Inverse modelling of SO ₂ emissions

Determining the emission rate and plume rise is crucial for accurate mod- elling of atmospheric dispersion from volcanic eruptions. The methods de- scribed in Chapter 2, i.e. (2.10) and (2.11), attempt to counter this but involve large uncertainties. An alternative method was developed based on (2.15) as part of Paper IV, where inverse modelling was used to determine the plume height and eruption rate of the Holuhraun eruption in 2014–2015. With the then available high-resolution (down to 0.1 degrees) global analysis data from ECMWF, it was also possible to compare results from FLEXPART-WRF and FLEXPART.

Dispersion simulations were configured with a number of potential emis- sion sources as illustrated in Figure 3.2. The potential emission sources were used to simulate a-priori plumes without any detailed information of plume rise or emission rate. The a-priori results were used in an inversion together with traverse measurements of the plume. The measurements were conducted using a mobile-DOAS instrument and were presented by G´ıslason et al. (2015).

Each traverse consisted of a series of around 1,000 measurements of the col- umn densities of SO ₂ , made by a mini-DOAS instrument mounted on a car.

The car traversed the plume along the ringroad on Iceland at a distance of 80–

240 km from the eruption site, depending on the current wind direction. The dispersion models (FLEXPART and FLEXPART-WRF) were used to calculate modelled vertical columns along the same transects of the plume as covered by the car traverses. Each individual measurment in each traverse was matched with a corresponding modelled value.

Two different inversion algorithms were tested: a non-negative least-squares

(NNLS) regression and modified version of the Eckhardt scheme. The main

2 4 6 8 10 12 6 (km)

5 4 3 2 1

(33)

(14) column

sources

surface sources

Figure 3.2. Illustration of potential emissions sources used for inverse modelling of emissions from the Holuhraun eruption (Paper IV). The setup consists of 33 stacked column sources located above the main fissure vent and 14 surface sources covering the lava field.

modification made was to the smoothing term, i.e. J ₃ in (2.15), which was adjusted to only apply to the column sources. Inversions were compared for time-averaged and time-interpolated model results. Finally, the posteri- ori emissions were used in extended simulations and compared with satellite instruments from the Ozone Monitoring Instrument (OMI) on NASA’s Earth Observing System, the AURA satellite (Levelt et al., 2006).

3.4 Model enhancements

In Paper III, we discovered problems with the dry-deposition of gases in FLEXPART-WRF. First, the land-use data supplied with the model was at a horizontal resolution of 0.3 ^◦ , while the output grid used a resolution of 2 km.

High-resolution land-use data was made available in FLEXPART-WRF as part of Paper III by adding support for reading this data from WRF. This so- lution did not completely solve the problem until it was discovered that an ill-formated table prevented roughness lengths, z ₀ , less than 1 m from being properly initialized. This in turn resulted in the aerodynamic resistance (r _A ) being undefined for all land-use types except water (where Charnock’s rela- tionship is used). Furthermore, a missing update of the table (which should have been included in version 3.0) caused the Charnock relation to be applied to a non-water land-use type.

Neither of the land-use problems were obvious in Paper II since the study

focused on particulate matter which can still be deposited through the settling

velocity according to (2.5). However, in Paper III the focus was on gases

which only have one deposition pathway through the surface layer according

to (2.4). After correcting z 0 , a comparison was made between one simula-

tion using the old land-use data and one using the new. Using the set-up

from Paper III, two simulations were made, each covering a one-year pe-

riod. The accumulated dry deposition from each of the land-use data sets were compared (Figure 3.3A). The dry deposition using the new land-use deviated locally from the original by ±50 %. Differences were also seen in the aver- age concentration of SO ₂ below 500 m altitude (Figure 3.3B), mainly due to different land types directly west of the emission source.

0.75 1.25

0.95 1.05

0.85 1.15 1.5

1.3

1.1 0.9 0.7 0.5

(B) (A)

g/m³ g/m³

g/m² g/m²

Figure 3.3. Difference between FLEXPART-WRF simulations using different land- use data (figure from Paper III). Both figures show average values from a simulation with land-use data from WRF divided by results from a simulation using the old built- in data. (A) shows differences in accumulated dry-deposition of SO 2 over one year.

(B) shows difference in average SO 2 -concentration 0-500 m a.g.l. during September–

November 2010.

4. Results

4.1 Ash hazard climatology

The probability of exposure to ash as a result of future volcanic eruptions on Iceland were studied for different eruption scenarios in Paper II. Results from the Askja scenario (see Table 3.1 for details) are shown in Figure 4.1;

the development over time is shown from top to bottom (four 24-hour pe- riods). Figure 4.1A–D shows the probability of exposure for eruptions oc- curring any time of the year, while E–H and I–L represent wintertime and summertime eruption, respectively. The area of highest probability of expo- sure is seen to propagate eastward over time. The hazard varies with seasons, with the strongest difference seen between summer (June–August) and winter (December–January). Wintertime eruptions have a higher probability of af- fecting most of the Scandinavian peninsula and areas around the Baltic Sea.

Since the polar front is weaker during summer, the westerly winds are weaker and less persistent, resulting in a higher probability of ash being transported in other directions. In winter, however, the strong polar front, results in strong westerly winds most of the time, transporting ash eastward in most of the cases. This transport is both more frequent as well as more efficient, resulting in high concentrations of ash at greater distances from the source.

A comparison was made between different eruption scenarios, an example of this is shown in Figure 4.2A-E, using a threshold value of 0.2 mg/m ³ ; this threshold corresponds to the no-flight condition in use prior to the eruption of Eyafjallaj¨okull in 2010 (Webster et al., 2012). Note that the Eyafjallaj¨okull eruption case (covering the most intense phases of the eruption) showed low probabilities (below 10 %) of exceeding the threshold values for most land areas beyond Iceland. The weakest scenario (based on Mt. Spurr) shows no exceedances for this period, partly due to the short duration of the eruption.

In general, ash at higher altitude is more consistently transported eastwards compared to lower levels. However, the longest transport is seen in the mid- level ( ∼6–11 km) where the jet stream is expected.

A comparison is also made with earlier results by Leadbetter and Hort

(2011) (Figure 4.2F), who used a similar approach. The scenario set up by

Leadbetter and Hort (2011) should, in terms of emission strength, correspond

to the weakest scenario in Paper II (i.e. Figure 4.2A). However, their results

are more similar to the Askja-1875 scenario (Figure 4.2C), suggesting that the

ash hazard is lower than previously estimated.



  



















































































    

    

    

    

    

    

    

    

    

    

    

    

















































Figure 4.1. Probability of ash concentration exceeding the no-flight threshold (2.0 mg/m ³ ) for an eruption similar to Askja-1875 (from Paper II). This scenario is de- signed to be roughly 10 times stronger than the Eyjafjallaj¨okull eruption in 2010. The probability is calculated for three different periods: yearly (A–D), winter (E–H) and summer (I–L). Different time periods are given, specifically (top to bottom) 0–24, 24–

48, 48–72 and 72–96 hours relative the onset of the eruption. All cases are for flight

levels FL200–FL350 (∼ 6–11 km a.s.l.).