A new simplified approach for simultaneous retrieval of SO2 and ash content of tropospheric volcanic clouds: an application to he Mt Etna volcano

A new procedure is presented for simultaneous estimation of SO2 and ash abundance in a volcanic plume, using thermal infrared (TIR) MODIS data. Plume altitude and temperature are the only two input parameters required to run the procedure, while surface emissivity, temperature, atmospheric profiles, ash optical properties, and radiative transfer models are not necessary to perform the atmospheric corrections. The procedure gives the most reliable results when the surface under the plume is uniform, for example above the ocean, but still produces fairly good estimates in more challenging and not easily modelled conditions, such as above land or meteorological cloud layers. The developed approach was tested on the Etna volcano. By linearly interpolating the radiances surrounding a detected volcanic plume, the volcanic plume removal (VPR) procedure described here computes the radiances that would have been measured by the sensor in the absence of a plume, and reconstructs a new image without plume. The new image and the original data allow computation of plume transmittance in the TIR-MODIS bands 29, 31, and 32 (8.6, 11.0 and 12.0 μm) by applying a simplified model consisting of a uniform plume at a fixed altitude and temperature. The transmittances are then refined with a polynomial relationship obtained by means of MODTRAN simulations adapted for the geographical region, ash type, and atmospheric profiles. Bands 31 and 32 are SO 2 transparent and, from their transmittances, the effective ash particle radius ( Re), and aerosol optical depth at 550 nm ( AOD550) are computed. A simple relation between the ash transmittances of bands 31 and 29 is demonstrated and used for SO 2 columnar content ( cs) estimation. Comparing the results of the VPR procedure with MODTRAN simulations for more than 200 000 different cases, the frequency distribution of the differences shows the following: the Re error is less than±0.5 μm in more than 60 % of cases; the AOD550 error is less than ±0.125 in 80 % of cases; thecs error is less than±0.5 g m−2 in more than 60 % of considered cases. The VPR procedure was applied in two case studies of recent eruptions occurring at the Mt Etna volcano, Italy, and successfully compared with the results obtained from the established SO 2 and ash assessments based on look-up tables (LUTs). Assessment of the sensitivity to the plume altitude uncertainty is also made. The VPR procedure is simple, extremely fast, and can be adapted to other ash types and different volcanoes.


