Digital Commons @ Michigan Tech Digital Commons @ Michigan Tech Determination of time-and height-resolved volcanic ash emissions Determination of time-and height-resolved volcanic ash emissions and their use for quantitative ash dispersion modeling: The 2010 and their use for quantitative ash dispersion modeling: The 2010 Eyjafjallajökull eruption

. The April–May, 2010 volcanic eruptions of Eyjafjallaj ¨ okull, Iceland caused signiﬁcant economic and social disruption in Europe whilst state of the art measurements and ash dispersion forecasts were heavily criticized by the aviation industry. Here we demonstrate for the ﬁrst time that large improvements can be made in quantitative predictions of the fate of volcanic ash emissions, by using an inversion scheme that couples a priori source information and the output of a Lagrangian dispersion model with satellite data to estimate the volcanic ash source strength as a function of altitude and time. From the inversion, we obtain a total ﬁne ash emission of the eruption of 8.3 ± 4.2 Tg for particles in the size range of 2.8–28 µm diameter. We evaluate the results of our model results with a posteriori ash emissions using independent ground-based, airborne and space-borne measurements both in case studies and statistically. Subsequently, we estimate the area over Europe affected by volcanic ash above certain concentration thresholds relevant for the aviation industry. We ﬁnd that during three episodes in April and May, volcanic ash concentrations at some altitude in the atmosphere exceeded the limits for the “Normal” ﬂying zone in up to 14 % (6–16 %), 2 % (1–3 %) and 7 % (4–11 %), respec-Correspondence tively, of the European area. For a limit of 2 mg m − 3 only two episodes with fractions of 1.5 % (0.2–2.8 %) and 0.9 % (0.1–1.6 %) occurred, while the current “No-Fly” zone criterion of 4 mg m − 3 was rarely exceeded. Our results have important ramiﬁcations for determining air space closures and for real-time quantitative estimations of ash concentrations. Furthermore, the general nature of our method yields better constraints on the distribution and fate of volcanic ash in the Earth system.

source parameters are required to obtain reliable results. Volcanoes exhibit a broad range of eruptive styles and variability (Woods, 1995), thus making theoretical attempts at predicting source parameters challenging. To date, only indirect methods are available to estimate ash emission rates (Mastin et al., 2009). For example, if eruption column heights are known, e.g., from weather radar measurements (Lacasse et al., 2004), empirical relationships may be used to estimate the mass flux of tephra (Sparks et al., 1997;Mastin et al., 2009). However, these relationships are loose , the vertical profile of ash emissions remains poorly quantified, and only a small but highly variable fraction of the tephra is fine grained enough (<30 µm) for long-range atmospheric transport.
In this paper, we present an objective method to determine the volcanic ash emission rate, which is based on inverse modeling constrained by satellite measurement data. The eruption of the Eyjafjallajökull volcano (19.61 • W, 63.63 • N, 1666 m a.s.l.) in Iceland in the year 2010 represents an ideal test case for our method, as this eruption released large amounts of ash while a wealth of measurement data have been collected during the event, which are available for model evaluation. After 18 years of intermittent seismic unrest (Dahm and Brandsdóttir, 1997), an effusive eruption of basalt on the eastern flank of Eyjafjallajökull occurred from 20 March to 12 April 2010, followed by an explosive eruption under the Eyjafjallajökull glacier on 14 April 2010. The interaction of magma and ice augmented explosive activity and generated large proportions of fine ash that were emitted into the atmosphere (Gudmundsson et al., 2010). Intense tephra discharge continued for several days and the prevailing meteorological conditions resulted in ash transport directly towards Europe, where air traffic was grounded for several days. The eruption strength increased again in May, leading to further air space closures. The Eyjafjallajökull eruption demonstrated how susceptible aviation is to volcanic eruptions in Iceland, as already suggested beforehand (Sveinbjörnsson, 2001).

Methods
The inverse method to determine volcanic ash emissions merges a priori information on the ash emissions, satellite observation data and simulations with a dispersion model to derive improved a posteriori ash emissions. In this section, we describe the various data sets and methods used for the inverse modeling.

Satellite data
Observational constraints on volcanic ash emissions were provided through infrared satellite retrievals of total column airborne ash loadings (Prata, 1989;Prata and Grant, 2001;Clarisse et al., 2010a) using data from two different satel-lite instruments, the geosynchronous Meteosat Second Generation (MSG) Spin-stabilised Enhanced Visible and Infrared Imager (SEVIRI) and the polar-orbiting MetOp Infrared Atmospheric Sounding Interferometer (IASI). These two instruments combine high temporal coverage (SEVIRI has a sampling time of 15 min but hourly averages were used here), with an enhanced sensitivity to ash (IASI has over 1000 spectral channels which can be used for ash detection).

