Volcanic SO 2 fluxes derived from satellite data: a survey using OMI, GOME-2, IASI and MODIS

Abstract. Sulphur dioxide (SO 2 ) fluxes of active degassing volcanoes are routinely measured with ground-based equipment to characterize and monitor volcanic activity. SO 2 of unmonitored volcanoes or from explosive volcanic eruptions, can be measured with satellites. However, remote-sensing methods based on absorption spectroscopy generally provide integrated amounts of already dispersed plumes of SO 2 and satellite derived flux estimates are rarely reported. Here we review a number of different techniques to derive volcanic SO 2 fluxes using satellite measurements of plumes of SO 2 and investigate the temporal evolution of the total emissions of SO 2 for three very different volcanic events in 2011: Puyehue-Cordon Caulle (Chile), Nyamulagira (DR Congo) and Nabro (Eritrea). High spectral resolution satellite instruments operating both in the ultraviolet-visible (OMI/Aura and GOME-2/MetOp-A) and thermal infrared (IASI/MetOp-A) spectral ranges, and multispectral satellite instruments operating in the thermal infrared (MODIS/Terra-Aqua) are used. We show that satellite data can provide fluxes with a sampling of a day or less (few hours in the best case). Generally the flux results from the different methods are consistent, and we discuss the advantages and weaknesses of each technique. Although the primary objective of this study is the calculation of SO 2 fluxes, it also enables us to assess the consistency of the SO 2 products from the different sensors used.


Importance of SO 2 fluxes in volcano monitoring
Volcanism is the surface expression of internal processes, driven by heat generated in the Earth's interior. During eruptions, solid, liquid and gaseous products are generated. The main driving force behind eruptions is exsolution of gas from magma during decompression, which drives ascent through the Earth's crust. Sulfur dioxide (SO 2 ) is one of the most abundant compounds among volcanic gases (e.g., Le Guern et al., 1982;Symonds et al., 1994;and Oppenheimer et al., 2011, among others). Since this gas is very soluble in water and is thermodynamically unstable at low temperature, the presence of SO 2 in volcanic plumes is characteristic of a high emission temperature. SO 2 plumes are the markers of volcanic activity that may be classified into two main types.
-Explosive activity: rapid exsolution of volcanic gases in the volcanic conduit generates an ensemble of particles (tephra) through fragmentation that is ejected explosively into the atmosphere forming a plume. Heat is derived from the erupted tephra and emitted gases, and atmospheric air is entrained, which increases buoyancy. Additional latent heat may be released as water condenses and freezes in the plume. Volcanic plume heights may reach altitudes well into the stratosphere and the

Published by Copernicus Publications on behalf of the European Geosciences Union.
maximum height depends on the mass flux rate (amount of material released as a function of time), the size distribution of erupted particles, and the local wind field.
-Effusive activity: driven by gas exsolution, although at lower rates than during magma fragmentation, molten rock products (lava flows) are erupted at the surface. This style of activity also results in particle generation but the "particles" (magma clots, vesicular scoria) are usually coarse enough to fall out close to the vent. Magmatic gases released are hot and may still produce a significant plume. Thermal energy is also available from the lava, however, these plumes tend only to reach midtropospheric levels except in some exceptional cases.
The main conclusion from the above is that the height at which magmatic gases are injected is not arbitrary but strongly dependent on the eruption regime. As we will show, the height at which SO 2 is present in the atmosphere also has a large influence on the techniques for detecting the amounts present, hence the height becomes a central issue for the characterization of plumes.
The most useful quantity for volcano monitoring is generally the SO 2 flux (expressed in kg s −1 or kT day −1 ), which, combined with other monitoring and petrographic data, yields unique insights into the dynamics of degassing magma. Changes in SO 2 flux are often used as an eruption precursor (e.g., Olmos et al., 2007;Inguaggiato et al., 2011;and Werner et al., 2011). The SO 2 flux is a marker of important volcanic processes as, in particular: (1) replenishment of the magmatic system with juvenile magma (e.g., Caltabiano et al., 1994;and Daag et al., 1996), characterized by a gradual increase of the measured SO 2 flux over a relatively long time period, (2) volatile exhaustion of a degassing or erupting magma body. These processes are marked by a gradual decrease in the SO 2 flux (e.g., Kazahaya et al., 2004) and, as the gases are the driving force of volcanic eruptions, generally precede the end of an eruptive period, and (3) plugging of the upper conduit system by crystallizing magma (e.g., Fischer et al., 1994). This process promotes gas accumulation under the plug, until its destruction or failure when its resistance is exceeded, which is believed to fuel long lasting "vulcanian" and "strombolian" activity (e.g., Iguchi et al., 2008). However, in certain circumstances, SO 2 measurements are less useful. This is the case when dissolution of the gas occurs in a hydrothermal system located between the magma and the surface (Symonds et al., 2001). This process, also known as scrubbing, can complicate the interpretation of SO 2 flux data, being not related to real magmatic processes. An overview on degassing from volcanoes can be found in Shinohara (2008) and Oppenheimer et al. (2011).
Measuring time series of SO 2 flux and integrating them over time is also crucially important for the quantification of the volcanic contribution to the global sulfur budget in the Earth system. SO 2 is, as a precursor of sulfate aerosols, important for air quality (Chin and Jacob, 1996) and climate (Haywood and Boucher, 2000;Rampino et al., 1988). Stratospheric injection of SO 2 especially can have a dramatic impact on the global climate (Robock, 2000;Solomon et al., 2011;Bourassa et al., 2012). SO 2 fluxes are also key for establishing a total gas emission inventory of volcanoes. Indeed, emission fluxes of all the other volcanogenic compounds into the atmosphere are usually calculated by scaling their concentration ratio [X]/[SO 2 ] to the SO 2 flux (e.g., Aiuppa et al., 2008). This inventory is important in volcano monitoring (e.g., CO 2 at Stromboli, and Aiuppa et al., 2011), for studying the volcanic impact on the local atmospheric chemistry (e.g., Oppenheimer et al., 2010) and on the environment (Delmelle et al., 2002), or for quantifying greenhouse gas emissions from volcanoes.

