References

An analytical inversion method has been developed to estimate the vertical profile of SO 2 emissions from volcanic eruptions. The method uses satellite-observed total SO 2 columns and an atmospheric transport model (FLEXPART) to exploit the fact that winds change with altitude – thus, the position and shape of the volcanic plume bear information on its emission altitude. The method finds the vertical emission distribution which minimizes the total difference between simulated and observed SO 2 columns while also considering a priori information. We have tested the method with the eruption of Jebel at Tair, Yemen, on 30 September 2007 for which a comprehensive observational data set from various satellite instruments (AIRS, OMI, SEVIRI, CALIPSO) is available. Using satellite data from the first 24 h after the eruption for the inversion, we found an emission maximum near 16 km above sea level (a.s.l.), and secondary maxima near 5, 9, 12 and 14 km a.s.l. 60% of the emission occurred above the tropopause. The emission profile obtained in the inversion was then used to simulate the transport of the plume over the following week. The modeled plume agrees very well with SO 2 total columns observed by OMI, and its altitude agrees with CALIPSO aerosol observations to within 1–2 km. The inversion result is robust against various changes in both the a priori and the observations. Even when using only SEVIRI data from the first 15 h after the eruption, the emission profile was reasonably well estimated. The method is computationally very fast. It is therefore suitable for implementation within an operational environment, such as the Volcanic Ash Advisory Centers, to predict the threat posed by volcanic ash for air traffic. It could also be helpful for assessing the sulfur input into the stratosphere, be it in the context of volcanic processes or also for proposed geo-engineering techniques to counteract global warming.

information on its emission altitude. The method finds the vertical emission distribution which minimizes the total difference between simulated and observed SO 2 columns while also considering a priori information. We have tested the method with the eruption of Jebel at Tair on 30 September 2007 for which a comprehensive observational data set from various satellite instruments (AIRS, OMI, SEVIRI, CALIPSO) is available.
Using satellite data from the first 24 h after the eruption for the inversion, we found an emission maximum near 16 km above sea level (asl), and secondary maxima near 5, 9, 12 and 14 km a.s.l. 60% of the emission occurred above the tropopause. The emission profile obtained in the inversion was then used to simulate the transport of the plume over the following week. The modeled plume agrees very well with SO 2 total 15 columns observed by OMI, and its altitude and width agree mostly within 1-2 km with CALIPSO observations of stratospheric aerosol produced from the SO 2 . The inversion result is robust against various changes in both the a priori and the observations. Even when using only SEVIRI data from the first 15 h after the eruption, the emission profile was reasonably well estimated. The method is computationally very fast. It is therefore

Introduction
Volcanic eruptions release gases (e.g., water vapor, carbon dioxide, sulfur dioxide, hydrogen sulfide, hydrogen chloride) and solid matter into the atmosphere (Robock and Oppenheimer, 2004). Solid matter contained in eruptions comprises a wide range from solidified pieces of lava with masses of many kilogrammes down to fine particles 5 in the submicron range. While large objects fall back to the surface close to the volcano, fine mineral particles, usually referred to as volcanic ash, can remain in the atmosphere for many days.
The injection height of both gases and ash into the atmosphere varies substantially. In eruptions of non-explosive volcanoes, the injection height is on the order of hundreds 10 of meters and is dominated by thermal plume rise, whereas explosive eruptions have a substantial initial exit velocity. The volcanic plume height also depends critically on the ambient atmospheric conditions (Oberhuber et al., 1998), in particular on the thermal stratification, humidity and wind profile. According to Halmer and Schmincke (2003), 80% of the plumes from explosive volcanic eruptions rise higher than 6 km, 60% higher 15 than 10 km, and 20% higher than 15 km. Material from very strong explosive eruptions such as the Pinatubo eruption in 1991 can reach altitudes of more than 30 km (McCormick et al., 1995).
Climatic effects of volcanic eruptions are dominated by sulfur dioxide (SO 2 ) emissions, as SO 2 is converted to sulfate particles which scatter sunlight and increase the Unfortunately, determining the vertical emission distribution of a volcanic plume is a challenge. In principle, it can be obtained from a high-resolution prognostic eruptioncolumn model such as the Active Tracer High Resolution Model (ATHAM), which uses the mass flux of pyroclastic material and the ambient meteorological conditions as boundary conditions (Oberhuber et al., 1998;Textor et al., 2003). However, getting ap- 10 propriate input data to these models can be problematic, especially in real time when little information on the pyroclastic mass flux and other details of the eruption are available. Depths of volcanic ash columns are often estimated by local observers but this is of unknown accuracy and many volcanic eruptions in remote areas are not observed by eye-witnesses at all. Furthermore, plumes containing SO 2 but no ash cannot be 15 seen directly (accompanying cloud features may be visible). The plume height can also be determined using aircraft (Mankin et al., 1992) or ground-based weather radar or lidar (Wang et al., 2008) but such observations are often not available. Satellite instruments, in principle, provide global coverage. The CALIOP (Cloud-Aerosol Lidar with Orthogonal Polarization) lidar on board of the CALIPSO (Cloud-Aerosol Lidar and 20 Infrared Pathfinder Satellite Observations) platform delivers global aerosol information with high vertical resolution  but poor horizontal sampling. Scanning thermal infrared sounders such as IASI (Infrared Atmospheric Sounding Interferometer) can provide some information on the vertical distribution of SO 2 in a volcanic plume but at very coarse resolution (Clerbaux et al., 2007). Most other satellite prod- 25 ucts (e.g., SO 2 retrievals from SEVIRI, AIRS or OMI, see below) have no or little vertical resolution.
In this paper, we derive the vertical emission profile from atmospheric transport patterns as observed from space by platforms delivering column information but no vertical 3764 resolution. At this stage, we consider only SO 2 , which is easier to retrieve than volcanic ash. This is important for geochemical and climate modeling in itself, and it may also serve VAACs as a proxy to volcanic ash. If the wind speed and/or direction changes with height in the atmosphere, and usually it does, the transport pattern derived from satellite imagery can be used to infer the initial emission profile. We do this by calculat-5 ing transport from many height intervals in the supposed eruption column. By applying an inverse modeling framework, we find the emission profile which leads to simulated spatial patterns of SO 2 column values matching best the observed ones over some assimilation time window. We explore this method in some detail for the recent eruption of Jebel at Tair.

