Simulation of the dispersion of the Eyjafjallajökull plume over Europe with COSMO-ART in the operational mode

Introduction Conclusions References


Introduction
After resting for 187 yr, the Iceland volcano Eyjafjallajökull woke up again on 20 March 2010.Starting on 14 April, massive emissions of volcanic ash occurred and finally led to a shutdown of civil aviation over entire Europe.
The emissions went on with variable strength until 23 May 2010.The volcanic eruption represents a unique field experiment for investigating atmospheric processes, such as transport, radiation, and cloud formation, on a large variety of scales by combining observations and numerical models.In addition to the academic interest, the huge economic costs of the shutdown of civil aviation raised the need for accurate forecasts of the temporal and spatial distribution of the ash plume.Quantitative numerical forecast of the spatial distribution of volcanic ash particles requires the knowledge of the source strength of the particles and gaseous precursors, the vertical distribution of the effective source heights, the size distributions at the source, and a forecast of these parameters.Introduction

Conclusions References
Tables Figures

Back Close
Full Almost parallel to the on-going eruption, numerical simulations were started with different numerical model systems by different research institutions and weather services.With the exception of the forecast done by the Volcanic Ash Advisory Centre, none of these model forecasts was quantitative.Most of these model systems were coupled offline, i.e. the aerosol transport models were driven by pre-calculated outputs of meteorological forecast models (e.g.Folch et al., 2011;Matthias et al., 2012;O'Dowd et al., 2012b).During specific episodes, aircraft measurements were carried out in addition to the remote sensing observation data taken by LIDARS, sun photometers, ceilometers, and satellites (Flentje et al., 2010;Ansmann et al., 2010;Emeis et al., 2011;Gasteiger et al., 2011;Schumann et al., 2011;Langmann et al., 2012;O'Dowd et al., 2012a).
After the shutdown caused by the volcano, many efforts have been made to transform the qualitative model results and the measurements into reliable quantitative results.Moreover, some of those model results in combination with the observations were even used to quantitatively calculate the temporal development of the source strength and the vertical emission profile (Stohl et al., 2011).
Already during the volcano eruption, we transferred the online coupled model system COSMO-ART into an operational version.After a few days of running the model in the hindcast mode, forecasts were launched every six hours for a period of at least 48 h with the same integration procedure as used for the weather forecast at Deutscher Wetterdienst (DWD).With this procedure, we produced a simple time-lagged ensemble that did not only give the spatial, temporal, and size distributions, but additionally allowed for a calculation of probabilities of a certain threshold concentration being exceeded.As quantitative measurements were published recently, we used those measurements to recalibrate our model results for subsequent calculation of absolute values of the volcanic ash number density and the mass concentration of individual size classes.
Additionally, we performed several sensitivity runs.
In the following sections we will describe the model system, the procedure of time integration, and some qualitative comparisons of model results with observations.Then, Introduction

Conclusions References
Tables Figures

Back Close
Full we will present our calibration procedure, comparisons of absolute values, and probability plots for areas exceeding the current aviation threshold values of volcanic ash.