Introduction
Worldwide volcanic activity is presently observed with an increasing variety of ground-and space-based instruments to study its influence on the Earth system on both regional and global scales.During volcanic eruptions large volumes of gases and ash can be injected into the atmosphere, forming plumes and dispersed volcanic clouds.These can reside in the atmosphere and be transported downwind, with a lifespan ranging from hours to years depending on the injection height, and with the scales and types of induced effects varying accordingly.The major hazard caused by volcanic ash clouds is for aviation safety (Casadevall, 1984) and timely alerts are needed to mitigate the risk.This kind of threat is monitored by nine Volcanic Ash Advisory Centres (VAACs) around the world, providing advice on the extent and location of ash clouds to local aviation regulators who are responsible for deciding when to impose air space restrictions.The most widely used algorithms for the detection and characterization of volcanic ash clouds are based on the brightness temperature difference (BTD) in the TIR spectral range (Prata, 1989a, b;Wen and Rose, 1994).They provide the spatial distribution of a volcanic ash cloud's total mass, mean effective radius, and aerosol optical depth (AOD) from thermal infrared (TIR) multispectral images.Because of the sporadic nature of volcanic eruptions and their large geographic extents, satellite remote sensing provides the most suitable technique for detecting and assessing volcanic emissions.The detection and assessment of volcanic ash clouds has been performed using different multispectral polar satellite instruments, including the Advanced Very High Resolution Radiometer (AVHRR, Prata, 1989a, b;Wen and Rose, 1994;Krotkov et al., 1999;Yu et al., 2002;Corradini et al., 2010), the Moderate Resolution Imaging Spectroradiometer (MODIS, Hillger at al., 2002;Ellrod et al., 2003;Watson et al., 2004;Corradini et al., 2008Corradini et al., , 2009Corradini et al., , 2010) ) aboard the Terra/Aqua NASA satellites, and multispectral geostationary satellite systems like for example the Geostationary Operational Environmental Satellite (GOES, Rose and Mayberry, 2000;Yu et al., 2002;Ellrod et al., 2003) and the Spinning Enhanced Visible and Infrared Imager radiometer (SEVIRI, Prata and Kerkmann, 2007) on board Meteosat Second Generation (MSG) satellites.
Along with ash, volcanic SO 2 emissions are monitored as an important indicator of volcanic activity (e.g.Allard et al., 1994;Caltabiano et al., 1994), as an ash proxy in the case of collocated SO 2 and ash clouds, and to study the impact on climate (Robock, 2000).The SO 2 algorithm based on TIR multispectral remote sensing data was first described by Realmuto et al. (1994) for the Thermal Infrared Multispectral Scanner (TIMS) airborne sensor, and it exploits the SO 2 absorption feature centred at 8.7 µm.This sensing scheme was then successfully applied to other airborne and satellite-borne multispectral imagers with a similar band in the TIR range and including MODIS, ASTER, and SEVIRI, later being extended to include the 7.3 µm SO 2 absorption feature (Pugnaghi et al., 2002;Watson et al., 2004;Prata and Kerkmann, 2007).
Further recent improvements involved the simultaneous assessment of SO 2 and ash.Since volcanic ash attenuates across the whole TIR atmospheric window (Watson et al., 2004;Corradini et al., 2009;Kearney and Watson, 2009), its presence in a plume interferes with SO 2 assessment based on the 8.7 µm centred band, and to a minor extent also for the 7.3 µm band.For this reason SO 2 assessment must be corrected for the ash contribution when both SO 2 and ash are present in the same volcanic cloud.This happens quite frequently and if the correction is not applied, the SO 2 reading can be grossly overestimated (Corradini et al., 2009;Kearney and Watson, 2009).While commonly used for quantitative assessments of SO 2 and ash clouds, these methods require a number of input parameters such as pressure, temperature and relative humidity atmospheric profiles, surface temperature and emissivity, specific volcanic ash optical properties, plume altitude and thickness, and they make intensive use of radiation transfer codes to evaluate the atmospheric corrections look-up tables (LUTs).Moreover, the results have some limitations due to the assumption of a single atmospheric model over the whole plume, and their reliability depends on well-known conditions that must be verified in the image (Prata et al., 2001).
The volcanic plume removal (VPR) procedure described here was developed in response to the need for a quick and reliable method offering global coverage.It is simple to apply and requires only two input parameters not directly derived from the multi-spectral TIR image itself, namely the average plume altitude and temperature.The procedure was applied to the Etna volcano test case using MODIS multispectral measurements.
The paper is organized as follows.In Sect. 2 the novel SO 2 and ash assessment procedure is described.In Sect. 3 the procedure is tested with the measurement of SO 2 and ash properties from simulated radiances, and in Sect. 4 the results obtained with the VPR procedure are compared to the results of the conventional LUT assessment approach (Corradini et al., 2009) in two Etna eruption case studies.A sensitivity analysis related to the plume altitude is also performed.Section 5 provides a summary of the conclusions.