SEVIRI
SEVIRI observes the earth disk over a total field of view of 70 • in 12 channels from the visible to the infrared. The SE-VIRI temporal sampling time is 15 min with a spatial resolution of ca. 10 km 2 at the below-satellite location which increases to ca. 100 km 2 near the edges of the scan. Only infrared channels were used in the analyses for a sub-region of the SEVIRI disk covering the geographical region 30 • W to 30 • E and 40 • N to 70 • N. The retrieval relies on measurements made at 10.9 µm and 12.0 µm, that are corrected for the effects of absorption and emission by atmospheric water vapour (Yu et al., 2001) and then inverted to determine optical depths and effective particle radius using established methods (Prata, 1989;Prata and Grant, 2001;Wen and Rose, 1994). A Mie scattering program (Evans, 1988) and a discrete ordinates model (DOM) (Stamnes and Swanson, 1981) were used to estimate the top-of-atmosphere (TOA) brightness temperatures at 10.9 and 12.0 µm, assuming a planeparallel cloud of andesitic spherical particles with uniform cloud temperature overlying a background of uniform surface temperature. Refractive indices of andesite (Pollack et al., 1973) as a function of wavelength were interpolated, convolved with SEVIRI response functions, and input to the Mie scattering code using a set of assumed modified-γ size distributions with different mean particle radii and standard deviations. The scattering code outputs a set of extinction coefficients, asymmetry parameters and single-scatter albedos for each particle size and wavelength. These values are then used in the DOM code to determine the TOA brightness temperatures. The cloud top and background surface temperatures are also required as input to the DOM code. Initial estimates of these are determined from the data by finding clear and opaque pixels and assigning these to the surface (T s ) and cloud top temperatures (T c ), respectively. As there is some error in this assignment, an ensemble of DOM calculations is performed for values of T c and T s that are ±10 K different for each initial temperature. Finally, each DOM calculation is performed in equal optical depth steps of 0.02 starting at 0 (clear field of view) and ending at 7.98 (opaque field of view). The result of these calculations is a table of simulated TOA brightness temperatures for the two SEVIRI channels and for each of the combinations of cloud top and surface temperature. Each entry in the table appropriate for the scene cloud top and surface temperature, is indexed by optical depth and effective particle radius.
Based on the water vapour corrected SEVIRI data, the table is interpolated in terms of brightness temperatures to find the closest entries corresponding to the retrieved optical depth and effective particle radius. This procedure uses both the 11 and 12 µm brightness temperatures. The retrieval provides the best fit effective particle radius (in steps of 0.02 µm) and infrared optical depth (in steps of 0.02). The mass loading (g m −2 ) can be calculated from, where ρ is the density of the ash, r e is the retrieved effective particle radius, τ λ is the optical depth and Q ext (λ) is the extinction efficiency, all functions of the wavelength, λ, leading to errors of 40-60 % in estimated mass loading (Wen and Rose, 1994). The retrievals have a preferential sensitivity to ash with particle diameters from 2 to 32 µm. The composition of the pixel (e.g. ash, water or ice cloud, clear) is identified prior to performing the ash retrieval. The main test requires that the water vapour corrected brightness temperature difference (BTD, which is a function of the atmospheric pathlength and so depends on the scan position) between the 10.9 and 12.0 µm SEVIRI channels must be <−0.8 K. In this case the pixel is assumed to contain volcanic ash (Prata, 1989). Other conditions applied include an opacity test where if the sensed brightness temperature at 12.0 µm is less than 230 K, the cloud is assumed to be optically thick and retrievals are not made.
For the inverse modeling, the 15 min pixel-by-pixel mass loading retrievals were binned into 0.25 • ×0.25 • grid cells and time-averaged to provide 240 longitude ×120 latitude values every hour. A parallax correction (Vicente et al., 2002) was applied to all ash-affected pixels assuming that the ash clouds were at 6 km height. This simplification results in a small error in geolocation, but is an improvement compared to using the data without a parallax correction.

IASI
IASI is a sunsynchronous polar orbiting infrared sounder (Clerbaux et al., 2009). With its high spectral resolution and low radiometric noise, it has proven very useful in monitoring a host of trace gases. Relatively little attention has been given to the sounding of aerosols with IASI and there are no current or planned operational products pertaining to aerosols. Infrared sounding of aerosols has, however, a number of distinct advantages, such as the availability of night time data and the high sensitivity to aerosol morphology. Due to its complexity, a sophisticated method for retrieving radius and mass loadings (Clarisse et al., 2010a), based on optimal estimation which iteratively fits synthetic spectra to an observed spectrum by varying radius and mass loading is not suitable for the treatment of large amounts of data. Another difficulty is that reliable results can only be obtained if the wavenumber-resolved refractive index of the observed aerosols is known. For the Eyjafjallajökull eruption, none of the few published refractive index data for the erupted ash matched the high resolution observations (for broadband sensors like SEVIRI this is not so much of an issue).
Recently, a sensitive method was presented for the detection of volcanic ash based on correlation coefficients (Clarisse et al., 2010b). We applied this method here for two months of IASI data for the region 60 • W to 50 • E and 45 • N to 90 • N. Mass loadings of sulfur dioxide in the upper troposphere and lower stratosphere were retrieved in parallel. Coincident measurements of sulfur dioxide provided good supporting evidence for the many puffs of ash detected. The ash detection algorithm does not provide mass loadings. As a measure for the total ash column, the BTD between the IASI channels at 1231.5 cm −1 and 1160 cm −1 (Clarisse et al., 2010b) was calculated (which is close to zero for clear or cloudy scenes and positive for ash scenes). Forward calculations based on basaltic ash optical properties show that this difference is approximately linearly proportional to the ash mass loading for all but the highest ash concentrations. The exact conversion factor depends mostly upon plume altitude and size distribution and was calibrated to match the SEVIRI retrievals of mass loadings from coincident measurements.