Measurements of volcanic SO 2 fluxes
Remote measurements of magmatic volatiles have focused a lot on SO 2 because it is arguably the most readily measurable by absorption spectroscopy, due to its low background concentration in the atmosphere (∼ 1 ppbv in clean air; Breeding et al., 1973), and the strong and distinctive structures in its absorption spectrum both at ultraviolet (UV) and thermal infrared (TIR) wavelengths (see the HITRAN database; Rothman et al., 2009). Measurements of volcanic SO 2 fluxes are widely carried out from the surface since the mid-seventies using the Correlation Spectrometer (COSPEC) instrument (Stoiber et al., 1983;Stix et al., 2008). More recently, consolidated measurements of total emission fluxes of SO 2 have also been made possible for a network of active degassing volcanoes (Network for Observations of Volcanic and Atmospheric Change -NOVAC; Galle et al., 2010 and references therein), with the advent of inexpensive and highquality scanning miniaturized differential optical absorption spectroscopy instruments (Mini-DOAS).
For unmonitored volcanoes or explosive volcanic eruptions, space-based measurements of SO 2 are more appropriate. Since 1978, the Total Ozone Mapping Spectrometer (TOMS) and follow-up instruments have been measuring SO 2 in the UV spectral range (Krueger, 1983;Krueger et al., 1995), although with a rather high detection limit, and helped to establish the long-term volcanic input of SO 2 into the high atmosphere (Carn et al., 2003a;Halmer et al., 2002). In the infrared, space-based sounding of SO 2 was also possible with TOVS (Prata et al., 2003), with data going back to 1978. Over the last two decades, SO 2 measurements from spacebased nadir sensors have undergone an appreciable evolution, owing to improved spectral resolution, coverage and spatial resolution (Thomas and Watson, 2010). The following instruments have been used with measurement channels that correspond to the infrared and ultraviolet SO 2 absorption bands: GOME (Eisinger and Burrows, 1998), SCIAMACHY (Afe et al., 2004), OMI (Krotkov et al., 2006), GOME-2 (Rix et al., 2009(Rix et al., , 2012, MODIS (Watson et al., 2004;Corradini et al., 2009), ASTER (Urai, 2004;Campion et al., 2010) (Prata and Kerkman, 2007;Corradini et al., 2009), AIRS (Prata and Bernardo, 2007) and IASI (Clarisse et al., 2008Walker et al., 2012;Carboni et al., 2012). Compared to ground-based measurements, daily satellite measurements provide a comprehensive cartography of volcanic emissions at a global scale, but only the strongest sources are picked up due to the limitations in the ground resolution and/or sensitivity of the current sensors. The standard output of satellite SO 2 retrievals is not a flux but a vertical column (VC). It represents the amount of SO 2 molecules in a column overhead per unit surface area (generally expressed in Dobson unit -1 DU: 2.69 × 10 16 molecules cm −2 ). From the SO 2 VCs, one can easily calculate the total SO 2 mass (knowing the satellite pixel size), which is a quantity useful to investigate volcanic activity. For example, Carn et al. (2008) monitored degassing of volcanoes in Ecuador and Colombia and showed that the corresponding SO 2 masses had some correlation with volcano seismicity and local observation of volcanic activity.
Estimating fluxes from emitted masses is notoriously difficult and implies more than applying (crude) scaling laws to the measured total SO 2 masses. Until now, there have been only a handful of studies addressing the problem of SO 2 flux estimates from space Urai, 2004;Pugnaghi et al., 2006;Campion et al., 2010;Merucci et al., 2011;Boichu et al., 2013;Carn et al., 2013). One difficulty comes from the fact that accurate flux calculations require additional information (not always well characterized) on the transport and loss of SO 2 following its release into the atmosphere. Another important issue arises from the limitations in terms of accuracy of the SO 2 VCs retrieved and subsequently used. Indeed, satellite SO 2 retrievals are not always in good agreement (see e.g., Prata, 2011). This is believed to be due to, for example, the impact of the calculated or assumed SO 2 plume height on the retrievals or to different sensitivities of the UV and IR bands to SO 2 , aerosols and clouds. Note that, in spite of its general interest from a remote-sensing point of view and for the calculations of SO 2 fluxes, a comprehensive intercomparison and validation study has hitherto not yet been carried out.
For all these reasons, strong and long-lasting volcanic eruptions are challenging as they produce a vertically distributed SO 2 mass profile, which is highly variable in time in response to the instantaneous eruptive power. Once emitted, SO 2 is dispersed along different transport pathways depending on injection altitude due to the vertical shear of the horizontal wind. Consequently, the SO 2 plume may be very inhomogeneous in space. Moreover, the SO 2 plume may extend over long distances (> 1000 km) and is then often only partly covered by satellite instruments that do not provide full daily global coverage. An additional difficulty comes from the fact that a multi-layered SO 2 plume is probed by a given satellite instrument with a measurement sensitivity strongly dependent on the altitude. Furthermore, if SO 2 is distributed in different altitude layers, each of them is characterized by different SO 2 loss rates.
The motivation for this collaborative study is an effort to estimate volcanic SO 2 fluxes using satellite measurements of plumes of SO 2 . We make use of the SO 2 products from the high spectral resolution OMI, GOME-2 (UV), IASI (TIR) and multispectral resolution MODIS (TIR) instruments. These are currently used in an automated mode to provide alerts for aviation safety (as a proxy for the presence of volcanic ash) or for volcano monitoring, in the Support to Aviation Control Service (SACS; Brenot et al., 2013), the European Volcano Observatory Space Services (EVOSS; Ferrucci et al., 2013) and the Support to Aviation for Volcanic Ash Avoidance (SAVAA; Prata et al., 2008) projects. We combine and compare four different approaches and investigate the time evolution of the total emissions of SO 2 for three volcanic events (different in type) occurring in 2011: Puyehue-Cordón Caulle, Chile (using IASI), Nyamulagira, DR Congo (using OMI and GOME-2) and Nabro, Eritrea (using IASI, GOME-2 and MODIS).
Although the main objective of this study is the determination of SO 2 fluxes, we also investigate the consistency of the SO 2 products from the different sensors used. Results from the OMI and GOME-2 UV instruments (with similar vertical measurement sensitivities) are compared, as the calculation of fluxes is in principle not affected by spatial resolution issues. Furthermore, SO 2 masses results from GOME-2, IASI -both on MetOp-A -and MODIS are compared for conditions where additional information on the altitude of the plume is available.
In the next section we give an overview of the instruments and algorithms used to retrieve SO 2 vertical columns. In Sect. 3 we describe the methods to derive SO 2 fluxes. Results are presented in Sect. 4 and conclusions are given in Sect. 5.

GOME-2
The second Global Ozone Monitoring Experiment (GOME-2) is a UV/visible spectrometer covering the 240-790 nm wavelength interval with a spectral resolution of 0.2-0.5 nm (Munro et al., 2006). GOME-2 measures the solar radiation backscattered by the atmosphere and reflected from the surface of the Earth in a nadir viewing geometry. A solar spectrum is also measured via a diffuser plate once per day. The GOME-2 instrument is in a sun-synchronous polar orbit on board the Meteorological Operational satellite-A (MetOp-A), launched in October 2006, and has an Equator crossing time of 09:30 LT (local time) on the descending node. The ground pixel size is nearly constant along the orbit and is 80 km ×40 km. The full width of a GOME-2 scanning swath is 1920 km, allowing nearly global coverage every day. In addition to its nominal swath width, GOME-2 performs once every four weeks measurements in a narrow swath mode of 320 km with increased spatial resolution.
The spectral range and resolution of GOME-2 allows the retrieval of a number of absorbing trace gases such as O 3 , NO 2 , SO 2 , H 2 CO, CHOCHO, OClO, H 2 O, BrO, as well as cloud and aerosol parameters. The SO 2 vertical columns are routinely retrieved from GOME-2 UV backscatter measurements of sunlight (Rix et al., 2009(Rix et al., , 2012van Geffen et al., 2008) using the Differential Optical Absorption Spectroscopy (DOAS) technique (Platt and Stutz, 2008). In a first step, the columnar concentration of SO 2 along the effective light path through the atmosphere (the so-called slant column) is determined using a nonlinear spectral fit procedure. Absorption cross sections of atmospheric gases (SO 2 , O 3 and NO 2 ) are adjusted to the log ratio of a measured calibrated earthshine spectrum and a solar (absorption free) spectrum in the wavelength interval from 315 to 326 nm. Further corrections are applied to account for Rayleigh and rotational Raman scattering and various instrumental spectral features. In a second step, a background correction is applied to the data to avoid nonzero columns over regions known to have very low SO 2 and to ensure a geophysical consistency of the results at high solar zenith angles (where a strong interference of the SO 2 and ozone absorption signals lead to negative SO 2 slant columns). In a third step, the SO 2 slant column is converted into a vertical column using a conversion factor (called air mass factor), estimated from simulations of the transfer of the radiation in the atmosphere. It accounts for parameters influencing photon paths: solar zenith angle, instrument viewing angles, surface albedo, atmospheric absorption, scattering on molecules, and clouds. For the air mass factor calculation, a SO 2 concentration vertical profile is also required. As the information on the SO 2 plume height is generally not available at the time of observation, three vertical columns are computed for different hypothetical SO 2 plume heights: 2.5, 6 and 15 km above ground level (see Sect. 2.5). It should also be noted that an important error source exists for high SO 2 column amount (> 50-100 DU), when the SO 2 absorption tends to saturate leading to a general underestimation of the SO 2 columns.

OMI
The Ozone Monitoring Instrument (OMI) is a Dutch/Finnish instrument flying on the AURA satellite of NASA (launched in July 2004) on a sun-synchronous polar orbit with a period of 100 min and an Equator crossing time of about 13:45 LT on the ascending node. OMI is a nadir-viewing imaging spectrograph that measures atmosphere-backscattered sunlight in the ultraviolet-visible range from 270 to 500 nm with a spectral resolution of about 0.5 nm (Levelt et al., 2006). In contrast to GOME-2, operating with a scanning mirror and onedimensional photo diode array detectors, OMI is equipped with two-dimensional CCD (charge-coupled device) detec-tors, recording the complete 270-500 nm spectrum in one direction, and observing the Earth's atmosphere with a 114 • field of view, distributed over 60 discrete viewing angles, perpendicular to the flight direction. The field of view of OMI corresponds to a 2600 km wide spatial swath. In this way complete global coverage is achieved within one day. Depending on the position across the track, the size of an OMI pixel varies from 13 km ×24 km at nadir to 13 km ×128 km for the extreme viewing angles at the edges of the swath. It should be stressed that due to a blockage affecting the nadir viewing port of the sensor, the radiance data of OMI are altered at all wavelengths for some particular viewing directions of OMI. This row anomaly (changing over time) can affect the quality of the products and hence reduce the spatial coverage of the data (see http://www.knmi.nl/ omi/research/product/rowanomaly-background.php). In this study, the OMI row anomaly pixels have been filtered by using the quality flag (QualityFlag PBL, Bit 11) in the OMSO2 L2 files.
From solar back-scattered measured radiances, SO 2 vertical columns are retrieved using the linear fit algorithm (Yang et al., 2007). The technique uses OMI measurements at several discrete UV wavelengths to yield total SO 2 and ozone columns as well as effective reflectivity (see also http://so2.gsfc.nasa.gov and http://satepsanone.nesdis.noaa. gov/pub/OMI/OMISO2/index.html). Radiative transfer calculations account for the observation geometry, Rayleigh and Raman scattering, atmospheric absorption, surface albedo, aerosols and clouds. Similarly to GOME-2 retrieval, the OMI SO 2 columns are provided for a set of layers with assumed center of mass heights at 0.9, 2.5, 7.5 and 17 km.

IASI
The hyperspectral Infrared Atmospheric Sounding Interferometer (IASI) was launched in 2006 onboard MetOp-A Hilton et al., 2012) for meteorological and scientific applications. Global nadir measurements are acquired twice a day (at 09:30 and 21:30 mean local equatorial time) with a small to medium sized footprint (from a 12 km diameter circle at nadir to an ellipse with 20 and 39 km axes at its largest swath). The instrumental characteristics of IASI are state-of-the-art with an uninterrupted spectral coverage from 645 to 2760 cm −1 , apodized spectral resolution of 0.5 cm −1 and instrumental noise mostly below 0.2 K below 2000 cm −1 .
Operational products of IASI distributed by EUMETSAT include temperature and humidity profiles, surface and cloud parameters and selected trace gases. IASI measurements can be used to retrieve a host of trace gases . Some on global daily scale (e.g., CO 2 , H 2 O, CO, O 3 , HCOOH, CH 3 OH, NH 3 , and CH 4 ), others on a more sporadic basis (e.g., HONO, C 2 H 4 , SO 2 , and H 2 S).
The spectral range from IASI covers three SO 2 absorption bands, the ν 1 band in the ∼ 8.5 µm atmospheric window; the Atmos. Chem. Phys., 13, 5945-5968, 2013 www.atmos-chem-phys.net/13/5945/2013/ strong ν 3 band at ∼ 7.3 µm and the combination band ν 1 + ν 3 at ∼ 4 µm. In this study, we make use of the SO 2 columns generated by the algorithm presented in Clarisse et al. (2012). The retrieval combines measured brightness temperature differences between baseline channels and channels in the SO 2 absorption ν 3 band (at 1371.75 cm −1 and 1385 cm −1 ). The vertical sensitivity to SO 2 is affected by water vapor absorption and is limited to the atmospheric layers above 3-5 km height. The conversion of the measured signal into a SO 2 vertical column is performed by making use of a large look-uptable and EUMETSAT operational pressure, temperature and humidity profiles. The algorithm calculates the SO 2 columns for a set of pre-defined hypothetical altitude levels: 5, 7, 10, 13, 16, 19, 25 and 30 km. The main features of the algorithm are a wide applicable total column range (over 4 orders of magnitude, from 0.5 to 5000 DU), low noise level and a low theoretical uncertainty (3-5 %).

MODIS
The Moderate Resolution Imaging Spectroradiometer (MODIS) is a multispectral instrument on board the NASA-Terra and NASA-Aqua polar satellites (Barnes et al., 1998, http://modis.gsfc.nasa.gov). Terra's descending node (from north to south) crosses the Equator in the morning at about 10:30 LT, while the Aqua ascending node (south to north) crosses the Equator at about 13:30 LT. MODIS acquires data in 36 spectral bands from visible to thermal infrared and its spatial resolution varies between 250, 500 and 1000 m. The sensor scans ±55 • across-track resulting in a 2330 km swath and full global coverage every one to two days. As IASI, MODIS covers all the three infrared SO 2 absorption bands, ν 1 , ν 3 and ν 1 + ν 3 . The ν 1 + ν 3 absorption feature (channel 23) lies in a transparent window but is very weak and affected by scattered solar radiation during the day. The strongest ν 3 signature (channel 28) is highly affected by the atmospheric water vapor absorption so that the SO 2 retrieval is accurate only for plumes in the middle troposphere/lower stratosphere. Finally the ν 1 feature (channel 29) lies in a relatively transparent spectral region and can be used to retrieve also the lower tropospheric SO 2 . The latter is, however, highly affected by the volcanic ash absorption and the retrieval is very sensitive to uncertainties on surface temperature and emissivity (Corradini et al., 2009;Kearney et al., 2009).
The SO 2 retrieval scheme for multispectral TIR measurements was first described by Realmuto et al. (1994) for the TIMS airborne spectrometer and was based on a weighted least squares fit procedure using instrumental measured radiances and simulated radiances obtained by varying the SO 2 amount. Further work carried out by several authors allowed the extension of SO 2 retrieval to different satellite sensors such as MODIS (Watson et al., 2004), ASTER (Realmuto et al., 1997;Pugnaghi et al., 2006) and SEVIRI (Corradini et al., 2009). Here we refer to the procedure scheme described by Corradini et al. (2009) in which the SO 2 column amount is retrieved on a pixel by pixel basis using the ν 3 signature and top of the atmosphere (TOA) simulated radiances.
The TOA radiances look up tables are obtained from MODTRAN 4 computations driven by the atmospheric profiles (pressure, temperature and humidity), the plume geometry (plume altitude and thickness) and varying the SO 2 amounts. The volcanic cloud top height has been determined assuming the cloud top temperature of the opaque parts of the cloud to be equal to the ambient temperature and deriving the height from the temperature vertical profile (Prata and Grant, 2001;Corradini et al., 2010).

Altitude-dependent sensitivity
One of the largest errors on the SO 2 columns is due to a poor a priori knowledge of the height of the SO 2 plumes. To assess this effect on the retrievals presented above, it is useful to consider the satellite SO 2 column averaging kernels. The averaging kernel is defined as the derivative of the retrieved vertical column with respect to the partial column profile (Jacobian), and basically reflects the changes in measurement sensitivity along the vertical axis. As an illustration, Fig. 1 shows examples of functions for UV (GOME-2 and OMI retrievals) and thermal IR (IASI and MODIS ν 3 retrievals) wavelengths.
The UV averaging kernel in Fig. 1 shows the typical behavior of the satellite measurement sensitivity to SO 2 located at different altitudes, with almost constant values in the upper troposphere and lower stratosphere and then a decrease below ∼ 10 km, essentially due to increasing Rayleigh scattering in the troposphere. One can see that SO 2 located at the surface will produce a signal an order of magnitude lower than if it was in the stratosphere. For the TIR measurement, the effect of altitude on the measurement sensitivity is modest in the stratosphere but is critical in the troposphere because of the sharp temperature gradient. The maximum value is reached at the tropopause (∼ 16 km), where the thermal contrast between the plume and the background radiation is highest. Below 5 km altitude, almost no sensitivity to SO 2 is observed because of strong water vapor absorption.
It should be emphasized that no conclusions can be drawn from Fig. 1 on the SO 2 detection limit as the averaging kernel calculations are independent of the signal-to-noise ratios (specific to each instrument).

Methods to derive SO 2 fluxes
As evocated previously, turning from SO 2 VCs or masses to flux data is nontrivial. In this section, we review several existing methods; a summary of which can be found in Table 1 (Sect. 5).
All the techniques are based on mass conservation in that a satellite SO 2 image obtained on a given day reflects the budget of SO 2 emissions and losses (through oxidation and/or deposition) since the start of the eruption. Hence, the retrieval of SO 2 fluxes from a sequence of consecutive SO 2 images generally requires a sufficiently good estimate of the loss term, although in some instances the loss of SO 2 may be considered negligible (e.g., close to the volcano; McGonigle et al., 2004). The simplifying assumption we made throughout this work and which is also supported by literature (review in Eatough et al., 1994), is that the kinetics of the SO 2 removal follow a first order law: where c is the SO 2 concentration and k is the reaction rate constant (s −1 ), which is also often expressed as k = 1/τ , where τ is the SO 2 e-folding time. A difficulty of using this assumption is that k is extremely variable, being sensitive to a number of poorly controlled and spatially variable factors such as plume altitude, cloudiness and atmospheric humidity (Eatough et al., 1994). A major factor is the availability of atmospheric oxidants in the gas (OH radical) and aqueous (H 2 O 2 ) phases. Wet and dry deposition, heterogeneous chemistry and sequestration of gases in ice are other important processes influencing the SO 2 loss rate. The reader is referred to the literature (e.g., Rose et al., 1995;Chin and Jacob, 1996;Graf et al., 1997;Jacob, 1999;Chin et al., 2000; for an overview of the physical and chemical processes of SO 2 removal. The literature on SO 2 reactivity in volcanic plumes (review in Oppenheimer et al., 1998) provides rate constants that span three orders of magnitude, from ∼ 10 −4 s −1 for an ash-rich plume in the tropical boundary layer (Rodriguez et al., 2008) to ∼ 10 −7 s −1 for a plume residing in the superdry and cold stratosphere (Read et al., 1993).

Box method
Arguably, one of the simplest methods to derive a volcanic flux is to consider the SO 2 mass contained within a circle or a box whose dimensions correspond to the total distance traveled by the plume in one day (Lopez et al., 2013) and which are usually determined using a trajectory model or radiosonde wind profile. To account for SO 2 loss, an age dependent correction e t/τ is applied, where t is the plume age defined as the ratio of the distance between a given pixel and the volcano and the wind speed. The daily flux is then simply calculated by dividing the mass within the box by 1 day. Note that this method is limited to cases where the age correction is well defined and accurate. It may be problematic at low altitude when the kinetics of the SO 2 reaction is fast and the plume quickly disperses below the detection limit.

Traverse method
The traverse method is derived from the approach used since the 1970s to operate the COSPEC on volcanoes from a mobile platform (Stoiber et al., 1983;Stix et al., 2008). It has been applied by several authors (Urai, 2004;Pugnaghi et al., 2006;Campion et al., 2010;Merucci et al., 2011) to high spatial resolution satellite images (ASTER and MODIS), as it enables one to retrieve fluxes from small plumes which can then be compared more easily to data obtained from groundbased measurements. Let us define traverses as plume cross sections, i.e., planes perpendicular to the surface and intercepting the transport axis of the plume. The general formulation of the SO 2 mass flux through a surface S is given by where c is the SO 2 mass concentration (kg m −3 ), v is the wind vector and n the unit vector normal to the surface S.
Integrating over the vertical and assuming a constant wind field, the SO 2 flux through a traverse can be calculated as where VC i is the column amount, l i the horizontal length of the i-th pixel of the profile, v the wind speed and θ the angle between the profile direction and transport direction. For simple wedge shaped plumes, the plume direction can be detected automatically by seeking the pixels around the volcano that have the maximal column amount, and an automatic traverse definition can also be easily implemented. A different procedure is based on the interpolation between the atmospheric wind direction profile, collected in the area of interest, and the plume altitude. Note that assuming a single wind speed over the whole plume area can lead to significant Atmos. Chem. Phys., 13, 5945-5968, 2013 www.atmos-chem-phys.net/13/5945/2013/ error, especially if the plume has a large dimension. Although not implemented in the present paper, a solution is to extend the approach to use a 2-D wind field gridded to the resolution of the satellite. Merucci et al. (2011) have demonstrated that considering multiple traverses at increasing distances from the volcano, one can reconstruct the flux history up to several hours before the overpass of the satellite, via the simple relation that exists between the transport speed, distance and duration. Further improvement can be obtained by applying a correction to the SO 2 flux for the loss rate. For a profile at a distance d we have Then the flux of SO 2 as a function of time is easily obtained via the relation t = t obs − t between the time elapsed since the start of eruption (t ), the time delay between the start of the eruption and the satellite overpass (t obs ) and the plume It should be noted that the traverse method in its simplest form (Eq. 4) is only applicable under the assumption of passive advection of the gas in a constant wind field.

Delta-M method
The "delta-M" method was applied to SO 2 measurements for the first time by Krueger et al. (1996), and has also been applied more recently to fire emissions by Yurganov et al. (2011). It relies on time series of the SO 2 mass obtained by successive satellite overpasses and on the mass conservation equation: In this equation, M is the total mass contained in the plume, F is the volcanic SO 2 flux and −k · M is the SO 2 loss term that takes into account actual chemical processes (gas phase oxidation, dry and wet deposition), loss from transport, but also the dilution of the outer part of the plume below the detection limit of the satellite. Here we assume the loss rate to be independent of position. Satellite observations provide discrete time series of the total SO 2 mass present in the atmosphere, with a time resolution that depends on the orbiting parameters and swath width of the instruments. The delta-M method inverts Eq. (5) to yield SO 2 fluxes from SO 2 mass time series as provided by the satellite. From this point, two slightly different approaches can be used. The simplest one is to solve analytically the differential equation (Eq. 5) assuming a constant flux over the time interval t between two consecutive mass estimates M i and M i−1 : A prerequisite for this method is that the whole plume must be covered. This might be an issue for very large plumes using data from a single satellite instrument with limited spatial coverage. As these series do contain some uncertainty, the resulting flux curves often display spikes that are likely not related to real source variation. The assumption of a constant k might also lead to systematic error, since it actually varies with a number of intrinsically variable parameters. Smoothing of the resulting flux series may be necessary to make these spikes disappear. A second approach is to fit an analytic function to the mass time series; the time-dependent flux is then straightforwardly obtained by applying Eq. (5) to the fitted curve.
In the following, the function chosen to fit the data has the form where K 1,2,3,4 are the fit parameters. This form was chosen to fit the basic skewed shape of the time series mass curve.
Hereafter, we will refer to the "skewed shape" method to designate this flux calculation technique. Note that fitting a function to a naturally variable dataset is not necessarily easy, and may lead to systematic errors especially for eruptions having a long duration and/or a complex behavior. Applying a lowpass filter to the data is sometimes necessary to have smooth enough time series. A major advantage of the delta-M and skewed shape techniques is that they are completely independent of the wind field, and hence applicable even if the plume is sheared in a complex wind field or if the plume is stagnating around the volcano in a low wind environment. The main drawback of these techniques is that they yield only first order estimates of the fluxes: the time-dependent fluxes are either too smooth or may contain spikes.

Inverse modeling method
The applicability of the methods described above is limited to cases where (a) SO 2 is injected at a single altitude so that errors related to the plume height ( Fig. 1) can be considered reasonably small, and (b) the SO 2 plume has a simple geometry and/or is well covered by the satellite sensors. Clearly, these conditions are not always met.
For multi-layered SO 2 plumes, a more elaborate technique to retrieve SO 2 fluxes has been designed, which exploits atmospheric transport patterns to derive time and height resolved SO 2 emission profiles (see also e.g., Hughes et al., 2012). To achieve this goal, we have followed the inversion approach that has already been applied in the past for the estimation of the vertical profile of SO 2 emission after the shortlasting eruptions of Jebel al Tair  and Kasatochi (Kristiansen et al., 2010), and more recently for the retrieval of time-and height-resolved volcanic ash after the Eyjafjallajökull eruption (Stohl et al., 2011). In this section, we only briefly summarize the concept of the method, the reader is referred to Sects. 4.1 and 4.3 for details on the inversion settings.
The inversion scheme couples satellite SO 2 column data with results from the FLEXPART Lagrangian dispersion model v8.23 (Stohl et al., 1998(Stohl et al., , 2005, publicly available at http://transport.nilu.no/flexpart). FLEXPART is used to simulate the transport of air masses. It is driven by the 3hourly European Centre for Medium-Range Weather Forecasts (ECMWF) wind field data with 1 • × 1 • resolution. Simulations extending over t days are made for a number n of possible emission scenarios, characterized by release intervals in altitude and time (emission cell). The output of the model is given every 3 h on a predefined altitude grid and on a latitude/longitude grid (dependent on the event studied).
For each emission cell, the simulation is performed with a unit mass source assuming no chemical processes (tracer run). The output of the model (expressed in ng m −3 ) is interpolated in space and time to closely match the satellite observations, and then converted to an atmospheric column by integrating the modeled profiles. The SO 2 columns (expressed in DU) are calculated using a simple unit conversion factor and by applying the satellite averaging kernels (Fig. 1) -estimated for each measured pixel -to the concentration profiles before integration. Moreover, a plume age dependent correction e t/τ is applied for each emission altitude to account for SO 2 losses, where t is the time since the release and τ the e-folding time.
From the resulting modeled columns, it is possible to simulate measurements using a simple linear combination: where y is the measurement vector containing the observed SO 2 columns (m × 1), x is the state vector gathering the SO 2 mass in each emission cell (n×1) and K is the forward model matrix calculated with FLEXPART (m × n), describing the physics of the measurements. The inversion of the emissions (x) from Eq. (8) is an illposed problem, and we have used the well-known optimal estimation method (Rodgers, 2000): a Bayesian method for retrieving a most probable state based on measurements, the expected measurements errors, best guesses of the target state vector prior to measurement, and the associated covariance matrix. The inversion is done by minimizing a cost function J : where S ε is the measurement error covariance matrix, x a is the a priori value of x and S a is the associated covariance matrix. The first term of Eq. (9) measures the misfit between model and observation and the second term is the deviation from the a priori values.
Following Rodgers (2000), the retrieved state vectorx is obtained from the maximum a posteriori solution: where • denotes a Schur product (element-by-element multiplication). F is a filter matrix filled (by 0 or 1) in such a way that the emission i at the time t i is constrained only by measurements at t j where t i ≤ t j ≤ t i + t days.
Generally the retrieval yields negative emissions for some grid cells. To avoid these unphysical results, several iterations (< 10) of the retrieval are necessary where S a is adjusted at each iteration for the elements of the vector with negative values (see Eckhardt et al., 2008 for details).
It should be stressed that, compared to the methods described in Sects. 3.1-3.3, the inverse modeling technique is more strongly dependent on accurate wind field data and is subject to possible misfits (e.g., due to the emission time step used).

Case studies
In the following we present the SO 2 flux results obtained with the different techniques presented in Sect. 3 for each of the three volcanoes studied: Puyehue-Cordón Caulle, Nyamulagira and Nabro.

Puyehue-Cordón Caulle
On 4 June 2011, an explosive eruption started at the Puyehue-Cordón Caulle volcanic complex in central Chile (40.59 • S, 72.12 • W; 2236 m a.s.l.), after a few months of increased seismic activity (Smithsonian Institution, 2012). The initial phase of the eruption was highly energetic and produced a sustained column that was bent by the high altitude winds at the tropopause (∼ 12-13 km; see e.g., Klüser et al., 2013). After four days, the eruption gradually decreased in intensity, but continued for several months. The ash and SO 2 plume circled the Earth three times between 50-80 • S (Klüser et al., 2013;Clarisse et al., 2012), before being dispersed below the detection limit of the satellites. We investigate this eruption using IASI, as it was the most successful at measuring SO 2 from this eruption. Indeed, the performances of UV sensors were limited because of low UV radiation level and strong ozone absorption typical for high solar zenith angle conditions (> 65 • ) during austral winter. The first IASI overpass occurred about seven hours after the beginning of the eruption and captured the concentrated initial plume. The next overpasses revealed lower SO 2 values in the vicinity of Puyehue-Cordón Caulle but the early plume could easily be identified moving eastward from the volcano.
As mentioned in Sect. 3, a prerequisite to any SO 2 flux calculation is the evaluation of the SO 2 loss term. For this purpose, it is useful to look at the retrieved total SO 2 mass as a function of date as displayed in Fig. 2 (left plot).
The time series of SO 2 , here assuming a plume height of 13 km, shows an increase in mass in the first three days after the start of the eruption and then exhibits a decrease as a result of the progressive oxidation of SO 2 . After 10 June, Atmos. Chem. Phys., 13, 5945-5968, 2013 www.atmos-chem-phys.net/13/5945/2013/ the eruption was continuing but as a result of their low intensity and altitude these emissions were below the detection limit of IASI retrievals. We fitted an exponential function to the decay observed at the tail of the total SO 2 mass curve resulting in an estimated e-folding lifetime of 6.8 days, typical for a dry atmosphere. Note that as the plume dilutes, parts of the plume will drop below the detection limit, so that the real e-folding lifetime is probably a bit longer than what we estimate here. From the SO 2 mass time series, we then applied the delta-M method (Sect. 3.3) yielding SO 2 fluxes averaged over a time interval t (Eq. 6). The results are displayed in Fig. 2 (right plot) for half-daily and daily mass estimates ( t = 1/2 or 1 day). Note that the flux estimates for the first half-day and first day have been calculated using a time interval of 7 and 19 h, respectively (instead of 12 and 24 h), as it is more in line with the actual time between the start of the eruption and the first observations. One can see in Fig. 2 that sometimes negative fluxes are returned and that the fluxes time series are quite noisy (especially for the half-daily estimates). This is because the method relies on accurate estimates of consecutive SO 2 masses, which can be problematic for large and dispersed plumes. The erratic behavior of the half-day fluxes on 17 and 18 June is for instance due to the method used for a.m./p.m. differentiation which for these days causes portions of the plume to be counted twice. This behavior disappears when looking at the daily flux.
In addition to the simple delta-M method, Puyehue-Cordón Caulle was found to be an ideal case study to test the traverse method (Sect. 3.2) because of the modest SO 2 loss term and because the SO 2 injection occurred in a strong and well constrained wind field leading to steady and almost one-dimensional plumes (see Fig. 3). The analysis was made for a sequence of four consecutive SO 2 images right after the start of the eruption; each time, the traverses were defined for a well-chosen portion of the plume (close to the volcano).
Following the traverse method formulation (Eq. 4), we calculate the SO 2 fluxes using the MERRA wind fields (available at http://disc.sci.gsfc.nasa.gov/giovanni/overview/) extracted at the location of the volcano and the time of observations. It allows us to reconstruct a high-resolution chronology of the SO 2 flux during the eruption . This is shown in Fig. 4, where the results from four consecutive IASI snapshots (4 June p.m. to 6 June a.m.) are displayed.
A striking feature is the fairly good overlap of the flux data for two successive overpasses, in spite of the constant wind field used in the calculation (for each IASI image). This clearly strengthens the adopted approaches. Note that for the first overpass, larger (and also more scattered) values are observed; this is probably due to the impact of ash (high concentrations in the early plume) on the retrieved values.
Compared to the delta-M method (also shown on Fig. 4, for t = 1/2 day), results from the traverse method provide more information on the temporal variations of the emissions, in particular the peak occurring about 6 h after the beginning of the eruption, not seen with the lower time resolution delta-M method. However, on average, the delta-M and traverse methods give values that are consistent with each other: ∼ 180 kT day −1 on the first day then dropping to ∼ 60 kT day −1 on the second day.
In order to consolidate the flux results of Fig. 4, we have also applied the inversion modeling method (Sect. 3.4) to the Puyehue-Cordón Caulle event. Since the injection heights are reasonably well known, we have used a small number of atmospheric layers but worked at a high temporal resolution. More precisely, emission grid cells have been defined as 3 layers of 1 km thickness from 11 to 14 km and 2-hourly intervals for the period 4-7 June 2011, starting at the onset of the eruption. The FLEXPART simulations extended over 3 days after each release period and the results are given on a 0.5 • × 0.5 • output grid. Following the description in Sect. 3.4, a plume age dependent correction is applied using an e-folding time τ of 6.8 days (Fig. 2). The measurement vector y (Eq. 8) is assembled with the SO 2 columns measured by IASI (> 1 DU) for the period between 4 and 9 June 2011 (daytime and nighttime overpasses) and evaluated on the output grid. The vertical column calculation is made assuming a plume height of 13 km. For the a priori emissions, a single value of 5 kT/2 h of SO 2 is used for each emission grid cell to form x a . The a priori covariance matrix S a has been chosen as diagonal with values corresponding to 250 % errors on the a priori emissions. For the observation error covariance matrix S ε , a diagonal matrix has been chosen and we assumed a conservative 50 % + 1 DU uncertainty on the measurements. The inversion results are displayed in Fig. 4 (grey bars), where the a posteriori emissions have been vertically integrated and scaled to yield SO 2 fluxes. Overall, the results of Fig. 4 are very positive as they show that, for the case of Puyehue-Cordón Caulle, it is possible to resolve the volcanic SO 2 emissions with a time step as short as ∼ 2 h. One can see that there is a very good correspondence between the results of the traverse and inversion modeling methods, which strengthens both approaches.
By integrating the fluxes of Fig. 4, we estimate the total SO 2 mass released by Puyehue-Cordón Caulle during the first two days of eruption to be about 200 kT. However, this estimate is a lower bound as it does not include the contribution of SO 2 emitted at lower altitude.
As an illustration of the inversion modeling method, Fig. 5 shows daily SO 2 column maps retrieved from IASI and sim- Fig. 4. SO 2 flux evolution through the first two days of the Puyehue-Cordón Caulle eruption obtained from the traverse method (analysis of four consecutive IASI images) and compared to the delta-M ( t = 1/2 day; pale orange bars) and inversion modeling methods (grey bars). ulated by FLEXPART using the retrieved emissions (K ·x), for 6 and 7 June 2011. As can be seen, a fair agreement is found between IASI and FLEXPART both for the dispersion patterns and the absolute SO 2 column values.

Nyamulagira
Nyamulagira (1.41 • S, 29.20 • E, 3058 m a.s.l.) erupted over 40 times in the last century, and is known to be one of the main sources of volcanic SO 2 worldwide . Based on TOMS measurements, the SO 2 output from Nyamulagira has been monitored during the 1979-2005 period, and a majority of eruptions produced more than 0.8 MT of SO 2 (Bluth and Carn, 2008).
Almost two years after the last eruption of 2-29 January 2010, Nyamulagira started erupting again, on 6 November 2011 (18:45 UTC) with observed lava fountains reaching up to 300 m and with large amounts of SO 2 easily detected from space. The eruption lasted five months, until 15 March 2012.
For several weeks, the plume of SO 2 was monotonically drifting westward from the volcano (as illustrated in the inset in Fig. 6); its altitude was between 4 and 8 km (based on wind direction data from ECMWF).
We investigated the 2011-2012 Nyamulagira eruption for the first month (when the signal was the strongest) using the UV sensors (OMI and GOME-2) as they are the most suitable to monitor this event owing to their good sensitivity down to the surface. A strong signal of SO 2 could also be detected with IASI especially at the beginning of the eruption, but the data are more difficult to use because of increasing errors on the columns with decreasing plume altitude (the plume of SO 2 was low in the atmosphere). . Phys., 13, 5945-5968, 2013 www.atmos-chem-phys.net/13/5945/2013/ Fig. 5. SO 2 columns measured by IASI (evening overpasses) and simulated by FLEXPART using the a posteriori emissions from the inversion for 6 June 2011 (left panels) and 7 June 2011 (right panels). The measured columns are averaged over the 0.5 • × 0.5 • FLEXPART output. The simulated columns are (1) calculated from the SO 2 profiles weighted by the corresponding altitude-dependent measurement sensitivity functions (averaging kernels) and (2) interpolated at the time of observations. The Puyehue-Cordón Caulle volcano is marked by a black triangle. Fig. 6. OMI daily SO 2 mass as function of plume age (see text) averaged for November 2011. The inset shows the monthly averaged SO 2 columns measured by OMI above Nyamulagira for November 2011.

Atmos. Chem
Since the plume altitude is not known exactly, we simply took the available middle tropospheric altitude (center of mass heights of 6 and 7.5 km for GOME-2 and OMI, respectively) for the SO 2 vertical column calculation. This assumption is a source of uncertainty but the effect on the SO 2 columns is modest, with errors lower than 30 % (Fig. 1). In order to isolate the volcanic SO 2 signal from the background noise, an SO 2 column threshold was assigned (VC > VC lim = VC+3·σ VC ), using the mean (VC) and standard deviation (σ VC ) of daily SO 2 VCs retrieved in the background region .
As already mentioned, the SO 2 plumes from Nyamulagira had a rather simple geometry owing to fairly stable meteorological conditions. This situation facilitates the calculation of the loss term (see e.g., Beirle et al., 2011) and of the SO 2 fluxes. The first loss term has been estimated from the decay observed in the total SO 2 mass downwind (westward) of the volcano. This is illustrated in Fig. 6 showing the averaged SO 2 mass for November 2011 as a function of plume age, expressed in hours.
The values are calculated using daily OMI SO 2 columns integrated along traverses perpendicular to the wind vector. The plume age is the ratio of the distance (defined as positive downwind) between a given pixel and the volcano and the wind speed. For each day, the wind vector is derived from the visual analysis of the transport direction of the plume and comparison with the 1 • × 1 • ECMWF 24 h averaged profile (closest to the volcano) of wind direction and speed. From  Fig. 6, the fitting of an exponential function (red curve) to the SO 2 mass data is straightforward and yields a mean SO 2 e-folding time of 22.5 h, in good agreement with the literature (Bluth and Carn, 2008). Moreover, we repeated the same analysis using GOME-2 data and obtained a consistent value of 24.9 h.
Having reasonably accurate values of the SO 2 loss term for both OMI and GOME-2, we can estimate the fluxes using some of the methods described in Sect. 3. We start to apply these methods to the data from OMI as it is the instrument with the most favorable spatial resolution and noise level.
We first used the delta-M method. Similarly to the analysis of the Puyehue-Cordón Caulle eruption, it is applicable to the daily SO 2 mass values (Fig. 7, grey bars) and does not require information on the daily wind fields.
We also tested the skewed shape method (Sect. 3.3). Since it implies the adjustment of a continuous and regular function (Eq. 7), we applied this method to the (smoothed) 5-day moving averages SO 2 mass time series (Fig. 7, black stems with the fitted curve in red).
The resulting fluxes are shown in Fig. 8 for both methods.
Also shown in Fig. 8, are the results from the box (Sect. 3.1) and traverse (Sect. 3.2) methods. For the box method, we selected the days for which the SO 2 plume had a simple geometry, was well covered by the satellite instrument and when the wind speed was not anomalously high or low. We calculated a daily flux (Fig. 8, blue bars) based on the data contained in a circle around the volcano, with a radius changing each day depending on the wind speed. For the traverse method, daily values (Fig. 8, pink squares) were calculated whenever possible by averaging the fluxes over the closest five to ten traverses (excluding the first one because of sub-pixel dilution issues). Note that the reconstruction of short-term variations of the flux, as for Puyehue-Cordón Caulle, was found difficult. Short and variable SO 2 lifetime, as well as a poorly constrained wind field far away from the volcano makes this kind of approach too uncertain in this case.
In general, the retrieved SO 2 fluxes from the different methods are found to be in reasonable agreement, all showing large values (40-230 kT day −1 ) for the first week of the eruption and then a flattening afterwards (10-30 kT day −1 ) reflecting a calmer eruptive regime. In total, we estimate the SO 2 amount produced by Nyamulagira during the first month of eruption to be of about 1 MT.
Compared to the other methods, the delta-M method shows more variable and uncertain results. This is mainly because the plume is often not completely covered by OMI for two consecutive days. Now comparing the results of the skewed shape, box and traverse methods, the differences observed can be explained by small differences in selection of data, details in the settings and uncertainties on the wind speed.
To consolidate the results, we also analyzed the GOME-2 SO 2 data using the same techniques as for OMI. The results are displayed in Fig. 9, together with the OMI fluxes. Only the data from the circle and skewed shape methods are shown for better readability.
A general consistency is found between OMI and GOME-2 results, although GOME-2 has a tendency to produce higher values than OMI (probably because of differences in the treatment of clouds in the retrievals and in meteorological conditions). As these results are in principle independent of Atmos. Chem. Phys., 13, 5945-5968, 2013 www.atmos-chem-phys.net/13/5945/2013/ Fig. 9. Time series of Nyamulagira SO 2 fluxes retrieved from the box (bars) and skewed shape (lines) methods applied to OMI and GOME-2 data.
the pixel size, they constitute however a valuable element of successful satellite-satellite comparison, complementing existing SO 2 flux comparison exercises (e.g., Campion, 2011;and Campion et al., 2012). Finally, it is important to note that, at the beginning of the event, the SO 2 fluxes (Figs. 8,9) could be underestimated because of SO 2 absorption saturation. However, from the inspection of the SO 2 column values, we believe the effect on the fluxes is small, with differences lower than 10 %.
In an attempt to ascertain a reasonable uncertainty on the derived fluxes, we conducted sensitivity tests by varying the different important parameters of the retrieval. Based on these sensitivity tests, a typical error for both OMI and GOME-2 estimates is of about 50 %, and includes the uncertainties on the SO 2 columns, cutoff values used, SO 2 losses and wind fields.

Nabro
Nabro (13.37 • N, 41.70 • E, 2218 m a.s.l.) is located in the African rift zone at the border between Eritrea and Ethiopia. In spite of no reported historical eruptions, Nabro erupted in the evening of 12 June 2011 around 20:30 UTC, expelling a powerful plume composed of enormous amounts of steam and SO 2 , with little amounts of ash. At the beginning of the eruption, SO 2 was massively injected at tropical tropopause levels where fast easterly jet winds (∼ 30 m s −1 ) transported the plume towards the Mediterranean Middle East and, later, to China. Less powerful, the next phase of the eruption injected SO 2 at lower altitudes but still the dispersion of the plume occurred over a large part of the African continent.
The eruption led to a ∼ 20 km-long lava flow in less than ten days (see http://earthobservatory.nasa.gov/). It lasted about one month; remnants of volcanic SO 2 could still be observed by the UV sensors on 17 July 2011. The comparison of satellite UV and TIR SO 2 images revealed inhomogeneous multi-layered SO 2 plumes dispersed over large distances, making the SO 2 fluxes calculation a complex problem to solve. Therefore we have used a more elaborated approach: the inverse modeling method (Sect. 3.4).
Note that Nabro is a difficult event to analyze as -unlike the Puyehue-Cordón Caulle (explosive) and the Nyamulagira (effusive) eruptions discussed above -it fundamentally was an effusive eruption releasing an exceptional gas jet that reached the tropopause prior to strong lava emission. As there is no evidence that this plume was ash laden as would be expected of a high altitude plume, Nabro remains an intriguing, atypical example in this respect. The transport mechanism of the initial Nabro eruption plume is also peculiar. Bourassa et al. (2012) have suggested deep convection associated with the Asian monsoon circulation to explain the strong impact of the eruption on the stratospheric aerosol loading. However, this mechanism is controversial after Vernier et al. (2013) and Fromm et al. (2013) provided satellite evidence for direct injection of the initial Nabro plume into the lower stratosphere.

Inversion settings
The inversion of the time-dependent SO 2 emission profiles is based on the maximum a posteriori solution described in Sect. 3.4 (Eq. 10). The emission grid cells are defined as: 16 layers of 1 km thickness from 2 to 18 km altitude and 64 6hourly intervals for the period 13-28 June 2011, plus an additional profile to cope with the plume for the very first hours of the eruption (12 June -20:30 UTC to 13 June -00:00 UTC). The inversion consists in estimating the SO 2 masses that have been emitted in each of these grid cells. In line with Eq. (8), we define a state vector x that gathers all the n = 1040 unknown emissions.
The FLEXPART simulations extended over 5 days after each release period and the results are given on a 1 • × 1 • output grid. As mentioned in Sect. 3.4, concentrations are interpolated in space and time with the satellite measurements (see below), multiplied by the corresponding SO 2 averaging kernels (to allow direct comparisons with satellite measurements) and vertically integrated to provide columns. For the present study, we have applied a correction for SO 2 loss using a τ value different for each emission layer. From the decay observed in the total SO 2 mass after the end of the first major injection episode (not shown), we derive an e-folding lifetime of 5 days for the high altitude plumes at 15-18 km (Vernier et al., 2013;Fromm et al., 2013). The e-folding time parameterization as a function of altitude was estimated by interpolation between the values at surface level, taken as 2 days, and the one at 15-18 km (5 days). The value of the SO 2 e-folding time at the surface has been varied until the best match was found between the measurements and the simulations. The observations used to build the measurement vector y (Eq. 8) are the SO 2 columns measured by IASI (daytime and nighttime overpasses) and GOME-2 for the period between 13 June and 3 July 2011, evaluated on the output grid. Note however that the GOME-2 measurements for 14 June 2011 have not been considered in the analysis because the instrument was operating in its narrow swath mode on that day.
As a baseline, the calculation is made for SO 2 columns larger than 1 DU and assuming a plume height of 16 km for both sensors. Hence the averaging kernels by definition equal 1 at this altitude.
To solve the inversion problem (Eq. 10), several other variables need to be specified: the a priori state vector (x a ), the error covariance matrices of the a priori (S a ) and of the measurements (S ε ). For the a priori emissions, we Atmos. Chem. Phys., 13, 5945-5968, 2013 www.atmos-chem-phys.net/13/5945/2013/ have estimated first-order SO 2 flux from the skewed shaped method (Sect. 3.3) applied to the global SO 2 total mass time series from GOME-2 (i.e., with a measurement sensitivity to SO 2 down to the surface). For each day, the inferred SO 2 mass was then distributed equally in each emission grid cell to form x a . The a priori covariance matrix S a has been chosen diagonal with values corresponding to 250 % + 2 kT errors on the a priori emissions, except in the layers below 6 km where a larger uncertainty was used (2000 % + 10 kT of the a priori emissions). Based on sensitivity tests, these settings were adopted to avoid an inversion over-constrained by the a priori values because of reduced sensitivity to SO 2 in the lower troposphere. For the observation error covariance matrix S ε , a diagonal matrix has been chosen and we assumed a conservative 50 % + 1 DU uncertainty on the measurements. Note that we tested the sensitivity of the inversion to the observation error covariance matrix S ε and it was found to be reasonably small.

Results
From inspection of the SO 2 column distributions, it appears that the GOME-2 and IASI data are quite similar, but the sequence of SO 2 images also shows differences that can be attributed (at least partly) to plumes of SO 2 at different altitudes. In principle, the inversion scheme exploits the instrumental complementarities to resolve the injection profiles in the altitude range extending from the surface to the lower stratosphere. At the beginning of the eruption, SO 2 is massively injected at high altitudes and the IASI retrieval is found to be optimal, owing to its large column applicability and 12 h revisiting time (as opposed to 24 h for GOME-2). Shortly after the first injection episode, SO 2 is more evenly distributed throughout the atmosphere, conditions for which GOME-2 provides information on the total atmospheric SO 2 column, including the contribution from the lower troposphere (not probed as efficiently by IASI). Note that after the first week of the eruption, the IASI SO 2 signal close to Nabro is quite small, indicating low altitude SO 2 plumes and hence most information on the emissions is brought by the GOME-2 data. Before inspecting the inverted fluxes, let us first compare the SO 2 columns retrieved from GOME-2 or IASI and simulated by FLEXPART using the retrieved emissions. As an illustration, Fig. 10 shows daily SO 2 maps for 15, 16, 18 and 20 June 2011.
Generally the simulated SO 2 vertical columns agree fairly well with the observations. However, it should be noted that the strong SO 2 plume released at the beginning of the eruption, which then moves towards East Asia, is not well reproduced by the simulations. This plume was injected at high altitudes in the first 15 h of eruption. The reasons for this discrepancy will be investigated in Sect. 4.3.3. Nevertheless, the overall good agreement between the observed and modeled SO 2 columns for the other plumes suggests that it is neither Fig. 11. Comparison of SO 2 flux estimates for the eruption of Nabro using different methods. Upper panel: vertically integrated a posteriori emissions from the inversion combining IASI and GOME-2 data (dark grey bars: > 6 km; light + dark grey bars: total flux), flux estimate of the early plume using the box method applied to IASI (orange bar; see text and Fig. 10a), results of the skewedshape method applied to GOME-2 data and used as a priori in the inversion (blue dotted line; z: 16 km), results of the delta-M method applied to IASI data (magenta and red lines for half-daily and daily mass estimates; z: 16 km) and results of the traverse method applied to MODIS (cyan squares). Lower panel: same as upper panel (but with a focus on the first four days of eruption) with flux results of the traverse method for three IASI images at the beginning of the eruption. a serious nor a general problem. From Fig. 10, one can also see that not all plume maxima are well captured. This points to limitations in the approach used (in terms of spatial and time resolution) to capture the complexity and details of the emissions .
The flux estimates are displayed in Fig. 11  vertically (from 2 to 18 km: sum of dark and light grey bars) and in time to yield daily SO 2 fluxes. As an illustration, Fig. 11 (upper panel) also shows the a posteriori fluxes for the 6-18 km altitude range (dark grey bars), the a priori fluxes (Sect. 4.3.1) and the results from other methods (see discussion below). From Fig. 11, several remarks can be made.
1. It is important to mention that, although it helps minimizing the cost function J (Eq. 9), working with an emission time step of 6 h may also be problematic in some cases, because the satellite revisit time is 12 h at best. This is particularly true for plumes at low altitude where the wind speed is modest and the SO 2 loss processes are fast. Consequently, the variability of the inverted emissions should not be strictly interpreted, as two consecutive values are likely too correlated. Therefore, we have recalculated the SO 2 emissions on a daily basis (more stable in time).
2. Accurate simulation of the transport of SO 2 was challenging for the plume released in the first 15 h of eruption. The corresponding inversion results are uncertain (because of misfit effects), and to estimate the SO 2 flux for the early injection we have applied the box method (see Sect. 3.1). Based on the FLEXPART simulations, it was possible to identify what part of the observed plume was released during the first 15 h (the plume in the orange box in Fig. 10a). By using moving boxes for three consecutive IASI overpasses (from 14 June a.m. to 15 June a.m. 2011), we calculated an averaged total SO 2 mass released in the first 15 h and hence an SO 2 flux ( Fig. 11 -orange bar) by applying a factor accounting for the SO 2 decay since the release.
3. In contrast to the fluxes from the skewed shape method (a priori fluxes), which are generally too low, the fluxes from the inversion modeling method (Fig. 11, upper panel) clearly reveal two volcanic regimes. In the first four days of eruption, large values of the fluxes (500-950 kT day −1 ) are observed and SO 2 is vertically distributed throughout the entire atmosphere (2-18 km range), except probably for the very early plume where heights are mostly in the 15-18 km range. On 17 June 2011 substantially lower fluxes (80-250 kT day −1 ) are observed and SO 2 plumes are mostly below 6 km. This second phase of the eruption is mostly effusive (see Sect. 1.1). It should be noted that the retrieved emissions in the lower troposphere for 13-14 June 2011 have higher uncertainties. This is essentially because the GOME-2 instrument was operating in its narrow swath mode on 14 June 2011 and no useful data was available on that day.
4. The delta-M method (Sect. 3.3) was applied to IASI columns (assuming a plume height of 16 km). One can see that the results (displayed in Fig. 11, magenta and red lines) capture nicely the large SO 2 flux at the start of the eruption and then agree well with the fluxes for the layers in the 6-18 km altitude range. After 20 June, the IASI signal becomes small and the delta-M method returns negative values, because of errors in the IASI retrievals and assumed SO 2 e-folding time. For the first two days of eruptions (13-14 June 2011), it was possible to apply the traverse method (Sect. 3.2) to portions of plumes measured by IASI. The results are shown in the lower panel in Fig. 11, which is a zoom on the first four days of the eruption. The results agree fairly well with the flux estimates in the 6-18 km range (where IASI has a good sensitivity to SO 2 ). Moreover, the results from the traverse method enable us to consolidate the flux results at the very beginning of the eruption (early plume).
5. The SO 2 fluxes have also been estimated by applying the traverse method (Sect. 3.2) to MODIS measurements. The pressure, temperature and humidity profiles, needed for the SO 2 retrieval computation (see Sect. 2.4), derive from ECMWF data considering a 3 • × 3 • grid around Nabro. The plume axes direction and the wind speed have been extracted from the ECMWF profiles at the plume altitude. By exploiting the high spatial resolution of MODIS, the mean SO 2 flux has been estimated in the region lying in the first 300 km from the volcanic vent (also to avoid the depletion of SO 2 since its release). Note that the MODIS TIR images for 13 and 16 June show a strong absorption due to condensed water vapor. Such particles absorb in the entire TIR spectral range and this leads to an overestimation of the SO 2 total amount. To account for this effect a correction for the water vapor content has been applied. The procedure is based on the same approach described by Corradini et al. (2009) to correct for ash but adapted here for water vapor optical properties (as input of the radiative transfer code). As can be seen, the MODIS fluxes do agree reasonably well with the estimates presented above.
6. From Fig. 11, we estimate the total SO 2 mass released by Nabro during the first fifteen days of eruption to be of about 4.5 MT.
7. We conducted test retrievals to estimate the uncertainty on the fluxes displayed in Fig. 11. A reasonable value of the total error is of about 50 %.

Transport of the early plume
We now explore the reasons for the differences between the observations and the FLEXPART simulated patterns. As mentioned above, the large SO 2 plume at the beginning of the eruption was poorly constrained by the inversion. This leads to an underestimation of the modeled SO 2 columns as illustrated in Figs. 10 and 12 (left and center plots). There are several possible explanations for this feature (see discussion by Eckhardt et al., 2008). One possibility is that errors on the high altitude wind fields propagate with time making the transport simulation uncertain. To test this hypothesis, we have run FLEXPART using different wind fields; namely the National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) analyses with 1 • × 1 • resolution. It emerged from sensitivity tests that the simulated dispersion patterns were close to the original data. Another option to explain the discrepancy between the observations and the simulations is to consider an effect from the eruption itself. It could be that for the early phase of the eruption the assumption of a passively advected plume is not fulfilled anymore. Discrepancies may arise from radiative heating of the plume (e.g., leading to self-lofting of the plume; Eckhardt et al., 2008), gas trapping in ice or because of water vapor conversion to ice (Textor et al., 2003). To test this assumption, we have made a sensitivity test illustrated in Fig. 12. Using the a posteriori emissions, one can see that the simulated columns strongly deviate from the corresponding IASI SO 2 columns for 14 and 16 June 2011. If the observed discrepancies are due to the action of the plume, it is reasonable to assume that several hours after the start of the eruption (12 June 2011 around 20:30 UTC) the physical and chemical processes that occurred in-situ have stopped, so that the plume can be considered as passive again. As a test, we have run a FLEXPART simulation starting on 14 June 2011 (i.e., more than 30 h after the start of the eruption) using the IASI SO 2 observations as emission cells, the source strengths being modulated by the corresponding measured SO 2 masses. Only the data with latitudes above 20 • N have been kept and the simulations assumed a uniform SO 2 profile between 15 and 18 km (Vernier et al., 2013;Fromm et al., 2013). The results are shown in Fig. 12 (right plot) on 16 June 2011 at the time of IASI measurements (morning overpass). One can see that, compared to the modeled columns using the a posteriori emissions, the simulated columns are much closer to the IASI observations. The fact that the pattern extending from West Mongolia to central China could be reasonably well reproduced indicates that the wind data are compatible with the real dispersion of the plume. This suggests that the differences between the observations and the original FLEXPART simulations are due -at least partly -to the action of the plume released in the first hours of the eruption.
More work is needed to study the transport of the Nabro SO 2 plume and to ascertain the importance of the Asian monsoon circulation in the uplift of this plume (Bourassa et al., 2012).

Intercomparison of GOME-2 and IASI
As shown in Sect. 4.3.2, the inversion modeling method applied to the eruption of Nabro provides good results. However, the technique implicitly assumes that the SO 2 columns from GOME-2 and IASI are compatible with each other, in the sense that differences in the column amount are only due to deviations from the actual plume height (Fig. 1). The aim of this section is to investigate whether this assumption is fulfilled.
In a first step, we used the information on the SO 2 plume heights provided by the a posteriori FLEXPART simulations. For this purpose, we calculate, at each time step and for each grid point, the height of the SO 2 center of mass using the optimized (non-zero) SO 2 concentration profiles. The retrieved SO 2 heights are then used to investigate the consistency of the SO 2 column datasets from GOME-2 and IASI. However, this analysis is only possible if the plume altitudes returned by the simulations are reliable enough. Therefore, we have made comparisons with independent altitude estimations obtained with collocated GOME-2 observations, to assess the accuracy of the modeled plumes heights. Here the effective height of the SO 2 plume is retrieved from the GOME-2 backscattered radiance measurements using a direct fitting approach Yang et al., 2010;Nowlan et al., 2011). An example of comparison in a region close to Nabro for 16 June 2011 is displayed in Fig. 13.
From the GOME-2 and FLEXPART SO 2 plume height images in Fig. 13, we found that both data are able to capture a complex situation with two SO 2 plumes at different altitudes: a plume at 3-5 km transported to the south and another plume at 8-12 km drifting westward of the volcano. More quantitatively, the modeled and observed plume heights are found to be in good agreement, with differences in the 1-2 km range. Although it does not constitute a proper validation, the comparison of Fig. 13 tends to consolidate the inverted emissions and related plume height estimations.
Next we compare the GOME-2 and IASI SO 2 data by accounting for the effect of altitude. This was achieved by interpolating the averaging kernels (estimated for each GOME-2 and IASI measurement; Fig. 1) at the corresponding FLEX-PART SO 2 plume heights. The calculation supplies conversion factors that are then applied to the respective vertical columns from both sensors. For the conditions where no information on the altitude was available, we took the plume height value from the closest grid point to the corresponding measurement. Figure 14 shows time series of SO 2 masses for the first fifteen days of the eruption of Nabro using the recalculated GOME-2 and IASI vertical columns. Fig. 13. Comparison of SO 2 plume height retrieved from GOME-2 spectra using a direct fitting approach (left) and simulated by FLEX-PART using the a posteriori emissions (right) for 16 June 2011 in the vicinity of Nabro.
The upper panel in Fig. 14 shows the GOME-2 total SO 2 masses as a function of time, for different selection criteria on the data. The black curve sums the SO 2 amount for all plume height scenarios from the surface to lower stratospheric layers, while the red curve only considers the data with SO 2 plume heights above 6 km. One can see that the agreement between GOME-2 (red curve) and IASI (blue curve) SO 2 masses is remarkable when the comparison is restricted to plumes in the altitude range where IASI is sensitive (> 6 km). The differences between IASI and GOME-2 masses are most of the time lower than 25 %. To rule out possible compensation effects, we have also compared the GOME-2 and IASI time series of SO 2 masses separately for different altitude ranges: 6-9, 9-12 and above 12 km. The results are displayed in Fig. 14 (lower panel). Again, the GOME-2 and IASI SO 2 masses data are in excellent agreement.
These results strengthen the inversion modeling approach presented in Sects. 4.3.1 and 4.3.2 and also confirm that the largest errors on the UV and TIR SO 2 columns is the error made by assuming a wrong SO 2 altitude.
One can argue that for individual scenes the discrepancies between GOME-2 and IASI can be much higher than what has been presented here, because of saturation issues. Indeed, the GOME-2 DOAS retrieval tends to underestimate the SO 2 amount for columns larger than 50-100 DU, whereas the IASI algorithm is designed to cope with nonlinear effects due to a large SO 2 amount . However, from the inspection of the SO 2 columns (measured by IASI), we believe this translates to small differences for the time series shown in Fig. 14 because this saturation effect only affects a small number of pixels, those with the highest SO 2 concentrations (at the beginning of the time series). Note that the plume may also contain large amounts of aerosols (notably ice) and gases (e.g., water vapor) all absorbing or scattering the photons differently depending on the wavelength. This can also cause differences between the UV and TIR SO 2 data (again mostly at the beginning of the time series), Atmos. Chem. Phys., 13, 5945-5968, 2013 www.atmos-chem-phys.net/13/5945/2013/ Fig. 14. Comparison of SO 2 masses measured by GOME-2 and IASI (morning overpasses) and constrained by the height estimates from the FLEXPART simulations, for the first 15 days of eruption of Nabro. Upper panel: total mass time series evaluated for all SO 2 plumes (GOME-2, black line) and for the plumes in the altitude range (> 6 km) where IASI is sensitive (GOME-2, red line; IASI, blue line). Lower panel: time series of SO 2 mass measured with GOME-2 (red line) and IASI (blue line) for SO 2 plumes in different altitude ranges (from left to right): 6-9, 9-12 and above 12 km.
depending on the local vertical distributions of aerosols and gases. Whether these differences are positive or negative is difficult to ascertain -a fortuitous compensation is not to be excluded. However, a comprehensive intercomparison study dealing with these issues is beyond the scope of this paper.

Conclusions
We have presented four methods for retrieving volcanic SO 2 fluxes, using SO 2 vertical column measurements from satellite nadir-looking instruments. These techniques rely on the mass conservation principle and, hence, a sequence of satellite SO 2 images carries information on the emission chronology. Generally the determination of SO 2 fluxes requires additional information on the transport and loss of SO 2 since its release.
The simplest methods we have used are applicable to eruptions with single injection height and are best suited for SO 2 plumes with simple geometry. For more complex conditions -such as those encountered during long-lasting eruptions characterized by multi-layered SO 2 plumes, e.g., -a more sophisticated technique has been used that combines satellite data and the FLEXPART dispersion model . The methods have been applied to data from UV-visible (GOME-2 and OMI) and thermal infrared (IASI and MODIS) sounders. The three eruptions discussed in this work are very different: Puyehue-Cordón Caulle belongs to the explosive regime, whilst that of Nyamulagira belongs to the effusive regime. We have seen that the third event, that of Nabro (Eritrea), is intriguing as it fundamentally was an effusive one which released -during its first phases -an exceptional gas jet reaching tropopause levels.
Although each method has advantages and weaknesses, overall, the flux results from the different methods are consistent and often complementary. We have shown that the satellite measurements are able to provide SO 2 fluxes with a sampling of a day or less (∼ 2 h in the best case), and with a total error below 50 %. A summary of the techniques and its main features is given in Table 1.
The analysis also enables us to compare properly the SO 2 data from the different sensors used and, in general, the time series of SO 2 flux/mass agree fairly well. More specifically, we have studied the complementarity of the (collocated) GOME-2 and IASI data. The SO 2 mass time series from GOME-2 and IASI are in excellent agreement when comparison is restricted to the altitude range where both instruments are sensitive to SO 2 : as this call for the plume altitude to be the key-factor of the satellite SO 2 data accuracy; further studies on space-borne SO 2 fluxes are expected to benefit from the advent of sensitive algorithms retrieving simultaneously the plume height and the vertical column of SO 2 (Yang et al., 2010;Nowlan et al., 2011;van Gent et al., 2013;Carboni et al., 2012).
Flux estimates will also be improved and better constrained by using (1) sensors with improved noise level, coverage and spatial resolution (e.g., forthcoming TROPOMI instrument onboard the Sentinel-5 Precursor mission; Veefkind et al., 2012), as well as revisiting time; (2) complementary data from GOME-2 and IASI on MetOp-A and MetOp-B platforms (tandem operations); and (3) a consolidated library of SO 2 e-folding times for different conditions (altitude, humidity, etc.) when it becomes available.
More generally, the methods described here are not limited to the retrieval of volcanic SO 2 emissions (point sources) and can also be applied to compute trace gas fluxes, for example, from fire plumes (Yurganov et al., 2011) or during events of long-range transport of pollutants.