The volcanic plume removal procedure
The VPR procedure computes the ash AOD, effective radius and mass, as well as the SO 2 column value from multispectral TIR images.After determination of the algorithm input parameters (computed once for each volcano and ash optical property; see below), the assessment only requires knowledge of the volcanic cloud altitude and temperature.
The core of the procedure consists in estimating what a satellite sensor would see if the plume were absent.Then a simplified model of a uniform plume of fixed altitude and temperature is applied.In the present work the multispectral radiometer MODIS on board the Terra and Aqua satellites (Barnes et al., 1998; http://modis.gsfc.nasa.gov/)was used and the procedure was adapted to the Mt Etna volcano (Sicily, Italy), using specific atmospheric profiles of pressure, temperature, and relative humidity (PTH) to compute the MODTRAN simulations needed to estimate the algorithm coefficients (see below).Finally, the present input parameterization of the procedure is valid only for a specific ash type (Volz, 1973).
MODIS observes the Earth's surface multiple times per day and on 36 different channels, from visible (VIS) to thermal infrared (TIR), with a spatial resolution of 1 km at nadir in this region (TIR).Even for those channels that are well within the atmospheric transmission windows, the radiance reaching the sensor is always attenuated by the atmosphere.Therefore, precise atmospheric corrections (Teggi et al., 1999) have to be implemented to detect and evaluate the volcanic plume content.The possibility itself of easily extracting useful information from remotely sensed TIR data without performing complex and time-consuming atmospheric corrections is of great interest for end users.It eliminates the need for ancillary data like local PTH atmospheric profiles and the need to use demanding, complex software like the atmospheric radiative models.This is particularly valuable during volcanic crisis, when a quick response is required.Clearly, these advantages are partially offset by the limited precision of the result due to the approximations introduced to simplify the procedure and reduce the computation time.Nevertheless, all the algorithms of this kind, like Split Window (Prabhakara et al., 1974;McMillin, 1975;Price, 1984), Dual Band (Crisp and Baloga, 1990;Dozier, 1981), and NDVI (Rouse et al., 1973;Roderick et al., 1996), were founded on this differential basis, and their widespread use is due to their simplicity, user friendliness, and speed.The novel VPR procedure presented here was designed with the aim of providing an easier, faster, but nevertheless reliable simultaneous assessment of volcanic ash and SO 2 columnar quantities.It is based on computation of the radiance that the sensor would have measured in the absence of the plume, and then the plume total transmittances in the MODIS bands 29, 31, and 32 (8.6, 11.0, and 12.0 µm).Band 31 and 32 are transparent to SO 2 , while band 29 is affected by both ash and SO 2 (Corradini et al., 2009).From the plume transmittances of bands 31 and 32 the effective ash particle radius (R e ) and the aerosol optical depth at 550 nm (AOD550) are derived.When R e and AOD are known, the ash mass can easily be computed (Wen and Rose, 1994).Band 31 plume transmittance is also used to compute radiance attenuation in band 29, which is caused only by the ash.This allows an estimation of the SO 2 component of the plume transmittance in band 29, and finally the SO 2 columnar abundance.Details of the algorithm are described below.

Radiance at the sensor without a plume
When a volcanic plume is imaged in the TIR, a dip in radiance is typically observed along a line normal to the plume axis.This absorption feature can change its shape and characteristics according to the wavelength considered and the plume composition.
A simple linear trend tangent to the two edges of the absorption feature (i.e. the plume) gives a good estimation of the radiance that would be detected by the sensor if the plume was absent.This is particularly true for a very uniform ocean surface, but works quite well even if the surface underneath the plume is the ground or a uniform meteorological cloud (see below).Vice versa, the estimation is less accurate if the two edges of the absorption feature are above two different surfaces -for example, ground-ocean, ground-cloud, or ocean-cloud.To simplify the following description a common wedge-shaped plume in a steady wind field is considered.
The radiance in the absence of the plume is obtained in three steps: (1) definition of a plume mask; (2) computation of the plume axis and image rotation so that the lines of the image are orthogonal to the plume; (3) obtaining the background radiance by linear interpolation of the radiance measured in the area surrounding the plume.The result of this process is an image similar to the original one but without the plume.
Figure 1a shows the following for the three considered MODIS bands (29, 31, 32) plus band 28 (7.3 µm): (1) the original measured radiance (circles) along a chosen transect normal to the plume axis.This image was measured by MODIS Terra during the Mt Etna eruption on 23 October 2011 at 21:30 UTC; the plume was over the ocean; (2) the thick coloured line represents the background radiance, i.e. the radiance obtained from the interpolation of the radiance measured in the area surrounding the plume.Only the part inside the plume edges (the cyan vertical bars, obtained from the plume mask) has been changed, while outside the edges the original measured radiances are maintained.Therefore, after this operation, the image with the plume removed is obtained; (3) the thin coloured line is the straight-line tangent to the absorption feature (plume); to determine this straight line, only the points between the blue and cyan bars are considered and only the ones satisfying a defined criterion.
Figure 1b shows another scan line normal to the plume of the same MODIS image.However, this time a location is chosen at which the left edge of the plume is over the ocean, while the right edge is over the ground (colder than the sea at 21:30 UTC).All three bands in the atmospheric transmission window (29,31,32) show this effect (difference) between the two edges.It follows that in this case the VPR procedure's results are less accurate.Vice versa, it is not present in band 28 because, due to the marked atmospheric water vapour absorption, this band is practically unaffected by surface characteristics.

The plume transmittance
The model used to determine the total transmittance in each band considered is schematically described in Fig. 2. The plume is located at a constant altitude Z p ; it has a constant thickness and a constant temperature T p , i.e. the temperature of the atmosphere at the altitude Z p .The VPR procedure requires both Z p and T p as input parameters.The plume altitude and temperature are clearly related, but a single variable is not sufficient due to the weather conditions.The plume transmittance is τ p , while τ and τ are the atmospheric transmittances below and above the plume, respectively.The total atmospheric transmittance when the plume is absent is as follows: τ = τ • τ .
In the absence of the plume and ignoring atmospheric scattering, the radiance at the sensor position is where ε is the surface emissivity, B(T s ) is the Planck function computed at the surface temperature T s , L d is the atmospheric downwelling radiance and L u is the atmospheric path or upwelling radiance.All the radiative transfer equation terms are wavelength dependent; this dependence has been omitted for clarity.
In the presence of a volcanic plume, we assume a radiance at the sensor given by As indicated in Fig. 2, τ p is the total plume transmittance (absorption and scattering) and L up is the path radiance due to the plume alone.
In Eq. ( 2) we have done the following: (1) ignored the increase of L d because of downward scattering within the plume; (2) assumed that the whole atmospheric path radiance (L u ) is attenuated by the plume transmittance τ p ; clearly this is only true for the layer below the plume; (3) assumed the layer of atmosphere above the plume to be completely transparent; that is, τ = 1.The effect of the increase in L d is strongly reduced by the high emissivity value of the surface, in particular for a sea surface.The approximations (2) and (3) in part compensate each other and they are increasingly true as the plume altitude increases.For Mt Etna, Eq. ( 2) is more realistic in the case of a volcanic eruption rather than during a simple quiescent degassing.
Both the plume transmittance and the plume path radiance towards the sensor depend on absorption/emission and scattering.It is assumed that where τ pa is the plume transmittance due to the absorption mechanism, while τ ps is the plume transmittance due to scattering by the plume particles (solid and/or liquid).

Absorption only
If only the absorption mechanism is present (e.g. a plume of only absorbing gases), then: τ ps = 1 and τ p = τ pa .In this case the radiance emitted by the plume towards the sensor is where ε p is the plume emissivity and T is the modified plume temperature.This temperature takes into account the different thickness of the layer of atmosphere above the plume From Eqs. ( 2) and ( 4), finally the plume transmittance can be calculated: (5)

Absorption and scattering
If both absorption and scattering mechanisms are present (e.g.plume with gases and aerosol), then τ ps ≤ 1.In this case the upwelling radiance due to the plume is where ε p is the plume emissivity and T the modified plume temperature seen above.From Eqs. ( 2) and ( 6): Clearly, the transmittance of the ash particles due to the scattering τ ps changes according to the slant path, the size of the particles and their concentration.In this procedure the plume transmittance is computed in two steps.

First step
The first step consists in determining an initial raw transmittance.This is computed using only two different values of vertical transmittance due to scattering: where τ psv = 0.965 is the vertical plume transmittance due to scattering.This value is tried first and is valid for optically thick pixels of the plume -that is, with a plume transmittance of τ p ≤ 0.75.The relation µ = 1 cos ϑ, where ϑ is the MODIS viewing angle for the pixel considered, is used to correct for the increase in the path length of the slant column density when compared to the vertical column.For the more or completely transparent pixels of the plume -that is, when a pixel's transmittance computed using τ psv = 0.965 is τ p > 0.75 -then Eq. ( 8) is computed again using τ psv = 0.98.These values of τ psv were empirically chosen using MOD-TRAN simulations.

Second step
The second step refines the transmittance obtained in the first step using Eq. ( 8).The result is adjusted, taking into account the aspects not adequately explained by the simple model applied above.
Using the optical properties of ash computed from the refractive indexes tabulated by Volz (1973), many different scenarios were simulated applying the MODTRAN radiative transfer model.Specifically considered were the following: the 12 monthly mean atmospheric profiles computed from the data recorded over the last 30 yr at Trapani (Sicily), the upper-air WMO station closest to the Mt Etna volcano; 4 plume altitudes (Z p = 4, 6, 8, 10 km) with a constant thickness of 1 km; 11 SO 2 columnar abundances (0 to 10 g m −2 , step 1 g m −2 ); 6 aerosol optical depths at 550 nm (AOD550 = 0.0, 0.078125, 0.15625, 0.3125, 0.625, 1.25); 6 effective radii (R e = 0.785, 1.129, 1.624, 2.336, 3.360, 4.833 µm); 12 viewing angles (0 to 55 degrees, step 5 degrees) for a total of 228 096 simulations for each MODIS channel considered.From all these simulations the radiance at the sensor when the plume is present (L p ) and the radiance without the plume (L 0 ) were computed for both Terra and Aqua MODIS response functions.Then, having established the plume altitude Z p and the air temperature T p at Z p (from the profiles), the plume transmittances for the three considered MODIS bands 29, 31, and 32 are derived using Eq. ( 8).All the transmittances from these procedures were then compared with the original MODTRAN transmittances, and a polynomial relationship was defined to adjust the modelled transmittances to the references (the MODTRAN transmittances).The final plume transmittance is therefore where a n are the fit coefficients derived from the best fit of a third order polynomial to the plot of MODTRAN transmittance vs τ p (see Fig. 3).
A final control is performed on the transmittance after the second step.When the plume is very transparent -that is, when the plume transmittance of band 31, computed with Eq. ( 9), is τ p,31 > 0.95 -then the total transmittance of band 29 (only) is recomputed using Eq. ( 5).In this case practically no ash is assumed.Of course, this generates a small SO 2 abundance overestimation, but the procedure seems to work better.
Table 1 reports the coefficients (a n ) of the cubic polynomial relationship of Eq. ( 9) for the Volz ash particles and for the three MODIS bands considered, estimated using all the included scenarios.The coefficients for the two MODIS radiometers aboard the Terra and Aqua satellites are slightly different because they have different spectral response functions.
Figure 3 shows the Terra MODIS sensor scatter plots (Aqua is the same) of the MODTRAN simulated plume transmittances for the three bands considered, versus the plots computed using the described procedure.The thick straight line represents a slope of unity.The three scatter plots exhibit close agreement (correlation coefficients r 2 greater than 0.995) between the MODTRAN plume transmittances and those derived with the present simplified model.Equation ( 9) shows that once the polynomial coefficients have been estimated, for a given MODIS image the volcanic plume transmittances are computed based on L p , L 0 , and the modified plume temperature (T ) (see Eq. 8).Because L p and L 0 are derived directly from the MODIS image (see Sect. 2.1), the only unknown is T , i.e. the volcanic plume altitude (Z p ) and temperature (T p ) (see Sect. 2.3.1).It is important to highlight that, by recomputing the polynomial coefficients, the VPR procedure can easily be extended to other ash types and applied to different volcanoes.

Plume ash abundance estimation
The total plume transmittance of the MODIS bands 31 and 32 is assumed to be due only to ash absorption and scattering and can be written as where µ takes into account the slant path and the AOD depends on the size and concentration of the ash particles.As can be seen in Fig. 4, which was obtained from MODTRAN simulations, the AOD of band 31 (AOD 31 ) depends linearly on the plume AOD550 (with null offset), but with a slope m 31 which is a function of the particle size.A similar linear relationship exists for band 32: From Eqs. ( 10) and ( 11) it is clear that the ratio of the logarithms for the transmittance values of bands 31 and 32 computed using Eq. ( 9) is equal to the ratio of the slopes and therefore is also a function of the effective particle radius: Figure 5 shows the trend of the ratio m 31 m 32 and also of m 31 versus the effective radius R e .In the illustrated range (0.785-4.833 µm), and for the considered type of ash (Volz), these are respectively decreasing and increasing monotonic functions.Having established m 31 m 32 from Eq. ( 12), it is possible to define R e , and then, from R e , the slope m 31 is also obtained; finally, the AOD550 is computed using Eq. ( 11).
Knowing the particle's effective radius R e and the AOD550 of a pixel, the ash mass is computed using the Wen and Rose (1994) simplified formula: where S is the pixel surface, ρ = 2.6 × 10 3 kg m −3 is the mean density of the ash particles and Q ext (R e ) is the extinction efficiency factor computed at 0.55 µm.This latter is a function of the effective radius R e .

Plume SO 2 abundance estimation
The total plume transmittance in the MODIS band 29 can be thought of as the product of two components: (1) ash, which attenuates the radiance both absorbing and diffusing radiation (here we consider only Volz-type ash, but obviously other ash types or other kinds of aerosol particles can be present in a plume), and (2) SO 2 , with a radiance attenuation effect due only to absorption: where τ pash,29 is the ash particle component, and τ pSO 2 ,29 the SO 2 plume component for the MODIS channel 29.The total  ) is computed using Eqs.( 9) or (5).Therefore, the SO 2 component (τ pSO 2 ,29 ) can be derived from Eq. ( 14) if the particle component (τ pash,29 ) is known.Arvani (2012) showed that for Volztype aerosols, the particle component of the transmittance of band 29 (τ pash,29 ) is highly correlated to that of band 31 (τ p,31 ).
Figure 6 shows this correlation; it was obtained considering, from among the abovementioned MODTRAN simulations, all cases without SO 2 in the plume.A third-degree polynomial was used to fit τ pash,29 from τ p,31 : The b n coefficients, for the two satellites Terra and Aqua, are reported in Table 2. Next, from the total plume transmittance of the MODIS band 29, the strictly SO 2 component of plume transmittance can be computed, and then the SO 2 columnar content abundance can be derived.In fact, the SO 2 component of the plume transmittance of band 29 was assumed to be where, as before, µ takes into account the slant path, c s is the SO 2 columnar content abundance (in g m −2 ), and β (in where T is the previously defined modified plume temperature in K and T 0 = 273.15K.

Assessment of SO 2 and ash properties from simulated radiances
To test the procedure, a comparison was performed between the assumed input data: the effective particle radius R e , the aerosol optical thickness AOD550, the SO 2 columnar content c s , and the results obtained using the procedure.By means of Eq. ( 2) and of MODTRAN simulations, 228 096 radiances at the sensor L p were computed assuming a constant sea surface emissivity ε = 0.98 and a monthly mean sea temperature.This temperature is the monthly average over 30 yr  of data (from NOAA, http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.ersst.html)for a wide Mediterranean Sea region near Mt Etna (latitude from 34 • N to 38 • N; longitude from 14 • E to 18 • E) (Arvani, 2012).For the same cases the radiances at the sensor without the plume (L 0 ) were also computed.Finally, using the assessment procedure described, SO 2 column densities and ash properties (AOD550 and R e ) were derived for the entire set of simulated radiances and compared to the MODTRAN input parameters used to generate these radiances.
Figure 7 shows the distributions of the differences between the procedure results and the input values for the three input variables.It can be seen that more than 60 % of the R e cases exhibit a difference lower than ±0.5 µm, about 80 % of the AOD550 differences are lower than ±0.125, and more than 60 % of the c s differences are lower than ±0.5 g m −2 .

Assessment of SO 2 and ash properties from two MODIS images
The SO 2 and ash masses and fluxes retrieved from the VPR procedure were compared with the results obtained applying a "standard" assessment approach based on comparison between measurements and the look-up tables computed from the radiative transfer model (LUT procedure, Corradini et al., 2009).Two Etna eruption events were considered for this comparison: the first was observed by MODIS-Aqua on 3 December, at 12:10 UTC during the third phase of the 2006 eruption (Andronico et al., 2009a, b), and the second event was observed by MODIS-Terra on 23 October, at 21:30 UTC during the 2011 lava fountain activity (see Fig. 8).The 2006 event was characterized by a low SO 2 plume located at an altitude of about 3.75 km (Merucci et al., 2011).The 2011 event was characterized by a higher plume altitude (about 5.5 km) and the presence of both ash and SO 2 .
The conventional approach uses the brightness temperature difference (BTD) applied to the MODIS channels 31 and 32 (Prata, 1989a, b;Wen and Rose, 1994;Corradini et al., 2008) to establish the ash parameters (AOD550, R e and ash mass).The chi-square procedure, applied to the MODIS channel 29, is used to retrieve SO 2 (Realmuto et al., 1994;Corradini et al., 2009Corradini et al., , 2010)).The LUTs were computed using MODTRAN code driven by the Trapani WMO atmospheric profiles closest to the two events (MODIS images) and Volz ash optical properties.The sea surface temperature was derived from the radiative transfer equation inversion, with a constant surface emissivity equal to 0.99 (Corradini et al., 2008(Corradini et al., , 2009(Corradini et al., , 2010)).The plume altitude for 3 December 2006 was derived from Andronico et al. (2009a, b), while the plume altitude for 23 October 2011 was evaluated from the satellite image by comparing the brightness temperature of the most opaque plume pixels with the atmospheric temperature extracted from the Trapani WMO profile (Prata and Grant, 2001;Corradini et al., 2009Corradini et al., , 2010)).The ash influence on MODIS channel 29 in the SO 2 assessment was corrected following the procedure described by Corradini et al. (2009).