A priori emissions
We compiled mean eruption column heights from six-hourly (VAAC, 2010) and daily  reports, and three-hourly radar data (Petersen and Arason, 2010). To determine the erupted mass flux, we used a one-dimensional model for convective volcanic plumes (PLUMERIA) (Mastin, 2007) which considers actual atmospheric conditions taken from ECMWF data. The model was run iteratively for each three-hour interval to estimate the mass flux corresponding to observed plume heights. This calculated mass flux was then vertically distributed according to model predicted magma densities and plume radius, yielding a time-height gridded inventory with 328 3-hourly intervals for the period 14 April to 24 May, and 19 layers of 650 m vertical resolution. We assumed that 10 % of the erupted mass was fine ash in the size range of 2.8-28 µm to which satellite measurements are sensitive, obtaining a total fine-ash emission of 11.4 Tg for the 41-day period considered. When put into the dispersion model, this leads to simulated vertically integrated atmospheric ash loadings that are in reasonable overall agreement with available satellite data. Notice that our a priori emission estimate should be more accurate than common operational methods, which use statistical relationships between plume heights and total ash emission (Sparks et al., 1997;Mastin et al., 2009). These methods do not consider actual atmospheric conditions and do not model the vertical distribution of the ash.
For determining the emission uncertainties, the timeheight emissions were first smoothed by applying a tri-cubic weight function. Subsequently, the emission uncertainties were set to 100 % of the highest weighted emission within a window of 1300 m and 6 h half-widths, respectively. Furthermore, minimum values were prescribed in order to avoid zero uncertainties for grid cells without ash emissions. Since we cannot objectively determine the true values of the emission uncertainties, they were chosen such as to allow the inversion to substantially change the emissions, while still being guided by the a priori estimate.