The eruption of Jebel at Tair
Jebel at Tair is a stratovolcano -a steep-sloped cone composed of alternating layers of ash, lava, and rocks produced by earlier eruptions. It is a roughly 4-km 2 island at the mouth of the Red Sea (15.5 • N, 41.8 • E), with a crater summit presently 240 m asl but more than 1500 m above the sea floor (see http://en.wikipedia.org/wiki/Jabal al- Tair   15 Island). Explosive eruptions were recorded in the past but for the last 124 years the volcano lay dormant. It erupted again in the afternoon of 30 September 2007(BGVN, 2007) (see also http://earthfromspace.photoglobe.info/spc jebel al tair.html). Several people died in a Yemeniti military base on the island. Photographic evidence exists that a deep cloud formed above the volcano but we are not aware of any estimates of 20 the height of this cloud.
There is little reliable information on the exact time of the eruption at Jebel at Tair. The most reputable source of information is the Bulletin of the Global Volcanism Network (BGVN). The October 2007 report (BGVN, 2007) provides a synopsis of information garnered from scientists, "eye-witness" accounts, news reports and historical Introduction Interactive Discussion reach high into the atmosphere (≈16 km as inferred from the IASI retrieval and SEVIRI 12 µm imagery: see next) and significant amounts of SO 2 (50-100 kt) were injected it is likely that the eruption was energetic and possibly phreatomagmatic with initial vertical velocities of 50 ms −1 or more (Carey and Bursik, 2000). We have used SEVIRI satellite imagery (12 µm channel) to look for indications of 5 the onset time of the eruption. Figure 1 shows the time evolution of the pixel with the smallest value of the 12 µm brightness temperature, which provides a very good indicator of cloud, within a small region centered over the volcanic island. The difference between 10.8 and 12 µm temperatures which is often a very good indicator of cloudiness or anomalies in an image is also shown. The 12 µm temperature drops rapidly 10 after 11:27 UT, when a small cloud is noticeable in the imagery. Temperatures continue to drop reaching values around 195 K at 12:12 UT and then staying close to this value afterwards. The behavior of the temperature record suggests that the initial eruption (there may have been several eruptions) occurred not later than 11:27 UT and that the cloud reached neutral buoyancy between 12:00-12:30 UT. The temperature difference 15 record shows a negative anomaly in the 11:57 UT image. This is often associated with the presence of volcanic ash but in this case is more likely indicative of overshooting and stratospheric penetration as very little ash was observed in any of the infrared or visible satellite imagery from this eruption. Later, the temperature difference rises and becomes strongly positive which is an indication of ice particles in the cloud. The mini-20 mum value of the 12 µm brightness temperature is 192.6 K at 12:57 UT. Figure 2 shows the ECMWF temperature profile at the closest grid point 1 . For an opaque cloud, a temperature of 192.6 K corresponds to either 16.5 or 17.1 km a.s.l., just below or above the cold point in the ECMWF temperature profile (192.1 K at 16.9 km). This indicates that the cloud penetrated the tropopause, found at an altitude of 15.3 km in the ECMWF 25 profile using the WMO standard definition. One hour earlier, at 11:57 UT, the 12 µm temperature is 203.4 K (−70 • C), which suggests cloud top heights of either 19.1 km or 14.4 km. The latter height seems more likely and hence we conclude that at this time the cloud was still undergoing vertical development.
In summary, the satellite observations combined with ECMWF and radiosonde temperature profiles suggest an initial eruption time no later than 11:27 UT, that the plume reached neutral buoyancy no earlier than 12:00 UT, had the minimum value of the 5 12 µm brightness temperature at 12:57 UT and that it reached a height of more than 16 km. It is probable that SO 2 emissions continued at a reduced rate, either in the form of smaller eruptions or by effusive activity. However, for the inverse modeling, we assume that the SO 2 observed by the satellites was all emitted during the main explosive eruption.

Satellite data
The eruption of Jebel at Tair occurred in the mid-afternoon of a cloud-free day. Several satellite-borne instruments were able to acquire data for this eruption. Of particular relevance to this study was the availability of very high temporal resolution imagery from the geosynchronous Meteosat Second Generation (MSG) Spin-stabilised Enhanced 15 Visible and Infra-Red Imager (SEVIRI). Data from the Atmospheric Infra-Red Sounder (AIRS) and the Ozone Monitoring Instrument (OMI) were also utilized. CALIPSO was able to detect the thin veil of sulfuric acid aerosol formed from the SO 2 erupted from Jebel at Tair. Measurements from the Microwave Limb Sounder (MLS), the Infrared Atmospheric Sounding Interferometer (IASI) and the second Global Ozone Monitoring 20 Experiment (GOME-2) were also available but we have not made use of these data.

OMI
OMI is an ultra-violet (UV) imaging spectrometer designed principally for measuring global ozone (columns and profiles) but with secondary goals of measuring other trace gases, including SO 2 . It measures solar back-scattered radiation in the UV between Dobson Units (DU; 1 DU=2.69×10 16 molecules cm −2 ) or less (Krotkov et al., 2006) but there is also a variable bias error of 0.2-0.5 DU, which depends on the reflectivity of meteorological clouds.
OMI SO 2 level 2 data products (OMSO2 V002) are available to the research community from the Distributed Active Archive Center (DAAC) maintained by NASA's Goddard 15 Space Flight Center (GSFC) and were used in our analyses. Although 0.25 • resolution data are available, we chose to use the swath product and re-sample these data to a common grid. The data were resampled to a grid of 0. gridded data. For input into the inversion scheme, data were averaged to hourly values. Not all of the SEVIRI measurements are useful for determining SO 2 concentrations and we use the retrieval scheme developed by  which relies only on the infrared measurements made at 6.2, 7.3 and 12 µm. The scheme makes use of the strong anti-symmetric stretch absorption feature of SO 2 centered near 7.34 10 µm. For sufficiently large SO 2 gas concentrations that lie above boundary layer water vapor, the top-of-the-atmosphere brightness temperature will be lowered compared to an atmosphere with no SO 2 gas. The decrease in brightness temperature depends on a variety of factors, including the amount of gas, its height in the atmosphere and the presence of interfering gases (water vapor) and clouds. During the first 48 h of the 15 Jebel at Tair eruption the skies were almost cloud-free and the atmosphere relatively dry. Since the gas cloud seems to have reached 16 km, water vapor effects are minimal and do not contribute a significant error to the retrieval. The retrieval method uses a look-up table based on off-line correlated-k radiative transfer calculations (Prata et al., 2003) to relate the band-averaged transmittance (for the SEVIRI channel centered at 20 7.3 µm) to SO 2 columns. The estimated accuracy for a single ifov is ±6 DU. The precision is difficult to ascertain because the major cause of interference is due to water vapor and this is not estimated during the retrieval process. However, for the cases where the SO 2 is sufficiently high in the atmosphere (above 6 km), the precision depends mostly on instrument noise, on the veracity of the off-line radiative transfer 25 and the opacity of the SO 2 gas cloud. We estimate the precision to be about ±6 DU. Introduction

AIRS
AIRS is a high resolution grating spectrometer operating at infrared wavelengths between 3.7 and 15.4 µm (Chahine et al., 2006) and is housed on the EOS-Aqua polar orbiting satellite with equator crossing times of 13:30 LT and 01:30 LT. AIRS scans a swath of ±49 • from nadir with an ifov of 1.1 • providing nadir pixels with dimensions 5 15×15 km 2 , increasing to 18×40 km 2 at the swath edge. Level 1b AIRS products (L1B-AIRS-IR-Rad-V005) were obtained via anonymous ftp from the GSFC DAAC. These data are supplied as granules of 90 pixels by 135 lines and were re-sampled to the same grid as the OMI and SEVIRI data. The first Aqua overpass at 23:47 UT on 30 September 2007, some 12 h after eruption onset, imaged an incomplete cloud be-10 cause the AIRS swath was insufficiently wide. However, the following day and up until 3 October, AIRS was able to provide good coverage of the SO 2 cloud. Like SEVIRI, AIRS has channels that cover the 7.3 µm SO 2 absorption feature, but at much greater spectral resolution (more than 100 channels). The SO 2 retrieval scheme developed by Prata and Bernardo (2007) was used to determine column abundances. 15 This scheme takes advantage of the high spectral resolution and has a better accuracy and precision than SEVIRI. For this case, with no cloud interference, little water vapor interference and a good estimate of the background reference the accuracy is estimated to be ±3 DU, or twice better than that of SEVIRI.

20
CALIPSO, with the CALIOP lidar on board was launched on 28 April 2006 to study the impact of clouds and aerosols on the Earth's radiation budget and climate (Winker et al., 2003. As part of the A-train, CALIPSO flies at 705 km altitude in a 98 which is in good agreement with the theoretical value of 7×10 −7 m −1 sr −1 (Vaughan et al., 2005). They found an altitude agreement between airborne lidar data and CALIOP 10 profiles to be within one CALIOP range bin (60 m).
We have used the total attenuated backscatter at 532 nm, which is a primary level 1 data product. The attenuated backscatter profile is the calibrated, range-corrected, laser energy normalized, baseline-subtracted lidar return signal (see Hostetler et al., 2006, for more details). Due to the better signal-to-noise ratio we have used night-time 15 transects only. The data were ordered and downloaded via ftp from the NASA Langley Atmospheric Science Data Center (ASCD, see http://eosweb.larc.nasa.gov/). To enhance the visibility of the faint layers originating from the eruption in the CALIOP backscatter plots, we have substracted a smoothed average (3 km vertically and 2-3 • horizontally) from a nearby region without any visible aerosol layer from the attenu-20 ated backscatter profiles. As we intend to use the result for a qualitative comparison between CALIOP profiles and FLEXPART simulations, this approach -contrary to a full retrieval of e.g. backscatter ratios (level 2 data are expected in 2008) -seems appropriate. Furthermore, the data have been median filtered over 300 profiles (ca. 100 km) horizontally and 300 m (5 range bins) vertically to decrease the CALIOP resolution to 25 that of the FLEXPART output.

Measurement comparison
Differences in the SO 2 column amount retrieved from AIRS, OMI and SEVIRI are found to be as much as 30%, with, in this case the IR retrievals giving larger columns than the OMI UV retrievals. Figure 3 shows an example of retrievals from OMI, SEVIRI and AIRS at around 10:30 UT on 1 October 2007, about 23 h after the eruption. The patterns of SO 2 distribution within the cloud are generally similar but there are noticeable differences. On average, AIRS columns are about 20% larger than OMI columns and SEVIRI columns are about 10% larger than OMI. Generally, there is northwestward transport of the SO 2 during the first day, which suggests transport with the southeasterly winds in the upper troposphere and lower stratosphere, which are found in the 10 ECMWF data (Fig. 2). The plume covers a relatively large area, which indicates that SO 2 was injected into the atmosphere at more than one altitude, thus allowing the vertical wind shear to rapidly disperse the plume. Spatial integration of SO 2 columns over the volcanic plume yields the total mass of SO 2 in the atmosphere -shown in Fig. 4 as a function of time. The data suggest an 15 emission of the order of 80 kt SO 2 but the mass varies substantially, both between platforms and with time. For instance, AIRS shows a decrease of SO 2 from 1 to 3 October (as it might be expected due to deposition and chemical conversion of SO 2 ), whereas OMI shows an increase during the same time period and a decrease only after that. The reason for the initial increase observed with OMI is not clear. SEVIRI values fluc-20 tuate strongly, showing first a steep increase, which is probably related to the fact that SEVIRI cannot see the entire SO 2 column when there is a lot of SO 2 , and then a steep decrease, which is related to the limited sensitivity of the retrieval. In the FLEXPART model calculations (for description, see Sect. 4), the total SO 2 mass decreases slowly if no cutoff is used, suggesting a lifetime of a few weeks. It drops steeply with the 25 6 DU threshold used for SEVIRI, demonstrating that the SEVIRI retrieval is not sensitive enough for obtaining useful information later than about 36 h after the eruption. The differences in the patterns and in the total mass retrieved from each of the sensors have repercussions for the inverse modeling. In order to remove at least the biases in the total mass, we normalized AIRS and hourly SEVIRI total masses for the first 24 h to the OMI total mass from the first overpass, which we assume to be most accurate.

Height sensitivity
The infrared retrieval schemes have a significant sensitivity to the height of the SO 2 5 cloud. Error in the retrieval of column abundance arises from error in the assumed height of the SO 2 cloud. For remote sounding within an infrared absorption band, neglecting the effects of clouds and other absorbers, the sensitivity to height can be assessed from the radiative transfer equation, 10 where I ν is the radiance emerging at the top of the atmosphere at wavenumber ν, B is the Planck function, and τ the transmittance. The Planck function varies with wavenumber and temperature, which varies with height. The transmittance depends on the absorber profile (q) and is also a function of wavenumber and height. Changing variables, the weighting function is: Weighting functions for a monochromatic channel situated near 7.3 µm are shown in Fig. 5  perturbed atmosphere peaks close to the location of the SO 2 injection. In practice the 7.3 µm channel is sensitive to the profiles of both water vapor and SO 2 and the transmittance of the atmosphere at this wavelength may be regarded as the product of the transmittances of the individual gases. The weighting function can then be written, where the subscripts i , j represent the two gases, in this case H 2 O and SO 2 . If the atmosphere contained only gas i then τ j =1, ∂τ j ∂z =0 and the weighting function is that due to gas j alone. This means that whenever SO 2 and water vapor are collocated it will be difficult to quantify the SO 2 . The AIRS retrieval uses an off-line radiative transfer model and a least squares estimation that reduces the error due to inaccurate 10 knowledge of the absorber height.
The OMI UV retrievals also have a sensitivity to the height of SO 2 , as explained by Yang et al. (2007). The weighting functions in the OMI retrievals are determined for distinct Umkehr layers with layer base altitudes of (approximately) 0, 5.5, 10.3, 14.7 and 19.1 km. The layer with base altitude 14.7 km is referred to as the 15 KM retrieval 15 and this is the SO 2 product used in our study. Averaging kernels (weighting functions) for the OMI retrieval peak between 5 and 15 km, are quite broad and depend on the surface characteristics as well as the profiles of SO 2 and the presence of clouds (see Fig. 7 of Yang et al., 2007). For this study we have used an averaging kernel for a clear atmosphere as illustrated in Fig. 5.
FLEXPART calculates the trajectories of tracer particles using the mean winds interpolated from the analysis fields plus random motions representing turbulence (Stohl 10 and Thomson, 1999). For moist convective transport, FLEXPART uses the scheme of Emanuel andŽivković-Rothman (1999) as implemented by Forster et al. (2007). Calculations were done for a SO 2 tracer, where dry deposition and reaction with the OH radical were considered as sinks. The dry deposition was calculated with the resistance method (Wesely and Hicks, 1977) using data from Wesely (1989) with updates.
Removal by reaction with OH radicals is a new feature in FLEXPART that uses monthly averaged three-dimensional OH concentration fields taken from a long-term simulation with the GEOS-CHEM model (Bey et al., 2001). Aqueous-phase chemistry reactions were not considered. Our reference inversion (see below) was set up to use data only from the first 24 hours after the eruption and, thus, is not critically sensitive to the 20 removal of SO 2 , since the lifetime of SO 2 in the upper troposphere and lower stratosphere is of the order of a few weeks. However, the removal is more important for the comparison of FLEXPART results with satellite data at later times (several days).
It was assumed that the emissions occurred in the column between the ground and 24 km above ground level (a.g.l.), or some subset of this column, above the volcano. 25 The total height range was divided into 160 layers each 150 m deep. For every layer, a simulation with 1 kg of tracer (carried by 15 0000 particles) released uniformly within the layer was performed. The simulations extended over the four days following the eruption. Concentrations were calculated on the same 0. the satellite data were re-sampled, with 9 layers of 2 km vertical resolution between 4 km and above 22 km agl, a single layer between the surface and 4 km agl, and another layer from 22-50 km a.g.l. Total atmospheric columns were calculated by weighting the concentrations in the 11 model layers with the weighting functions (averaging kernels) shown in Fig. 5 (see Sect. 3.6). These model-derived column values represent source-5 receptor relationships, since they were obtained with a unit mass as source. The actual mass released at each level is determined through the inversion. Following the inversion, a single longer simulation over 14 days was made releasing 3 million particles according to the estimated emission profile. The output from this simulation was produced at higher vertical (500 m) but coarser horizontal (1 • ) resolution. 10 It was used for comparisons with independent data.

Inversion method
The estimation of the SO 2 emission profile (SO 2 sources) is based on the analytic inversion method of Seibert (2000Seibert ( , 2001. It has been improved by allowing an a priori for the unknown sources, a Bayesian formulation considering uncertainties for the prior 15 and the observations, and an iterative algorithm for removing negative emission values. The n=160 unknowns (source elements) are put into a vector x, while the m observed values are put into a vector y o , where the superscript o stands for observations. Modeled values y corresponding to the observations can be calculated as 20 where M is the m×n matrix of source-receptor relationships calculated with FLEXPART. One may expect to obtain x by means of multiple linear regression which minimizes the sum of the squared deviations model-observation. However, with the fine resolution of the source that has been introduced, observations do probably not contain sufficient information to constrain well all elements of the source vector, making the problem ill-conditioned. Therefore, regularization or, in other words, additional information is necessary to obtain a meaningful solution. Seibert (2000) has shown that simple Tikhonov regularization, which constrains the squared length of the solution vector ||x|| 2 ≡ x T x in combination with a further term requiring the solution to be smooth can yield useful solutions for inversions of atmospheric trace substances even without ex-5 plicit prior source values. Implicitly, Tikhonov regularisation means zero prior values. Including an explicit a priori source vector x a , we can write and as an abbreviation Mx ≈ỹ. 10 The inversion is then done minimizing a cost function J=J 1 +J 2 +J 3 with the three contributions J 1 measures the misfit model-observation, J 2 the deviation from the a priori values, and J 3 the deviation from smoothness. σ o is the standard error of the observations, 15 and σ x the vector of standard errors of the a priori values. The operator diag(a) yields a diagonal matrix with the elements of a in the diagonal. D is a tridiagonal matrix with elements on the main diagonal equal to −2 and elements of the diagonals above and below equal to 1 (discrete representation of the second derivative), and ǫ is a regularisation parameter determining the weight of this smoothness constraint compared to 20 the other two terms. The standard errors of the observations could be made specific for each receptor element, as done for the prior source vector. However, here we only specify average standard errors for each of the three satellite data sets used. If the three satellite data sets are used together in one inversion, the first part of the cost function becomes where the index k refers to the three data sets. The above formulation implies normally distributed, uncorrelated errors, a condition that we know to be not fulfilled. Observation errors (also model errors are subsumed in this term) may be correlated with neighboring values, and deviations from the prior sources are likely to be asymmetric, with overestimation being more likely than under-5 estimation as zero is a natural bound. The justification for using this approach is the usual one: the problem becomes much easier to solve, detailed error statistics are unknown anyway, and experience shows that reasonable results can be obtained.
Minimization of J leads to a linear system of equations (LSE) to be solved forx (Menke, 1984): The LSE is solved with the LAPACK 3 driver routine SGESVX, based on LU factorisation with calibration of rows and columns (if necessary) and iterative refinement of the solution.
In the case of the inversion with all the three satellite data sets used together, we  Inaccuracies in model and data will in general cause such a method to find solutions containing unphysical negative emissions. In the linear framework this cannot be prevented directly as positive definiteness is a nonlinear constraint. A possible workaround that has been adopted here is to repeat the inversion after reducing the standard error values for those source vector elements that are negative, thus binding the solution 25 closer to the prior values at these heights. This procedure is iterated until the sum of all negative emissions is less than 1‰ of the sum of the positive emissions. During the iteration, which converges quickly, previously negative source elements may change their sign to positive. In this case, the tightening of the value towards the a priori is reduced. The standard errors are correspondingly recalculated as For the practical application, x a , σ x , σ o and ǫ need to be assigned proper values. Regarding the a priori emissions, Clerbaux et al. (2007) have reported a coarseresolution IASI SO 2 profile in the vicinity of the volcano, which shows a broad upper tropospheric/lower stratospheric maximum. We have taken a similarly shaped profile 10 and a total SO 2 emission of 80 kt as our standard a priori (Fig. 6). The uncertainties are taken as proportional (a factor of 2.5 larger) to the respective emission value, except for the lowest 5 km where we choose a larger uncertainty because of the limited sensitivity of the satellite to SO 2 in the lower troposphere. The magnitude of the uncertainty was determined by trial and error, and was chosen to allow substantial corrections to the 15 initial profile. We tested the sensitivity of the inversion to the a priori emission profile by replacing our standard a priori with a constant emission profile and a zero emission profile (see Fig. 7).
The standard error of the observations σ o should be specified for each receptor element and they should contain not only the measurement error but rather be a stan-20 dard misfit between the observations and the model results. Lacking detailed information, we only specify three standard errors: σ o1 =6 DU for all SEVIRI measurements, σ o2 =3 DU for all AIRS measurements, and σ o3 =2 DU for the OMI measurements. The SEVIRI and AIRS standard errors are the actual measurement uncertainties, assuming that here the relatively large measurement uncertainty dominates the measurement- Interactive Discussion error of only 0.5 DU. Our assumed standard error of 2 DU is four times larger, allowing for some variable biases in the OMI retrieval and, furthermore, assuming that for the comparison with OMI the larger part of the misfit stems from the model simulations and the data re-mapping. The weight of the smoothness condition ǫ was determined subjectively as ten times 5 the average standard error of the a priori values. This value was chosen in order to retain robust fine-scale features of the inversion but remove some of the fine-scale variation in less well constrained parts of the profile.

10
Real-time applications (e.g., in VAACs) require a rapid response to volcanic ash hazards. Thus, emission profiles needed by forecast models should be available as soon as possible after an eruption. Less time-critical studies could take advantage also of observations taken at later times, but errors in the satellite retrievals (relative to the decreasing SO 2 values in the plume) and in the model simulation grow in time. To 15 minimize the impact of such errors on the inversion and to make the inversion a realistic example for a real-time context, we use only data from the first 24 hours after the eruption for our so-called "reference" inversion (we use more data in a sensitivity experiment). During the first 24 hours, hourly data from SEVIRI as well as data from single overpasses of AIRS and OMI -both about 23 hours after the eruption (Fig. 3)

20
-are available. We did not use SEVIRI data from the first nine hours after the eruption, since SEVIRI has problems seeing the entire SO 2 column when there is much SO 2 present; contamination of the retrieval by eruption-induced clouds and particles is also most problematic during the first few hours. As described in Sect. 2, the period of active vertical development of the plume is 25 framed by the eruption time, about 11:30 UT, and the time the coldest cloud top tem-3780 perature was observed, about 13:00 UT. Since the active plume development is not simulated by FLEXPART, it is not clear which starting time within this period is most appropriate for the model. We tested three intervals: 11:30-12:00 UT, 12:00-12:30 UT, and 12:30-13:00 UT, during which particles were released at a constant rate. The inverted vertical emission profile was rather similar for these three intervals but the cost 5 function was minimal and correlation between the model and the measurements was greatest for the last interval, so we consider 12:30-13:00 UT as the optimum release time. Figure 6 shows the results from our reference inversion that used data from SEVIRI, AIRS and OMI, as well as results from inversions that used the data from only one 10 instrument at a time, during the hours 10-24 after the eruption. The reference profile (red line in Fig. 6) shows a strong and highly localized emission peak at about 16 km, and secondary peaks at 14, just below 12 km and at 5 km. Smaller emissions are found up to almost 20 km, resulting in 60% (10%) of the total mass being emitted above the local tropopause at 15.3 km (above the cold point at 16.9 km) as determined from 15 the ECMWF data. The sharp decrease of emissions around the cold point could be recovered well by the method because of the strong change in the winds at this altitude (Fig. 2). It is in excellent agreement with the minimum observed cloud top temperature (see section 2), which also indicates a plume top at 16.5 or 17.1 km a.s.l. The emission maxima are collocated with layers of enhanced stability in the atmospheric temperature 20 profile (see Fig. 2), which is in agreeement with the expectation that detrainment of air from the convective updraft was responsible for the injection of SO 2 into the ambient flow.
The estimated emission profile is remarkably robust. Inversions done separately for all three platforms (SEVIRI green line, OMI orange line, AIRS blue line in Fig. 6) yield 25 results that are generally very similar to our reference result. The largest difference occurs for the inversion using AIRS data, which broadens the 16-km peak obtained with the other data sets and shifts it upward by 1-2 km. Large differences occur also below 6 km, where the SEVIRI results show a broad and strong peak which is much weaker in the AIRS and OMI results. These large differences in the lower troposphere result from the decreasing sensitivity of the SO 2 retrievals with decreasing altitude and, thus, a relatively poor constraint on emissions there. Overall, however, Fig. 6 demonstrates that data from a single platform would have sufficed to obtain an emission profile very similar to our reference profile. The results using the SEVIRI data only are particularly 5 encouraging, since SEVIRI data are available every 15 min and can most easily be used in real time. Encouraged by this, we made another inversion using SEVIRI data only from 10 to 15 h after the eruption (violet line in Fig. 6). Even this profile is reasonably close to our reference profile such that a relatively quick estimate of the emission profile could have been made in a real-time situation. 10 In Fig. 7, the result of an inversion using only OMI data until 4 October is shown (blue line). This results in a generally similar profile but a reduced peak at 16 km and increased emissions at 17-18 km, which is somewhat similar to the result using the AIRS data from the first 24 h (blue line in Fig. 6). This inversion yields a much better agreement with OMI data up to 10 days after the eruption but it is possible that the 15 higher altitude of the emission peak is an artifact of the inversion which compensates for growing errors in the transport simulation. For instance, there might have been lofting of the plume en route, which was not properly simulated by FLEXPART. This could have been assisted by radiative heating caused by SO 2 and ash. However, given that we found no evidence of significant amounts of ash, the heating by SO 2 alone 20 would probably have been less than 2 K/day (see Fig. 4 of Gerstell et al., 1995), even on the first day when SO 2 columns where largest.
All the a posteriori emission profiles deviate strongly from our a priori estimate. To further explore the sensitivity to the a priori profile, we show results of inversions using SEVIRI data with vertically constant (orange line) and zero a priori (red line) emissions. 25 For the constant a priori profile, the total mass is the same as used previously (green line in Fig. 7, repeated from Fig. 6 for convenience). At most altitudes, the results are quite similar and, thus, not very sensitive to changes in the a priori profile. Relatively large changes occur below about 5 km and above 23 km where the results are less well 3782 constrained by the measurements and, thus, are bound tighter towards the prior than at other altitudes.
We also explored the sensitivity of our results to the normalization of the total mass to the OMI total mass. Removing this normalization leads to strong fluctuations of the total mass observed by SEVIRI from hour to hour (see Fig. 4), which in turn produce 5 weaker correlations between observed and simulated SO 2 . Nevertheless, the resulting emission profile (violet line in Fig. 7) is still similar to the normalized case (green line), although the total emitted mass is somewhat reduced.

Comparison with independent OMI data
Next we compare the results of a FLEXPART simulation using the reference emission 10 profile as input with independent OMI data from the period 1 to 6 October. The daily SO 2 maps shown in Fig. 8 are composites of data from several overpasses occurring over a period of a few hours. The FLEXPART results have been sampled in the same way, i.e., at the hours of the OMI overpasses and using the OMI weighting function. On 1 October (Fig. 8a), roughly 23 h after the eruption, the SO 2 cloud already covers 15 a relatively large area to the northwest of Jebel at Tair (see also Fig. 3). There is excellent agreement between OMI and FLEXPART, which is perhaps not surprising because these OMI data were part of the input used in the inversion. Nevertheless, the good agreement shows that the ECMWF winds are compatible with the actual dispersion of the volcanic plume and that FLEXPART can handle the transport situation 20 very well. The western part of the simulated plume is the highest (up to 20 km altitude), whereas the eastern part contains contributions mainly from 10-16 km, with smaller contributions also from lower altitudes.
Over the next 24 h, the plume changes travel direction and heads eastward. On 2 October (Fig. 8b), both OMI and FLEXPART show a filamentary two-tailed plume 25 stretching over more than 25 • longitude. On 3 October (Fig. 8c), the plume already stretches over more than 40 • longitude. This filamentation is due to vertical wind shear, with the eastern part of the plume being located at about 15-16 km asl and the western 3783 Introduction part being located at about 17-18 km asl. FLEXPART still reproduces the overall plume shape well, including the plume's two long tails. However, it appears that FLEXPART has too much SO 2 in the plume's head near 75 • E (at altitudes of about 15-16 km) and too little in the northern tail of the plume -in particular, the OMI maximum near 33 • E is not reproduced. The alternative inversion result which used OMI data until 4 October 5 reproduces the SO 2 distribution much better as it has lower emissions at 15-16 km and higher emissions at 17-18 km (not shown). On 4 October (Fig. 8d), both the observed and simulated plume stretch over more than 70 • longitude. The leading part of the plume has almost reached the eastern seaboard of Asia. The southern tail, located near 10-12 km, starts disappearing in 10 the model and has already nearly disappeared in the observations. On 5 October, the plume's head has reached Japan, and the southern tail has now disappeared in both the model and the measurements. The model strongly overestimates SO 2 in the plume head which may partly be due to a too slow removal of SO 2 in the model in this tropospheric part of the plume -notice that aqueous-phase chemistry was ignored 15 in the model calculations but might have been important in this part of the plume. On 6 October (Fig. 8f), the satellite measurements become scattered, whereas the model still suggests a continuous plume. By this time, the satellite retrievals might have difficulties seeing the full plume, as SO 2 columns in parts of the plume have become quite low and clouds obscure part of the plume. Nevertheless, it seems the 20 model overestimates SO 2 columns south of 35 • N -again, this part of the plume is in the troposphere and conversion to sulfate may have been quicker in reality than in the model. This is supported by CALIPSO observations on 7 October (see later), which show an aerosol cloud at 14 km asl in a region where OMI sees little SO 2 . Overall, the agreement between the model-simulated and the observed plume trans-25 port is quite good, even though not all plume maxima are well captured. The outstanding discrepancy is that the trailing part of the plume -which originated in the model from near 17-18 km -is underestimated, whereas the leading part -which originated in the model from near 15-16 km -is overestimated. Indeed, the alternative inversion 3784 using OMI data until 4 October redistributes the emissions to higher altitudes compared to our reference inversion (Fig. 7). Whether this really indicates an initially higher emission, a self-lofting of the plume to greater altitudes en route, or other errors in the model transport is not clear.

5
Aerosols formed by the conversion of SO 2 to sulfate cause enhanced backscatter. In the following, we compare CALIPSO profiles of total attenuated backscatter at 532 nm with SO 2 concentrations simulated by FLEXPART. The comparison is qualitative as we compare two very different quantities. FLEXPART does account for the oxidation of SO 2 by OH radicals but has no tracer for the oxidation product, sulfate, which causes the backscatter. Sulfate can be removed by precipitation in the troposphere such that tropospheric features found in the FLEXPART results may not always be seen by CALIPSO. Nevertheless, a qualitative comparison of plume features is sufficient for our purpose of evaluating the altitude of the simulated plume. Enhanced backscatter can also be caused by clouds. At the altitudes where we 15 find the volcanic plume (about 14-18 km), there may be ice clouds (cirrus) which can be clearly identified in the CALIPSO data by their depolarization signal and "normally" much higher backscatter. The suspected volcanic plume features are so faint that even though the scattering layers can be identified unambiguously, a beyond-doubt identification of these layers as sulfate aerosol is difficult. However, we can rule out 20 alternative aerosol sources, since the CALIPSO backscatter features are found only in the region where the volcanic plume was observed by the other satellite instruments. In addition, the altitudes of 14-18 km are very seldomly reached by normal convection in the subtropics or middle latitudes, suggesting a violent injection of the aerosol into this height range. We did not attempt to identify aerosol backscatter features at lower Interactive Discussion made after more than one week, even though the dilution of the plume counteracts sulfate formation. CALIPSO starts observing the volcanic plume on 2 October (Fig. 9). The lidar profile cuts through the plume's head and observes a thin veil of enhanced backscatter between 28-30 • N at about 13-14 km, approximately 2 km below the tropopause. The 5 plume simulated by FLEXPART stretches further south and is tilted, reaching down to 10 km near 27 • N. The thickness of the aerosol layer is overestimated by FLEXPART, in qualitative agreement with the overestimation of SO 2 in the plume's head compared to OMI observations in that part of the plume. Still, the approximate plume position is reasonably well captured. 10 We show two other examples where CALIPSO observed two different parts of the plume over the Pacific Ocean on 7 (Fig. 10) and 8 October (Fig. 11) On 7 October (Fig. 10), CALIPSO cut through the leading part of the plume and detected a 1-2 km thick aerosol layer at latitudes of [27][28][29][30][31][32][33][34][35] • N and at about 14-15 km a.s.l., just below the tropopause. The observed aerosol layer is again located at the top of the simulated 15 plume, which is also thicker (extending down to 11 km) and stretches further south (to 24 • N). This is one of the strongest backscatter enhancements seen by CALIPSO in the volcanic plume. Noteworthy is the fact that the OMI retrieval shows very little SO 2 in this part of the plume ( Fig. 8f shows OMI data from one day earlier), probably suggesting that a substantial fraction of the SO 2 was already converted to sulfate at 20 these altitudes. On 8 October (Fig. 11), CALIPSO cut through the trailing part of the plume and found a thin veil of volcanic aerosol near 17 km, 1-2 km above the tropopause, and at latitudes of 24-32 • N. The altitude and the thickness of the simulated plume is in good agreement with the observations but it stretches over a larger latitude range. There In summary, the plume simulated by FLEXPART is thicker and extends further to the 3786 north and south than the regions of enhanced backscatter found in the CALIPSO data. The CALIPSO observations are always located near the simulated plume top. SOME MORE TEXT....

Conclusions
We have developed an inverse modeling technique for estimating the vertical profile of 5 SO 2 emissions from a volcanic eruption, using total column measurements of SO 2 from satellites and a Lagrangian particle dispersion model. The method was applied in a case study of the explosive eruption of Jebel at Tair in the Red Sea on 30 September 2007. The good coverage of the Jebel at Tair event by satellite observations under excellent, almost cloud-free conditions and the subsequent long-range transport made 10 it an ideal test case. Important conclusions from our work are as follows: -From total column measurements of SO 2 by a suite of satellite instruments (AIRS, OMI, SEVIRI), we estimate a total emission of 80 (±20) kt of SO 2 into the atmosphere. The dispersion of the SO 2 plume could be observed by some of these instruments for more than a week. Starting from two days after the eruption, 15 highly resolved vertical profiles of aerosol backscatter (sulfate aerosols are produced from the gaseous SO 2 ) were available from CALIPSO.
-Our reference inversion used total-column data from AIRS, OMI and SEVIRI from the first 24 hours after the eruption and yielded an emission maximum at about 16 km asl, and secondary maxima near 5, 9, 12 and 14 km. According to this 20 inversion, 60% of the mass of SO 2 was injected above the tropopause located at 15.3 km, and 10% above the cold point in the temperature profile located at 16.9 km. The sharp decrease of emissions around the cold point agrees well with the cloud top height of 16.5 or 17. -Sensitivity experiments showed that data from a single platform (either AIRS, OMI or SEVIRI alone) and from the first 24 hours after the eruption would have sufficed to produce an emission profile in good agreement with our reference profile. Even using SEVIRI data only from the hours 10-15 after the eruption gave comparable results. This is particularly important since SEVIRI data are operationally avail-5 able in real time every 15 minutes. Sensitivity experiments have also shown that the results are robust against changes in the a priori emission distribution that was used in the inversion, including an a priori zero emission profile.
-Using the emission profile from the reference inversion, the overall plume dispersion as observed by OMI, including transport first to the northwest, than to the 10 east, creation of a two-tailed elongated plume stretching over several dozen degrees of longitude, and transport across Asia and over the Pacific Ocean, could all be simulated well over the course of about a week. However, quantitatively, the relative SO 2 distribution within the plume was not so well simulated. On different days, CALIPSO observations showed thin veils of stratospheric aerosol that 15 were well collocated with the FLEXPART plume. The observed plume tended to be thinner than the simulated one. An inversion experiment using OMI data until 4 days after the eruption shifted the emission maximum from 16 km to 17-18 km and brought the simulation in closer agreement with both the OMI and the CALIPSO observations. However, this may not actually be due to an emission at 20 higher altitude but may instead compensate for errors in the simulated transport (probably due to radiative heating and self-lofting of the plume) en route.
-Using our method, the emission altitudes of volcanic eruptions can be estimated with great accuracy, thus facilitating the understanding of the climatic impacts of stratospheric SO 2 injections by volcanic eruptions. An improved such understand-ference with the Earth system.
-Our analytical inversion method is computationally very efficient. Once the underlying dispersion model calculations are completed, the inversion only takes a few seconds on a normal personal computer. This makes it suitable also for real-time applications in Volcanic Ash Advisory Centers (VAACs).

5
-Aviation requires information of the volcanic threat at designated flight levels. Our results here offer a great improvement over current practice which advises the closing of the entire airspace from ground level to the uppermost flight level.
Knowing that most of the ash or SO 2 is above a flight level provides an opportunity for aircraft to safely fly below the hazard.

10
Further improvements of our inversion method could include some of the following: -For operational application in VAACs, the procedure could be reimplemented with volcanic ash aerosol mass instead of SO 2 if appropriate observations (e.g., of aerosol optical depth) are available. Alternatively, even when using the emission profile obtained from the inversion for SO 2 , a model simulation including the 15 gravitational settling of aerosol could be done subsequently by assuming that the emission profiles for SO 2 and ash are the same except for a vertically constant conversion factor.
-Since the eruption time (or the time when the emissions were effectively injected into the atmosphere) is often not known accurately, a straightforward extension 20 of our inversion procedure would be to consider several emission intervals. The inversion algorithm could then optimize both the vertical and temporal emission distribution at the same time.
-Considerable improvement of the satellite retrievals could be expected when the actual vertical SO 2 distributions from FLEXPART are used instead of standard 25 profiles. Since the inversion uses the satellite data, an iterative scheme alternating between the satellite retrieval algorithm and the inversion algorithm would be 3789 Introduction needed. Even larger improvements would be possible by assimilating not the retrieved SO 2 columns but satellite radiances. This would require the addition of a radiative transfer scheme to the dispersion model, however.
-VAACs need to track and forecast volcanic ash clouds for several days. Even with a perfect source term, atmospheric transport model output will be subject to 5 growing errors because of errors in the underlying wind fields, interpolation errors, self-heating of the plume, etc. On the other hand, new satellite information becomes available every day or, with SEVIRI, even every 15 minutes. This calls for a data assimilation procedure where the horizontal position of the ash cloud and its vertical mass profile in each grid column is regularly reassessed on the basis Vertical Mass Distribution [t/m] A priori SEVIRI SEVIRI 15h OMI AIRS ALL Fig. 6. Inversion results when using only SEVIRI data (green line), only OMI data (orange line), only AIRS data (blue line), and all data combined (our "reference" case, red line), during the hours 10 to 24 after the eruption. Also shown is an experiment that only used SEVIRI data from hours 10 to 15 after the eruption (violet line). The thick black line shows the a priori profile, the thin black line its assumed uncertainty.   Fig. 6; SEVIRI constant (orange line) is using a flat a priori profile (shown by the thick black line) and uncertainty (thin black line); SEVIRI zero is using a zero emission a priori (red line) with the same flat uncertainty range as in the constant a priori scenario; SEVIRI w/o normalization (violet line) is using the data without normalization of the total mass to the OMI total mass, and OMI to 4 October (blue line) is using OMI data until 4 October 2007. Introduction  Fig. 8. Comparison of SO 2 columns measured by OMI and simulated by FLEXPART using the emission profile from our reference inversion for (a) 1 October 09:00-12:00 UT, (b) 2 October 09:00-12:00 UT, (c) 3 October 7-13 UT, (d) 4 October 04:00-12:00 UT, (e) 5 October 04:00-11:00 UT, (f) 6 October 02:00-07:00 UT. The satellite data are shown by the color shading and the FLEXPART results are shown as isolines for 1 mg m −2 (thick black line) and 30 mg m −2 (thick grey line). Continental outlines are shown by thin red lines. Notice that the individual panels show different regions -axes are labelled with longitudes and latitudes, respectively. The location of the volcano is marked with a red triangle in the first four panels.  Fig. 11. Same as Fig. 9 but for 8 October at 17:00 UT.