Results
Figure 9 shows the SO 2 and ash maps, as well as fluxes extrapolated from the two MODIS images using the VPR and LUT procedures.The fluxes for the two images were computed by multiplying the columnar abundance along a transect perpendicular to the plume axis by the wind speeds (Pugnaghi et al., 2002;Corradini et al., 2003Corradini et al., , 2009;;Theys et al., 2012).The latter, computed by interpolation from the plume altitude and the Trapani WMO wind speed profiles, were 5 and 12 m s −1 for the 3 December 2006 and 23 October 2011 images, respectively.The maps show the same structures and the fluxes have common trends.In Table 3 the total SO 2 and ash masses and mean fluxes assessed using the two procedures are compared.

Plume altitude sensitivity analysis
The main advantage of the VPR procedure is that SO 2 and ash assessments can be obtained from the remote sensing data themselves, requiring in addition only the mean plume altitude and temperature.Uncertainty in these parameters clearly affects the extrapolation errors.In this section a sensitivity study is conducted to investigate SO 2 and total ash mass evaluation errors due to uncertainty of the plume altitude input.Table 4 shows the plume altitude variation considered and the corresponding plume altitude (Z p ) and temperature (T p ) for both test cases.T p temperatures were obtained by interpolating the Trapani WMO atmospheric temperature profiles for the days in question.
Table 5 shows the SO 2 and total ash mass extrapolated using the novel procedure with different plume altitudes and temperatures.By considering an altitude uncertainty of ±500 m, the evaluation errors lie within 15 % for the 3 December 2006 event, and within 33 % for the 23 October 2011 test case.The evaluation error increases significantly when the plume altitude uncertainty is ±1000 m, becoming 30 % and 50 % for the 2006 and 2011 events, respectively.
As the mean altitude decreases (and temperature increases) the radiance emitted by the plume towards the sensor increases, and the layer of atmosphere above the plume is increasingly less negligible.For this reason (see Table 5) the uncertainty of SO 2 and ash abundances is greater when the plume altitude is lowered by an altitude variation, compared to when it is increased of the same amount.Almost the same relative variations with plume altitude were obtained using the LUT procedure.