Model simulations
To simulate the dispersion of volcanic ash, we used the Lagrangian particle dispersion model FLEXPART (Stohl et al., 1998(Stohl et al., , 2005. The simulations accounted for gravitational particle settling (Naeslund and Thaning, 1991) and wet and dry deposition (Stohl et al., 2005), but ignored ash aggregation processes. FLEXPART was driven with three-hourly meteorological data from two different sources, namely the European Centre for Medium-Range Weather Forecasts (ECMWF) analyses with 0.18 • ×0.18 • resolution and 91 model levels, and the National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) analyses with 0.5 • ×0.5 • resolution and 26 levels. Comparisons between simulation results using these two alternative data sets were used to quantify model uncertainty (see below).
To improve the a priori emissions by the inversion algorithm, it was necessary to run the dispersion model for each one of the 6232 (328 times×19 layers) emission grid cells. Each one of these 6232 scenario simulations explored the sensitivity of total atmospheric ash columns in our modeling domain (30 • W to 30 • E and 40 • N to 70 • N) to the ash emissions in a single emission grid box. The simulations extended over six days and carried 360,000 particles in nine different size bins (4,6,8,10,12,14,16,18, 25 µm diameter; particle density 3000 kg m −3 ). Per simulation, a unit mass of ash weighted with the assumed size distribution (see below) and multiplied by the sensitivity of the satellite retrievals was distributed over all particles. The model output consisted of total atmospheric columns of ash for comparison with satellite observations at a spatial resolution of 0.25 • ×0.25 • . No vertically resolved ash concentrations were determined, to keep the model output at a manageable size.
To obtain vertically resolved ash concentrations, we also ran single model simulations based on the gridded a priori and a posteriori emission fluxes. These simulations produced output at high vertical (250 m) resolution and for a global domain and extended over the full period of the eruption. For these simulations, the ash size distribution was extended beyond the range measured by satellite, using a total of 25 size classes from 0.25-250 µm and releasing 24 million particles in total. Results from these simulations were used for model validation and to quantify the ash concentrations over Europe.  The blue lines are a fit through the initial size distribution required to match the ground sample distribution with the model (thick blue line) as well as shifted distributions used for sensitivity analyses (thin blue lines). The red line shows the sensitivity of the SEVIRI and IASI satellite retrievals to ash particle size. The yellow and grey boxes indicate the size range used for the inverse modeling and for other modeling, respectively. Only the ash mass in the yellow size range (ca. 70 % of the total mass for the fitted reference distribution) is constrained by the satellite measurements.

Ash particle size distribution
To define the ash particle size distribution, we forced the dispersion model, using the a priori source term, to reproduce the measured deposit size distribution (see Fig. 1) at a location close to the volcano. This was achieved by adjusting the emitted size distribution. We then fitted two lognormal curves to the optimized emitted size distribution to specify our initial ash size distribution dM d logD p , where M is mass and D p particle diameter. We obtained a distribution with a primary mode at 10 µm diameter (geometric standard deviation 1.3) and a secondary mode at 180 µm (geometric standard deviation 0.35), and used further distributions with the primary mode shifted towards smaller (7 µm) and larger (13 µm) particle size modes for sensitivity analyses (Fig. 1). Notice that with the chosen initial reference model size distribution (thick blue line in Fig. 1), the size distribution of the deposited ash (black asterisks in Fig. 1) can be matched nearly exactly by the model. The initial size distributions are broadly consistent with measurements of airborne ash in the 18 May 1980 Mount St. Helens eruption (Hobbs et al., 1982). Modeled size distributions with maximum modal diameter >13 µm produced results in poorer agreement with satellite observations at long distances from the volcano due to increased gravitational settling and correspondingly reduced ash atmospheric residence time. Alternatively, modeled particle size distributions with a maximum modal diameter below 7 µm are inconsistent with downwind in situ measurements of ash particle size distributions in the Eyjafjallajökull ash cloud, which show larger modal diameters (Schumann et al., 2011).

Inversion algorithm
In previous studies (Eckhardt et al., 2008;Kristiansen et al., 2010), we developed an inversion algorithm to calculate the vertical distribution of sulfur dioxide emission rates for instantaneous volcanic eruptions, using only total column observations of sulfur dioxide. The algorithm extracts emission height information from the horizontal dispersion patterns, which depend on altitude because of the vertical shear of the horizontal wind. By matching the observed plume with many simulations initialized at different altitudes, a best-fitting vertical emission profile is obtained as a linear combination of all emission height scenarios. The algorithm is based on the theoretical work of Seibert (2000) and was used in modified form also for determining the spatial distribution of greenhouse gas emissions (Stohl et al., 2009). A full description of the algorithm was given previously (Eckhardt et al., 2008;Stohl et al., 2009) and a method to derive a posteriori uncertainties by propagating the uncertainties in the a priori emissions, the observations and the model calculations through the inversion algorithm was developed by Seibert et al. (2011). Here, we do not describe the algorithm again but explain the few modifications that were necessary for this study. Sensitivity tests and evaluation with measurement data were also reported in these previous studies but some more will be presented here. In this paper, we use this algorithm for the first time to yield volcanic ash emission rates. This did not require any changes in the inversion algorithm as differences in the transport and loss processes are accounted for in the FLEXPART simulations. Furthermore, it is the first time we derive emissions both as a function of altitude and time for a six-week long eruption. The few necessary changes to the algorithm are explained below.
To establish the sensitivity of atmospheric ash loadings to spatially and temporally resolved emissions, 6232 different emission scenarios were established. For each scenario, a unit amount of ash was emitted in one of 19 vertical layers stacked up to 12.3 km altitude and during one of 328 3-h time intervals between 14 April and 24 May 2010. For each scenario, FLEXPART was run to evaluate the atmospheric ash total column loadings. The model results for all scenarios were matched (i.e., ensuring spatio-temporal co-location) with about 2.3 million 0.25 • ×0.25 • gridded satellite observations, which were available hourly for the geosynchronous platform and twice daily for the polar orbiter. The matched data set was then fed into the inversion algorithm which optimally merges satellite observations, a priori emissions and simulated sensitivities (Eckhardt et al., 2008;Kristiansen et al., 2010). The result of the inversion is vertically and temporally resolved a posteriori emissions on the original emission grid, obtained as a linear combination of all scenario source terms, which optimize the agreement with both the a priori emissions and the satellite observations when considering the uncertainties of both data sets.
In our previous studies (Eckhardt et al., 2008;Kristiansen et al., 2010;Seibert, 2000), we have determined emissions only as a function of height. However, adding time did not require any changes in the inversion algorithm other than organizing the two-dimensional emission information in a onedimensional vector so that we were able to use our existing computer code. The only coding change necessary was to ensure that the vertical smoothness condition (Eckhardt et al., 2008) is indeed applied only in the vertical and not in time.
An important improvement over our previous work (Eckhardt et al., 2008;Kristiansen et al., 2010;Seibert, 2000) is that the model error is approximated for every individual grid cell as the difference between a posteriori model results based on GFS and ECMWF meteorological input data, respectively. These model results were obtained from simulations over the full study period using, by means of an iterative loop, the a posteriori emissions. While a larger model ensemble would be preferable to characterize the model error, this is an improvement over assuming an arbitrary constant model error. The model results were also used to identify regions where ash older than six days contributed more than 20 % of the total ash. Data from such locations were not used for the inversion because the sensitivity calculations extended only over six days and older ash would contaminate the inversion. All other satellite observations detecting ash were used for the inversion but 75 % of the observed zero concentration values were removed, applying a random data thinning scheme which was weighted towards keeping cases providing a strong emission constraint. In total, 2.3 million observation cases were used for the inversion.
While the inversion method formally propagates stochastic errors in the input data to the calculated a posteriori emissions, the overall uncertainty of the emissions is driven mainly not by stochastic errors but by the 40-60 % errors of the satellite ash retrievals, which are partly systematic. Therefore, we report all ash emission and concentration errors as 50 %. In particular, low-altitude ash clouds with loadings <0.5 g m −2 are often below the detection limits of the two satellite instruments used and, thus, emissions during episodes with less intense eruptions may be biased low. On the other hand, errors for thick high-altitude ash plumes may be smaller than 50 %. Notice the switch from a linear to a logarithmic scale above 10 t s −1 (yellow area). All heights throughout the paper are given in meters above sea level.

A posteriori emissions
From the inversion, we obtain a total fine ash (diameter 2.8-28 µm) emission of 8.3±4.2 Tg, about 73 % of the a priori estimate. Extrapolating the emissions to the size range of 0.25-250 µm using the size distribution shown in Fig. 1 yields a total ash emission of 11.9±5.9 Tg but we must keep in mind that the mass outside the 2.8-28 µm range is not actually constrained by satellite data. On average, the a posteriori emissions are shifted towards higher altitudes (Fig. 2a) and the temporal evolution of the emissions is considerably different from the a priori values (Fig. 2d). For example, high-intensity emissions between 16-18 April are reduced by more than a factor of four, whereas emission intensity on 12-13 and 15 May is substantially increased, due to satellite observations that are inconsistent with the a priori emissions. The results for the GFS-and ECMWF-based inver-sions are very similar, both with respect to the temporal and vertical distribution. Generally, the a posteriori emissions are more variable than the a priori emissions, both in time and altitude (Fig. 2b and c). The a posteriori emissions are released mainly in a few strong pulses, typically close to or even above the top of the a priori eruption column. During these pulses, little ash was emitted below 4 km above sea level.
The inversion reduced the root-mean-square (RMS) error between the gridded ash total columns in the ECMWF-datadriven simulation and the satellite data by 28 % when comparing a priori and a posteriori results. Visual inspection suggests that the a posteriori ash dispersion patterns are more consistent with those observed by the satellites and a large fraction of the remaining RMS error is due to noise in the satellite data, as shown below.

Emission uncertainty reduction
The a priori uncertainties of the volcanic ash emissions are assumed to be 100 % of the highest emission value in the vicinity of a space-time grid cell. These errors are reduced by the inversion due to the incorporation of observation information. The relative error reduction (Fig. 3) is typically largest where ash emissions (Fig. 2) -and, thus absolute values of a priori uncertainties -are highest. The reason for this is that the signal-to-noise ratio is larger in grid cells containing a lot of ash. Error reductions can be close to 100 %, meaning that the formal methodological a posteriori errors are very small. However, real errors are likely much larger, due to probable systematic satellite retrieval errors of 40-60 % and error correlation, which is not taken into account by the inversion method. Furthermore, the inversion scheme assumes errors to be uncorrelated, thus yielding a too large error reduction. Still, Fig. 3 provides valuable information as it indicates where the satellite data provide a good constraint on the emissions and where such constraints are probably weaker (e.g., during the period 6-11 May).

Sensitivity to ash particle size distribution
While dispersion model results are sensitive to the assumed ash particle size distribution, the impact on the estimated source term is remarkably small for the tested size distributions (Fig. 4). One reason for this is that the satellite retrievals are only sensitive to fine volcanic ash in a rather limited range of particle sizes (Fig. 1). Thus, only differences in the shape of the size distribution in this range can affect the results. Absolute differences in sedimentation velocities are relatively small for particle sizes smaller than about 10 µm. Furthermore, the inversion is guided most strongly by ob-servations of large volcanic ash columns in the vicinity of the volcano, which further limits the effects of different sedimentation rates for the chosen size distributions.

Sensitivity to changing the satellite data set
We performed ECMWF-based inversions also for subsets of the satellite data used, namely for either SEVIRI or IASI data alone. The a posteriori total fine ash emissions when using only SEVIRI (IASI) data were 7.9 (10.4) Tg, 7 % less (22 % more) than the 8.5 Tg obtained for the ECMWF-based inversion when using both data sets. The inversion using IASI data only is closest to the a priori emission of 11.4 Tg. The main reason for this is that the number of gridded IASI observations is about an order of magnitude smaller than the number of SEVIRI observations, thus providing a weaker constraint on the emissions, which therefore remain closer to the a priori values.
The emission changes made by the inversion are, however, relatively consistent for the two data sets separately and for the combined data set. For instance, the high a priori emissions from 16-18 April are considerably reduced in all cases, even though the reduction is smallest when using IASI data only (Figs. 2 and 5). All inversions also lead to substantial emission increases for 12-13 May and to a general shift of the ash emissions to higher altitudes.

Model evaluation
Model results were evaluated against a large number of independent observational datasets, which are presented in the following.

Analysis of observed plume top heights
Eruption column top heights were estimated from archived images of two webcams viewing the Eyjafjallajökull summit when the view was not obscured by clouds. Both webcams are made by MOBOTIX and were installed by Mila (http://www.mila.is). Webcam 1 was located at a distance of approximately 10 km and has a field of view of 8.2×6.1 km. The maximum visible altitude is about 4.8 km. Webcam 2 was positioned roughly 15 km from summit and has a 17.25×12.3 km field of view. The maximum visible altitude is about 7.8 km. When the plume was clearly visible, plume top altitudes were estimated from geometrical principles, taking into account the camera characteristics. Errors in the plume top altitudes can occur when the plume is tilting from or towards the cameras. For a few cases, plume top heights in the vicinity of the volcano were also inferred from stereographic observations made with the Multi-angle Imaging Spectro-Radiometer (MISR) (http://www-misr.jpl. nasa.gov/) and thermal emissivity observations from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER). Furthermore, a few observations were made from a research aircraft (Schumann et al., 2011). Differences between estimated plume heights from the various data sources are substantial (Fig. 6), which is partly due to strong temporal variability of eruption column heights. There are also biases between the various data sets, for example webcam 2 gives consistently higher plume tops than webcam 1 (notice, however, that webcam 1 cannot observe plume tops above 4.8 km). The a posteriori ash emissions are broadly consistent with the observed plume tops. With few exceptions, high plume tops are observed only for periods for which the inversion resulted in relatively large ash emissions, and the modeled plume tops are within the range of observed values. In particular, there is often very close agreement between the modeled and webcam 2 plume tops. Furthermore, in general, the a posteriori plume heights are in better agreement with the webcam observations than the a priori plume heights.
One period with very different a priori and a posteriori emissions is 12-13 May, for which eruption intensity was reported to have declined slightly from previous days with plume tops at 4-5 km . Dense low clouds obscured the eruption plume from the ground on 12 May, making reports of visually observed plume tops inher-ently uncertain. The inversion adjusted both the emission rate and height to higher values (Fig. 2), which is confirmed by plume top heights from independent satellite data (5.5 km) and a short period of webcam observations unaffected by clouds (6-7 km). Due to the emission changes on 12-13 May, the a posteriori simulation reproduces an ash cloud observed by satellite over Great Britain on 14 and 15 May, which the a priori simulation misses nearly completely (see movies described in the Appendix).

Evaluation with lidar data
For evaluation of the vertical structure of the modeled ash cloud, we used data from the CALIOP (Cloud-Aerosol Lidar with Orthogonal Polarization) lidar on the CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations) platform (Winker et al., 2007), which has a vertical resolution of 30-60 m. We analyzed total attenuated backscatter at 532 nm, which is a primary Level 1B profile product (version 3.01). This signal responds to aerosols (including volcanic ash) as well as water and ice clouds which in many cases can be distinguished.  Volcanic ash emitted between 6-10 May wrapped around an anticyclone centered south of Iceland and formed several ash bands extending throughout much of the Eastern North Atlantic (Fig. 7). The satellite observations show detailed horizontal distributions but also some important differences between the SEVIRI and IASI data sets due to retrieval problems, for instance in the vicinity of the volcano and near 37 • N where ash loadings are uncertain, partly because of clouds. CALIPSO passed over the ash cloud and scanned along a vertical cross-section through several ash bands, allowing an evaluation of the modeled ash vertical distribution for ash aged from a few hours to several days. Since the measured total attenuated backscatter cannot be easily converted to ash mass concentrations, the comparison with the model is qualitative and further complicated by water and ice clouds, which produce similar backscatter signals as ash clouds. Nevertheless, the fresh plume (age ≈4 h) observed emanating from the volcano (near 60 • N) is clearly identifiable as volcanic ash and modeled at the correct position, height and vertical extent (2-4 km). The cross section through aged (≈4 days old) ash from 49-57 • N shows modeled ash and measured backscatter maxima at various heights between 3 and 9 km. Not all features are in perfect agreement but the modeled plume maximum near 51 • N is very close to the observed backscatter maximum. The ash band from 46-50 • N (age ≈40 h) is well simulated, with the ash cloud extending from 3-6 km and tilting upwards to the south. Finally, the ash band near 38 • N below 5 km (age ≈3 days), with faint features tilting northward up to 9 km, is also reproduced by the model.
In the Appendix, we show a comparison between the a priori and a posteriori model results for the overpass discussed above (Fig. A3) and four more comparisons with other CALIPSO overpasses. While water and ice clouds often complicate the comparisons, qualitatively we find that the a posteriori results are in better agreement with the CALIPSO data (Fig. A1-A5). We also evaluated our results against quantitative vertical ash concentration profiles obtained from lidar measurements over Europe (Fig. A6-A8) and find that both the modeled a priori and the a posteriori ash concentrations are relatively similar to the observed concentrations, considering the uncertainties in both the model results and the measurements.

Evaluation with aircraft data
The volcanic ash cloud was also probed by the DLR Falcon research aircraft carrying instrumentation for aerosol and other measurements. We use the ash mass concentrations derived from the measured particle size distribution obtained for the most realistic assumption of the ash particles' optical properties, as described in detail by Schumann et al. (2011). The first research flight was performed over Germany on 19 April 2010 (Fig. 8). At that time, the ash concentrations over Europe were already much lower than from 15-17 April 2010 (see movies described in the Appendix) but most of the European air space was still closed because the ash concentrations over Europe were not known at the time. However, both measured and modeled total ash concentrations along the flight track were considerably lower than the current threshold (<0.2 mg m −3 ) for the "Normal" flying zone (European Commission, 2010). Simulations using either meteorological input data reproduced the observed ash layers near Leipzig, Stuttgart and Munich. Despite model over-and underestimates for sections of the flight, the overall bias is small and within systematic measurement uncertainties (Schumann et al., 2011). The a posteriori mass concentrations are furthermore closer to the measured values than the a priori mass concentrations. Higher ash concentrations were observed by the Falcon on other flights. In particular, extensive ash layers were probed over the North Sea on 17 May (Fig. 9) and over the North Sea and Germany on 18 May (Fig. 10). Besides near Iceland, these two flights detected the highest ash concentrations of the entire airborne campaign (Schumann et al., 2011). The model captured these ash layers and there is relatively good quantitative agreement between the a posteriori model results and the measurements. Plumes of volcanic ash were encountered in 34 events on 10 flights in different regions over Europe and the North Atlantic, allowing a statistical comparison between the model and the measurements. Schumann et al. (2011) reported mean ash concentrations for 12 ash encounters in their Table 3 and referred to the remaining 22 ash encounters in the text. One case was excluded from our analysis, namely the descent into the top part of a dense ash plume relatively close to Iceland on 2 May. The descent was abandoned for safety reasons when high ash concentrations were encountered and, thus, the aircraft spent only three minutes in that part of the plume (Schumann et al., 2011). The model produces orders of magnitude too low ash concentrations at the exact location of the ash encounter. However, this is not surprising given the sampling strategy and, therefore, this case was excluded from the comparison.
The average ash concentration during the remaining 33 plume encounters was 64.3 µg m −3 . The ECMWF-based model sampled along the flight legs with ash encounters produces mean a priori concentrations of 29.5 µg m −3 and mean a posteriori concentrations of 50.9 µg m −3 (Fig. 11), even though, overall, the total a posteriori ash emissions are lower. The improved agreement with the observations comes from the more accurate a posteriori plume positions, which is also reflected by an increase in the Pearson correlation coefficient from an insignificant value of 0.18 to a significant (at better than the 0.05 % level) value of 0.61. Notice that a small underestimation by the model is expected and does not necessarily indicate a model bias. The reason for an expected low bias is that not the entire observation data set could be used because the measurements are also sensitive to water, ice and other large particles. Instead, the data for the comparison were selected by screening the entire observation data set (including gas-phase measurements) for volcanic plumes. An unbiased but imperfect model underestimates the observations in such a comparison. Slight displacements of the modeled plumes relative to the observed plumes lead to sampling lower-concentration parts of the plume in the model as compared to the observations during the observed ash encounters. This is compensated by modeled ash present at other times of the flight, typically in the vicinity of the observed plumes, but these flight sections were not used in the statistics, thus yielding a low bias for the model. This methodological bias could have been avoided only when using the entire measurement data set for comparison, which was not possible. Using the GFS data, there is no improvement in the correlations by the inversion but the negative bias is reduced substantially (Fig. 12).

The area over Europe affected by volcanic ash
Based on the a posteriori simulations, we determine the area over Europe (defined here as the longitude-latitude box 10 • W-30 • E, 36 • N-60 • N) where volcanic ash was present somewhere between the surface and 13 km altitude at concentrations exceeding thresholds set (ex post facto) by the European Commission (2010) (Fig. 13). During three episodes in April and May, volcanic ash concentrations at some altitude in the atmosphere exceeded the limits for the "Normal" flying zone in up to 14 % (6-16 %), 2 % (1-3 %) and 7 % (4-11 %), respectively, of total area over the European domain, which triggered partial closures of European air space. For a limit of 2 mg m −3 only two episodes with fractions of 1.5 % (0.2-2.8 %) and 0.9 % (0.1-1.6 %) occurred, while the current "No-Fly" zone criterion of 4 mg m −3 was rarely exceeded. Clearly, the area affected by volcanic ash depends strongly on the choice of the threshold value. While the European Commission (2010) has set thresholds during the eruption, the International Civil Aviation Organization states in the year 2007 (ICAO, 2007): "Unfortunately, at present there are no agreed values of ash concentration which constitute a hazard to jet aircraft engines." Notice also that the exceedance statistics were derived from our gridded model output. High ash concentrations would occur more frequently with increased model resolution. However, it is unclear whether short exposures of an aircraft to high ash concentrations in very small areas could cause engine damage.

Conclusions
In this paper we have, for the first time, determined the ash emissions of a volcanic eruption as a function of time and altitude. For this, we have used an inverse method which optimally merges a priori information on the emissions based on observed plume heights and an eruption column model, satellite observations of total atmospheric ash columns and sensitivity calculations with a dispersion model. We applied the method to the eruption of Eyjafjallajökull in April-May 2010, which caused severe problems for aviation over Europe because of its extensive ash emissions. We evaluated the model simulations using webcam observations of the eruption column, ground-based and space-based lidar observations, and aircraft observations and found that the inversion consistently improved the quality of the model simulation.
From the inversion, we obtain a total fine ash emission of 8.3±4.2 Tg for particles in the size range of 2.8-28 µm diameter. Extrapolation to the size range of 0.25-250 µm yields a total ash emission of 11.9±5.9 Tg but this value depends on the shape of the assumed size distribution and its uncertainty is difficult to quantify and may be much larger than the assumed 50 %. It is likely that our ash emission rates are lower estimates of true emissions, since some of the ash is deposited locally and never observed by satellite. While the model in principle accounts for this, it ignores processes in the immediate vicinity of the vent (e.g., high turbulence, ash aggregation, local precipitation formation), which will enhance local ash deposition. Thus, our estimates should best be viewed as the ash emissions that are subject to long-range transport.
We find that during three episodes in April and May, volcanic ash concentrations at some altitude in the atmosphere exceeded the limits for the "Normal" flying zone in up to 14 % (6-16 %), 2 % (1-3 %) and 7 % (4-11 %), respectively, of the European area. For a limit of 2 mg m −3 only two episodes with fractions of 1.5 % (0.2-2.8 %) and 0.9 % (0.1-1.6 %) occurred, while the current "No-Fly" zone criterion of 4 mg m −3 was rarely exceeded.
Our methodology has broad applicability. It is efficient enough for real-time application and could supply ash forecasting models (Witham et al., 2007) with an objectively derived quantitative source term. Improved forecasts would then allow more effective emergency response. For retrospective analysis, more accurate knowledge about the spatial distribution of volcanic emissions in the atmosphere would improve the quantification of their radiative (Robock, 2000) and environmental (Durant et al., 2010) impacts. For in-stance, improved estimates of ash deposition into the ocean would allow a better quantification of ocean fertilization , which could be relevant especially for the Icelandic Sea which may be iron-limited (Nielsdóttir et al., 2009).

Appendix A Comparison with space-borne and ground-based lidar data
Figures A1 to A5 show four additional comparisons between cross-sections of CALIPSO total attenuated backscatter and simulated ash concentrations based on both the a priori and a posteriori emissions. Figure A3 is the same case as shown in the main part of the paper (Fig. 7) but here also the a priori model results are shown, allowing to assess the model improvement by the inversion.
The model simulations based on our a posteriori emissions were compared with observations from ground-based lidars in Germany (Leipzig, Munich, Hohenpeissenberg) on 16-17 April. Figures A6-A8   centrations. The left panel shows the a posteriori, the right panel shows the a priori model results. Notice that many features seen in the CALIPSO data (e.g., the extensive layer near 1 km altitude) are water or ice clouds, not ash. 33 Fig. A1. Comparison between satellite-retrieved and modeled ash on 8 May at 04:00 UTC. Top row: Ash total columns retrieved from geosynchronous satellite observations between 04:00-05:00 UTC (left panel) and a posteriori (middle panel) and a priori (right panel) modeled ash total columns. The red line shows the track of CALIPSO overflying the ash plume. Middle row: Vertical cross-sections through the a posteriori (left panel) and a priori (right panel) modeled ash concentrations. Bottom row: Vertical cross-section along the red line shown in top panels of total attenuated backscatter at a wavelength of 532 nm obtained from CALIPSO (colored) superimposed with 20 µg m −3 (black line) and 200 µg m −3 (white line) isolines of modeled volcanic ash concentrations. The left panel shows the a posteriori, the right panel shows the a priori model results. Notice that many features seen in the CALIPSO data (e.g., the extensive layer near 1 km altitude) are water or ice clouds, not ash.     data. The ash cloud mostly extends from about 2-6 km altitude, and a clearly simulated maximum occurs over Leipzig on 16 April at 15:00 UTC (Fig. A6b, c) as also observed by the Leipzig lidar ( Fig. 1 in Ansmann et al., 2010). Especially the shape of the distinctly tilted ash cloud is in very good agreement with the lidar-observed ash cloud. Furthermore, for the ECMWF-based simulations the simulated ash concentrations (≈500-1200 µg m −3 ) over Leipzig and Munich (Fig. A7) agree with the mass concentrations retrieved from the lidar observations (Figs. 4 and 5 in Ansmann et al., 2010), where the indicated mass concentrations were on the order of 1000 µg m −3 in the densest ash layer and short-term peak values were close to 1500 µg m −3 . Other reported  maximum mass concentration of 1.1 mg m −3 (0.65 to 1.7 mg m −3 ) are also in agreement with the simulated concentration values.

Movies of ash dispersion
Movies comparing ash total columns retrieved from SEVIRI with FLEXPART a priori and a posteriori simulations are available from http://zardoz.nilu.no/ ∼ sabine/APRIL.gif and http://zardoz.nilu.no/ ∼ sabine/MAY.gif. These movies show that the a posteriori results are always in better agreement with the SEVIRI observations than the a priori observations, although sometimes differences are relatively small, showing that our a priori emissions are also of high quality.