The model system, procedure of integration, and input data
We transferred the comprehensive fully online-coupled model system COSMO-ART (Vogel et al., 2009) used for research purposes into the operational forecast mode within the experimental forecast system at DWD.Here, online-coupled means that the Eulerian prognostic equations for e.g.volcanic ash are solved and integrated time step by time step in the same way as it is done for the moisture scalar variables.COSMO-ART is the extension of the operational weather forecast model of DWD (Baldauf et al., 2011).Operational mode means that for a forecast period of up to 78 h forecasts were launched every six hours with the same model setup as applied for the numerical weather forecast.The forecast runs started from the operational analyses for which current observations of standard meteorological variables were assimilated via the nudging technique.Six individual size distributions with diameters between 1 and 30 µm were simulated.Deposition, sedimentation, and below cloud scavenging were taken into account.The vertical distribution of the source strength and the source strength itself are crucial parameters.Source heights were taken as published by the volcanic ash advisory centre London (VAAC, http://www.metoffice.gov.uk/aviation/vaac/),UK that is responsible for the official forecast of ash coming from volcanoes in Iceland according to international agreements.
Running the dispersion model in a forecast mode also allows for studying the predictability of the plume.As we obtained a time-lagged ensemble of up to 12 members for each specific date, it is possible to present not only the current forecast, but also a measure of the reliability of the forecast.Introduction

Conclusions References
Tables Figures

Back Close
Full 2.1 The model system COSMO-ART (Vogel et al., 2009;Stanelle et al., 2010;Bangert et al., 2011Bangert et al., , 2012;;Knote et al., 2011;Athanasopoulou et al., 2013;Knote and Brunner, 2013;Lundgren et al., 2013) is a recently developed comprehensive online-coupled model system.The model system was and still is developed to quantify the interaction of gaseous and particulate air constituents with the state of the atmosphere on the regional scale.As the model system so far has been used to carry out simulations during time periods in the order of weeks, it is currently transferred into a model suitable for climate simulations.Within COSMO-ART, the treatment of natural (e.g.mineral dust and sea salt) and anthropogenic (soot, secondary organic and inorganic) aerosol particles is usually based on the so-called modal approach.In case of the volcanic ash, we used a sectional approach.In practice that means that we described the size distribution of the number density of the volcanic ash particles by six individual size bins that were centered at 1, 3, 5, 10, 15, and 30 µm.Coagulation between the size bins was not taken into account.Chemical reactions that might change the chemical state of the ash particles were neglected.The atmospheric concentrations of particles emitted by the volcano are modified by advection, turbulent diffusion, sedimentation, deposition, and washout by below cloud scavenging.

Emission data
Crucial input parameters to calculate the volcanic ash plume are the emission data.In our case, we need the plume height at the source, the vertical profile of the emission strength, and the size distribution of the ash particles at the source.The only parameter that was determined with certain accuracy was the plume height at the eruption site.We used the observations that were published on the VAAC London website.The information on the plume height was updated every six hours.Figure 1 shows the temporal development of the plume height used for our simulations.For the source strength, we made the assumption of a constant emission rate for the number density on all vertical Introduction

Conclusions References
Tables Figures

Back Close
Full model levels.The absolute value of the source strength was scaled linearly depending on the plume height.The actual plume height was kept constant until a new observation was available.

Procedure of integration
The operational production of 78 h forecasts at 00 and 12:00 UTC and 48 h forecasts at 06:00 and 18:00 UTC started on 23 April 2010, i.e. nine days after the first significant eruption.For a good representation of the volcanic ash field as a starting point of our forecasts, the nine days' history period was simulated in the hindcast mode.This means that we simulated the whole period in a single run and to obtain the best possible meteorological simulation, we used own operational analyses of COSMO-EU as boundary data.Since the synoptic situation was characterised by a north-westerly flow, we were able to use the standard domain of COSMO-EU with Iceland located close to the north-western boundary.
In the forecast mode the model is driven by boundary data of the operational GME (global model of DWD).The individual model runs were started from the operational analysis data.The operational data assimilation cycle consists of 3-hourly runs of the COSMO model, where the basic meteorological fields are nudged to the observations.These assimilation runs are complemented by own surface, sea surface, temperature, and snow analyses as well as by an additional soil moisture analysis which aims at improving the 2 m-temperature forecasts.The forecast runs are started from own main run analyses with shorter cut-off times for the observations.
Since the prognostic volcanic ash fields are not part of the operational data assimilation cycle, the 6-hourly forecasts of volcanic ash were handed over from one simulation run to the next.For studying the pure forecast mode, we repeated the model runs four times a day starting with 14 April 2010, i.e. not starting from a hindcast after nine days.Introduction

Conclusions References
Tables Figures

Back Close
Full

Model results and sensitivity studies
In the following section we will present the uncalibrated hindcast model results and a qualitative comparison with observations.During the first days of the eruption, volcanic ash was injected into the atmosphere up to a height of 11 km.Consequently, it was transported rapidly at higher levels towards Europe.A comparison of the simulated number concentration of the 1 µm bin at 7500 m above surface with the satellite pictures shows that the model captures the horizontal distribution of the ash plume quite well (Fig. 2).The simulated number concentrations are given in arbitrary units.The volcanic ash that was located above a narrow band of clouds is reproduced.Nevertheless, the exact location of the plume is hard to judge from these and other satellite data, as the volcanic ash is obscured by clouds.Profiles of LIDAR backscatter coefficients were measured during the eruption at several sites.In the following section, we will compare our model results with the LIDAR measurements taken close to Munich (Gasteiger et al., 2011).
We performed a set of sensitivity runs to assess the effects of different uncertainties or physical processes on the model results.This includes runs to quantify the impact of an increased vertical resolution, the contribution of washout by below cloud scavenging, the effect of sub grid scale convection, and the effect of different vertical emission profiles.
Very often, the observations of volcanic ash plumes show thin vertical layers of a few hundred metres.Atmospheric Eulerian models, such as COSMO-ART, usually use a telescoping grid in the vertical direction.As the vertical grid size increases with height, it is difficult to resolve these thin layers.In order to study the impact of the vertical resolution on the model results, we increased the number of vertical layers from 40 used operationally to 80.This already leads to changes in the meteorological variables and, consequently, in the concentration distributions.Additionally, artificial vertical mixing caused by the vertical grid size is reduced when the vertical resolution is increased.
Figure 3 shows vertical profiles of the simulated particle number concentration for the Introduction

Conclusions References
Tables Figures

Back Close
Full individual model runs.Figure 3a shows the results for the reference run (40 vertical layers).The model results with 80 vertical layers (Fig. 3b) show a decrease of the vertical thickness of the ash plume that corresponds to the decrease of the thickness of the vertical model layer at the corresponding height.To give one example, the thickness of the plume decreases by almost 800 m at 03:00 UTC on 17 April 2010.By increasing the vertical resolution, it is therefore possible to reproduce the vertical thin layers corresponding to the observation.Both simulations calculate particle concentrations greater than zero not only at higher elevations, but also in the boundary layer.The model results can be compared to the backscatter coefficient measured at the location of the Munich LIDAR (Gasteiger et al., 2011).These observations are given in Fig. 4.
The simulated ash layer above the boundary layer starts to decrease around midnight on 17 April 2010 and finally reaches 2000 m at the end of the given time interval.Based on the results of the model run with 80 vertical levels (Fig. 3b), a vertical subsidence of 8.5 cm s −1 is calculated between 01:30 UTC and 08:00 UTC and of 2.35 cm s −1 between 08:00 UTC and 21:00 UTC of 17 April 2010.Calculating the same numbers from the LIDAR profile of Gasteiger et al. (2011) gives 5.5 cm s −1 between 01:30 UTC and 08:00 UTC and 1.7 cm s −1 between 08:00 UTC and 17:00 UTC (Fig. 4).
In order to study the impact of modified emission profiles, another model run was performed.In this run we changed the emission profile from constant emissions with height to a more mushroom-like vertical emission pattern following Oberhuber et al. (1998).Figure 3c shows the results for this simulation.The largest differences occur in the boundary layer where the concentrations are much lower.
Figure 5 shows the time series of the uncalibrated number density for the 1 µm bin of the individual sensitivity runs at the surface at a grid point that corresponds to the location of the observatory Hohenpeißenberg, Germany (47 • 48 N, 11 • 01 E).The results show that for that period the parameterization of deep convection does not change the result with the exception of small variations on 21 April 2010.Larger differences in comparison to the reference run occur when washout is neglected especially between Introduction

Conclusions References
Tables Figures

Back Close
Full From the practical point of view, there is great interest in knowing not only the areas that are affected by the ash cloud of a volcanic eruption, but also the absolute values of the concentrations.During the Eyjafjallajökull eruption event, such quantitative forecasts were hardly possible.The most important reasons were the unknown source strength and the shape of the source profile as well as lacking observations of absolute concentrations.However, the measurements during the event allow for a recalibration of the model results.The procedure applied to do this recalibration is described below.
Measurements with the DLR Falcon aircraft were carried out during several time intervals (Schumann et al., 2011).On 2 May 2010, the aircraft flew over the North Atlantic into the top part of the fresh ash cloud.Among other quantities, size distributions of the particles were measured (Fig. 6, left).We used the measured size distributions of the number density to calculate the recalibration factors of our simulated size distributions.
The simulated number densities of the individual size bins were multiplied by these factors.On 17 April 2010, the Falcon measured vertical profiles of the number density close to Leipzig, Germany (Fig. 6, right).We used the measured number densities for particles with a diameter greater than 1.5 µm and summed up the normalised number concentrations of our size bins 3, 5, 10, 15, and 30 µm for that location and time.For heights above 3000 m we determined a factor of roughly 100 for the final calibration of the simulated number densities (Fig. 6, right).Below 3000 m there is still a large difference between measured and simulated values.These differences can be explained by non-volcanic aerosol which is not covered by the simulations.This is confirmed by CO measurements that are also presented in Schumann et al. (2011).By this procedure, Introduction

Conclusions References
Tables Figures

Back Close
Full we finally transferred our results into quantitative ones.In order to validate the recalibrated number densities, we used observations that were carried out by DWD at the observatory Hohenpeißenberg which is located at roughly 400 km distance to Leipzig. Figure 7a shows the temporal development of the measured number density of particles with a diameter of 3.5 µm and the simulated (calibrated) number densities for the size bin 3 µm.The results of the individual forecast runs are shown.With respect to the level of agreement, three time periods can be distinguished (labelled A, B, and C, Figure 7a).During time period A, the model results differ substantially from the observations.We found that the quality of the forecast increases for shorter forecast times.The most excessive overestimation seen in period A is related to the final hours of a 78 h forecast run.This shows the huge impact meteorology has on ash dispersion.During the three days of period B, the agreement of the individual forecast results with the observations is much better and the spread is reduced.Again, it is obvious that especially the forecast results with the greatest time lag are most different from the observations.During the beginning of period C and especially during the first three days, the model results differ totally from the observations.Although some temporal patterns might be visible, no quantitative agreement is reached anymore.We tracked this deficiency back and found that the extremely high simulated number concentrations have their origin at the volcano three to four days earlier.
Starting with that point of time, the height of the plume decreased drastically from 9 km to 5 km.Due to our initially assumed linear relation between source height and source strength, the emitted mass decreased linearly, by about 45 %.Based on the findings of Mastin et al. (2009), we modified our source function as follows: h is the top of the plume height at the location of the volcano in m.This formulation guarantees that the emissions are unchanged when h is 8500 m, otherwise they are modified according to Mastin et al. (2009).With this new parameterisation of the 13448 Introduction

Conclusions References
Tables Figures

Back Close
Full emissions, we repeated our simulations.Figure 7b shows the results for the grid point Hohenpeißenberg for comparison with Fig. 7a.Within the periods A and B, the results differ only slightly.Starting on 21 April 2010, the results with the modified source strength start to differ significantly from the old ones.Within period B, an improvement of the model results compared to the observations is not always found.However, the model runs with the shortest time lags give better results, thus showing the influence of meteorology again.Although the difference between model results and observation decreases during on 21 April, it is still unacceptably large.On 22 April and following days, the results are in almost perfect agreement with the observations.The results for 21 April again questioned our source function.Stohl et al. (2011) derived emission data for the 2010 Eyjafjallajökull eruption using inverse modelling techniques.Figure 2 of the original paper compares the a priori estimated emission strength based on a plume model of Mastin (2007) with the a posteriori emission from the inversion algorithm.From 16 April to 19 April 2010, huge differences up to a factor of ten between both emissions are found.Based on this finding, we reduced our emission during the time period with the biggest differences, i.e. from 17 April 00:00 UTC until 18 April 12:00 UTC by a factor of 5.0 and repeated our simulations.Figure 7c shows the results of this simulation.The large deviation of the model results from the observations on 21 April disappeared totally.During time period C, the model results are somewhat lower than the observation for this model run.Overall, the agreement with the modified source function is quite good.

Simulated mass densities and their probabilities
Civil aviation and official air traffic control authorities depend on the ash forecasts for flight planning and decision making.Numerical forecasts of volcanic ash advisory centres usually deliver spatial distributions of ash concentrations at different flight levels.
Those ash forecasts are associated with large uncertainties.We explained the reasons of those uncertainties and presented examples of differences due to those uncertainties.Based on our findings, we propose to give probabilities of exceeding a certain 13449 Introduction

Conclusions References
Tables Figures

Back Close
Full threshold value instead of absolute concentrations in case of volcanic eruptions.A method to determine such probabilities is described below.The four daily model runs with 78 h forecasts at 00:00 and 12:00 UTC and 48 h forecasts at 06:00 and 18:00 UTC result in up to twelve realisations for a specific target date.These different realisations form a time-lagged ensemble that allows us to calculate certain probability measures.
In Fig. 8 the probability of the maximum mass concentration of volcanic ash between flight level 200 and 350 (approx.6100-10 675 m) exceeding a threshold of 2000 µg m −3 is given for 16 April 2010 at 12:00 UTC.For this date, the ensemble consists of 10 members.To improve the statistical basis for the calculation of the probability, an additional upscaling using a 6 × 6 neighbourhood of each grid point was performed, resulting in 360 values for all points.The shown quantity, i.e. the maximum mass concentration in certain height ranges bounded by specific flight levels, is a standard product of the VACCs.Using the given ensemble approach we are able to provide chosen statistical measures in addition to the simple deterministic ones.

Conclusions
The research model system COSMO-ART was used in an operational forecast mode at DWD already during the on-going eruption.By these operational forecasts, we produced a time-lagged ensemble.For the first simulations, a linear dependence of the source strength on the plume height at the source was given.As long as no quantitative measurements were available, only qualitative concentration forecasts of volcanic ash could be produced.Sensitivity studies to quantify the impact of source profile, washout, and vertical resolution were carried out to quantify the uncertainties due to these input parameters and processes.Thin vertical layers of ash as reported by observations can only be represented when the vertical resolution is compared to the value currently used in the operational weather forecast of DWD.A comparison of the qualitative results with backscatter ratios measured by a LIDAR system close to Munich shows a descending ash plume in good agreement.The rate of descent was used to calculate Introduction

Conclusions References
Tables Figures

Back Close
Full vertical velocities from the model results and from the observations, which differed by a factor in the order of 1.5 only.
Based on quantitative measurements with the DLR Falcon that became available after the volcanic eruption, we calibrated our model results using observations of the particle size distribution close to the volcano and vertical profiles at Leipzig.These calibrated model results were then compared with surface observations at Hohenpeißenberg.The comparison reveals a reasonable agreement as long as the plume height of the volcano did not change.Big differences occurred after the plume rise decreased by a factor of 2, indicating that the assumption of a linear relation between the source strength and the plume height is far from reality.We repeated the simulations taking into account the nonlinear relation between source strength and plume height by Mastin et al. (2009).Although this new relation gave much better results, there was still one day left, on which our ash concentration forecast exceeded the observed values by a factor of 15.Using a source strength correction factor based on the inverse modelling work by Stohl et al. (2011), we again corrected the source function.After this correction, our model results were in nearly perfect agreement with the observations.Based on our time-lagged ensemble, we calculated probabilities of a certain threshold value being exceeded.This information is of high value for aviation in addition to the already available absolute concentration distributions published by the Volcanic Ash Advisory Centres.
In cases of future volcanic eruptions, several requirements should be fulfilled to produce optimal volcanic ash forecasts for decisions on the closure of the airspace.In situ measurements with aircrafts close to the volcano are necessary to obtain information about the size distribution.Together with concentration measurements further downstream, those observations are necessary for the calibration of the model results.
Both require operational measurements with well-equipped aircraft which should be ready for operation within a short term.Finally, remote sensing data from LIDARS and ceilometers should be used for data assimilation.The latter requires model-specific developments.Introduction

Conclusions References
Tables Figures
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 19 and 20 April 2010.Increasing the number of vertical layers leads to changes in the appearance of the concentration maxima on 17 and 21 April 2010.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 1 :
Figure 1: Temporal development of the height of the top of the plume in April 2010 at Eyjafjallajökull as used in the simulation.Data was obtained from the reports of VAAC London.

Fig. 1 .Figure 2 :
Fig. 1.Temporal development of the height of the top of the plume in April 2010 at Eyjafjallajökull as used in the simulation.Data was obtained from the reports of VAAC London.

Fig. 3 .
Fig. 3. Vertical cross-section of the uncalibrated ash concentration (arbitrary units) at Munich for different model configurations.(a): reference run with 40 vertical layers, (b): model run with 80 vertical layers, and (c): results with the modified vertical emission profile indicated in the inlay of (c).

Figure 4 :Fig. 4 .Figure 4 :
Figure 4: Logarithm of range-corrected signal of MULIS at = 1064 nm over Maisach from 16 April 2010 17:00 UTC to 17 April 2010 17:00 UTC and from 0 to 10 km above ground; white areas denote periods without measurements (taken from Gasteiger et al., 2011).The two arrows indicate the descent of the plume during two time intervals.

Figure 5 :
Figure 5: Temporal development of the uncalibrated simulated ash concentrations (arbitrary units) at Hohenpeissenberg for the model configurations indicated in the legend.

Fig. 5 .
Fig. 5. Temporal development of the uncalibrated simulated ash concentrations (arbitrary units) at Hohenpeissenberg for the model configurations indicated in the legend.

Figure 6 :
Figure 6: Left: measured size distributions close to the volcano on 2 May, 2010.Right: measured (red) and simulated (blue) number concentrations for particles greater than 2 µm at Leipzig on 17 April, 2010.Measured data taken from Schumann et al., 2011.

Fig. 7 .
Fig. 7. Temporal development of the simulated (thin) and observed (thick) number density of the 3 µm particles.(a) linear relation between source height and source strength, (b) relation between source height and source strength following Mastin et al. (2009), and (c) correction of the source strength according toStohl et al. (2011).The thin lines give the results of the forecast runs that are started with a time lag of 6 h.