Conclusions
A volcanic plume removal (VPR) procedure is described, involving computation of effective particle radius, aerosol optical depth of ash, and columnar content of both ash and SO 2 directly from remotely sensed TIR data.
The model applied is very simple: a plume at a constant altitude with a negligible layer of atmosphere above it.The procedure is reasonably independent of the surface underneath the plume as long as this is sufficiently uniform.Clearly, the layer of atmosphere above the plume is increasingly less negligible as plume altitudes reduce.Plume attenuation due to scattering mechanisms depends on particle size and concentration (aerosols optical depth).The proposed model accounts for these effects by using a modified plume temperature and two empirical values (one for dense plumes and one for transparent plumes) for vertical plume transmittance due to scattering.Finally, the results are refined by means of a polynomial relationship specific to the satellite sensor, the ash type, and the volcano.The polynomial coefficients are dependent on the sensor response functions, the optical properties of the ash, and the monthly averaged atmospheric profiles for the area of the volcano.The modified plume temperature and the polynomial coefficients were derived from a large number of MODTRAN simulations.
The polynomial parameters used in this study are valid for TERRA and AQUA MODIS data (bands 29, 31, and 32 at wavelength 8.6, 11.0, and 12.0 µm, respectively) collected over the Mt Etna volcano, and were tested for the usual tropospheric plume altitudes of the volcanic emissions.Moreover, they are applicable for the ash optical properties described by Volz (1973), which refer to an ash type considered typical of Mt Etna's emissions.Finally, the monthly mean profiles measured over the last 30 yr at Trapani (the WMO station closest to Mt Etna, and located in the western tip of Sicily) were used in the MODTRAN simulations.The parameters are stored in a secondary file and can easily be changed when recomputed for a different aerosol type (e.g.andesite) and/or sensor (e.g.SEVIRI).An important element in this methodology is the relationship between the ash transmittances of bands 31 and 29.This enables estimation of the SO 2 abundance from the total transmittance of band 29.The method is effective only if the foreseen type of ash is present.If more than one kind of particle is present in non-negligible percentages, the total effect of the aerosols is assumed to be due to the considered ash type alone, and the procedure gives erroneous results.
The comparison between the data used as input for the MODTRAN simulations and the results of the VPR procedure is quite positive (see Fig. 7).This agreement was in part expected because the procedure was adapted to these data, while the agreement between the VPR and LUT procedures was less predictable.Figure 9 very clearly shows the agreement between the two procedures, and the resulting maps exhibit the same structures.In particular, they both show the same plume characteristics in the maps of SO 2 and ash emitted on 23 October 2011.The SO 2 peak is shown at the beginning of the eruption (in the furthest downwind part of the plume), while the ash peak is temporally in the final portion of the phenomenon (in the portion of plume closest to the vents).
As demonstrated in the last part of the paper, the VPR procedure is quite sensitive to plume altitude (and temperature).A plume altitude estimate error of 1000 m (about 7 Kelvin degrees) can generate an SO 2 abundance percentage difference of 50 %.On the other hand, an altitude variation of 1000 m is not small for a plume at 5500 m and with a normal minimum altitude of 3500 m (top of Mt Etna).
It is worth noting that although the VPR procedure is described using specific case studies of Mt Etna eruptions, it can easily be adapted to any other volcano or ash type, provided new polynomial parameters are computed using local meteorological data sets and different ash optical properties.
After the initial derivation of the polynomial parameters for a specific volcano and meteorological conditions with a radiative transfer model, VPR requires no further atmospheric profiles (PTH) or subsequent radiative transfer model runs.Rapid plume information is retrieved using only the plume mean altitude and temperature as inputs.
It is simple to apply, applicable worldwide, and extremely fast, and could therefore prove very useful during a volcanic crisis, when rapid responses are particularly important.
Edited by: R. Schofield Fig. 1.(a) Radiance of the MODIS bands 28, 29, 31, and 32 in a transect perpendicular to the plume axis over the ocean.Circles: measured radiance.Thick lines: radiance without plume.Thin lines: straight-line tangent to absorption feature (plume) obtained considering only the points between the blue and cyan bars.(b) Radiance of the MODIS bands 28, 29, 31, and 32 in a transect perpendicular to the plume axis over both ocean (left) and land (right).Circles: measured radiance.Thick lines: radiance without plume.Thin lines: straight-line tangent to absorption feature (plume) obtained considering only the points between the blue and cyan bars.

Fig. 3 .
Fig. 3. From left to right: scatter plots of the MODTRAN simulated plume transmittances versus those computed using the procedure described, for channels 29, 31, and 32, respectively.The thick straight line represents a slope of unity.The r 2 correlation coefficients are, respectively, 0.9954, 0.9988, and 0.9982.

Fig. 4 .
Fig. 4. AOD of the MODIS TERRA band 31 vs AOD at 550 nm (input MODTRAN parameter).The red points contain all the MOD-TRAN simulations (12 monthly profiles, 4 plume altitudes, and 12 viewing angles).The coloured lines are the linear interpolation for the different effective ash particle radius (R e ).

Fig. 5 .
Fig. 5. Trend of TERRA m 31 m 32 and m 31 slopes versus the effective radius R e .

Fig. 6 .
Fig. 6.Plume transmittance of TERRA band 29 vs plume transmittance of TERRA band 31 considering only ash in the volcanic plume.Green line is a cubic polynomial fit.

Fig. 9 .
Fig. 9. SO 2 and ash assessments using the VPR (upper plates) and LUT (middle plates) procedures for the two Etna events considered.The lower plates show the flux comparison.

Table 2 .
Cubic polynomial relationship coefficients (see Eq. 15) computed for the MODIS-Terra and MODIS-Aqua channels.

Table 3 .
SO 2 and ash total masses and mean fluxes (v is the wind speed used) for the two MODIS test case images.All the results were obtained considering the plume only over the sea.

Table 4 .
Plume altitude variation, plume altitude, and plume temperature for the two test cases.

Table 5 .
SO 2 and ash total masses computed using the VPR procedure for different plume altitudes and temperatures.Total Mass [t] SO 2 Total Mass [t] Ash Total Mass[t]