Satellite-derived sulfur dioxide (SO₂) emissions from the 2014–2015 Holuhraun eruption (Iceland)

The 6-month-long 2014–2015 Holuhraun eruption was the largest in Iceland for 200 years, emitting huge quantities of sulfur dioxide (SO₂) into the troposphere, at times overwhelming European anthropogenic emissions. Weather, terrain and latitude made continuous ground-based or UV satellite sensor measurements challenging. Infrared Atmospheric Sounding Interferometer (IASI) data are used to derive the first time series of daily SO₂ mass present in the atmosphere and its vertical distribution over the entire eruption period. A new optimal estimation scheme is used to calculate daily SO₂ fluxes and average e-folding time every 12 h. For the 6 months studied, the SO₂ flux was observed to be up to 200 kt day⁻¹ and the minimum total SO₂ erupted mass was 4.4±0.8 Tg. The average SO₂ e-folding time was 2.4±0.6 days. Where comparisons are possible, these results broadly agree with ground-based near-source measurements, independent remote-sensing data and values obtained from model simulations from a previous paper. The results highlight the importance of using high-resolution time series data to accurately estimate volcanic SO₂ emissions. The SO₂ mass missed due to thermal contrast is estimated to be of the order of 3 % of the total emission when compared to measurements by the Ozone Monitoring Instrument. A statistical correction for cloud based on the AVHRR cloud-CCI data set suggested that the SO₂ mass missed due to cloud cover could be significant, up to a factor of 2 for the plume within the first kilometre from the vent. Applying this correction results in a total erupted mass of 6.7±0.4 Tg and little change in average e-folding time. The data set derived can be used for comparisons to other ground- and satellite-based measurements and to petrological estimates of the SO₂ flux. It could also be used to initialise climate model simulations, helping to better quantify the environmental and climatic impacts of future Icelandic fissure eruptions and simulations of past large-scale flood lava eruptions.

1 Introduction

Sulfur dioxide (SO₂) is one of the most important magmatic volatiles for volcanic geochemical analysis and hazard assessments due to its low ambient concentrations, abundance in volcanic plumes and spectroscopic features. Tropospheric volcanic SO₂ and its conversion products can affect the environment, human health, air quality and the radiative balance of the Earth (Schmidt et al., 2015; Gíslason et al., 2015; Schmidt et al., 2012; Gettelman et al., 2015; Ilyinskaya et al., 2017; Boichu et al., 2016; McCoy and Hartmann, 2015; Malavelle et al., 2017). Measurements of SO₂ from volcanic eruptions are vital, both to understanding the underlying volcanic processes and also the wider-scale environmental impacts of volcanism. The Icelandic Holuhraun eruption lasted from 31 August 2014 to 28 February 2015 and produced the largest lava volume in Iceland for more than 200 years (Gíslason et al., 2015). During September 2014, Holuhraun's average daily SO₂ emission exceeded daily SO₂ emissions from all anthropogenic sources in Europe by a factor of 3 (Schmidt et al., 2015 and references therein). The weather conditions and terrain made regular ground-based plume measurements extremely challenging during the winter months (Pfeffer et al., 2018). The high latitude of the eruption meant there was insufficient sunlight to reliably detect the volcanic plume using UV satellite sensors beyond the end of the October 2014. Due to insufficient sunlight, ground-based UV instruments did not measure SO₂ during the darkest 7 weeks of winter. Under these circumstances satellite-based thermal infrared spectrometers are an optimal source of high temporal resolution SO₂ amount and altitude.

The Infrared Atmospheric Sounding Interferometers (IASI) on board the Metop satellite platforms provide several observations of Holuhraun each day. The plume altitude and the SO₂ column amount are retrieved from the measured top-of-atmosphere spectral radiance (Carboni et al., 2012). For the first month of the Holuhraun eruption, previous studies have shown good agreement between IASI measurements and those from the Ozone Monitoring Instrument (OMI), ground-based and balloon-borne measurements and atmospheric dispersion model simulations (Schmidt et al., 2015; Vignelles et al., 2016). In this work IASI measurements are used to produce the first time series of the Holuhraun SO₂ plume. Retrievals of SO₂ amount from the Metop-A satellite are binned and averaged for successive 12 h periods to give global coverage twice a day for the entire period of the eruption. The time series of the SO₂ mass present in the atmosphere is used to calculate SO₂ fluxes and an average SO₂ e-folding time, under the assumption that the flux is constant over a 12 h period. The results are compared with ground-based Brewer measurements of the SO₂ column amount and with measurements of near-source plume altitudes and fluxes from the Icelandic Meteorological-Office (IMO). The data set presented can be used for comparisons to other ground- and satellite-based measurements and to petrological estimates of the SO₂ flux and to initialise, for instance, climate model simulations, helping to better quantify the environmental and climatic impacts of volcanic SO₂.

2 IASI SO₂ iterative retrieval scheme

IASI is an infrared Fourier transformer interferometer on board the Metop-A and Metop-B satellites. It measures in the spectral range 645–2760 cm⁻¹ with spectral sampling of 0.25 cm⁻¹ and has global coverage every 12 h. The IASI data set used in this study was the level 1c (apodised) data set from the EUMETSAT and CEDA archive.

The details of the retrieval scheme are summarised briefly below. For more details see Carboni et al. (2012, 2016). An SO₂ retrieval is performed for all IASI pixels that present a positive result in the SO₂ detection scheme (Walker et al., 2011, 2012). The detection scheme uses all the channels in the ν₃ band (1300–1410 cm⁻¹), while the iterative retrieval uses all the channels in both spectral ranges, ν₁ and ν₃ (1000–1200 and 1300–1410 cm⁻¹). The strongest SO₂ band is the ν₃, around 7.3 µm and is contained within a strong water vapour (H₂O) absorption band so it is not very sensitive to emission from the surface and lower atmosphere. However, above the lower atmosphere, this band contains valuable information on the vertical profile of SO₂. Fortunately, differences between the H₂O and SO₂ absorption spectra allow the signals from the two gases to be decoupled in high resolution measurements. The SO₂ ν₁ band around 8.7 µm (1000 to 1200 cm⁻¹) is within an atmospheric window. This allows the radiation from the surface to reach the satellite from deep within the atmosphere, enabling the retrieval of SO₂ down to the surface.

The detection scheme is a linear retrieval with one free parameter, the column amount of SO₂. In particular we assume the vertical distribution of SO₂ and the atmospheric vertical profiles (temperature and trace gases). We do not take into account negative thermal contrast so that regions with negative thermal contrast give (i.e. where the surface is colder than the atmosphere) zero (or negative) values of the SO₂ column amount. In our scheme we consider detection to be positive if the output of the linear retrieval is greater than a defined positive threshold (0.49 effective DU, following Walker et al., 2012). The detection limits for a standard atmosphere (with no thermal contrast) are estimated to be 17 DU for a SO₂ plume between 0 and 2 km, 3 DU between 2 and 4 km and 1.3 DU between 4 and 6 km (Walker et al., 2011). The detection scheme can miss part of an SO₂ plume under certain circumstances, such as low-altitude plumes, conditions of negative thermal contrast and where clouds are present above the SO₂ plume, masking the signal from the underlying atmosphere. Due to these uncertainties the estimated mass of SO₂ in this paper should be regarded as a minimum.

Figure 1 Map showing the maximum of the SO₂ column amount (in Dobson units, DU) retrieved within the considered area from 30 to 90^∘ N (black rectangle) from September 2014 to February 2015.

We perform the iterative retrieval for the pixels that give positive detection results. All the channels in the ranges 1000–1200 and 1300–1410 cm⁻¹ (the 7.3 and 8.7 µm SO₂ bands) are simultaneously used in the iterative optimal estimation retrieval scheme to obtain the SO₂ amount, the altitude of the plume and the surface temperature. The scheme iteratively fits the forward model (simulations) with the measurements, through the error covariance matrix, to seek a minimum of a cost function. The forward model (FM) is based on RTTOV (Saunders et al., 1999), extended to include SO₂ explicitly, and uses ECMWF profiles interpolated to the measurement time and location. The error covariance matrix used is the global error covariance matrix in Carboni et al. (2012). It is defined to represent the effects of atmospheric variability not represented in the FM, as well as instrument noise. This includes the effects of cloud and trace gases which are not explicitly modelled. The matrix is constructed from differences between FM calculations (for clear-sky driven with ECMWF profiles) and actual IASI observations for a wide range of conditions, when we are confident that negligible amounts of SO₂ are present. Only quality-controlled pixels are considered: these are values where the minimisation routine converges within 10 iterations, the SO₂ amount is positive, the plume pressure is between 0 and 1100 mb and the cost function is less than 10. A comprehensive error budget for every pixel is included in the retrieval. Rigorous error propagation, including the incorporation of FM parameter error, is built into the system, providing quality control and comprehensive error estimates on the retrieval results. The IASI SO₂ retrieval is not affected by underlying cloud. If the SO₂ is within or below a cloud layer its signal will be masked and the retrieval will underestimate the SO₂ amount.

3 Temporal evolution of SO₂ mass and SO₂ vertical distribution

The retrieval results from the Metop-A orbits during the period from September 2014 to February 2015 and from 30 to 90^∘ N are combined (twice a day, i.e. morning and afternoon overpasses) to produce maps of retrieved SO₂ amount and altitude. The animation in the Supplement (S1) shows the evolution of the plume for each day. The Holuhraun eruption is the main source of SO₂ over that period. Other minor sources include SO₂ emitted from intermittent volcanic activity on the Kamchatka Peninsula, Etna volcano (28 December 2014, 1, 2 and 21 January 2015) and anthropogenic SO₂ emissions from Beijing, China. Satellite observations at the pixel level do not provide sufficient information to distinguish between SO₂ from Holuhraun and SO₂ from other sources. For example, the elevated SO₂ near Beijing on 21 December 2014 appears to be from an anthropogenic source, but the elevated SO₂ in the same area on 31 December 2014 is from the Holuhraun eruption. In this study all the SO₂ measured from 30 to 90^∘ N between September 2014 and February 2015 is referred to as Holuhraun SO₂.

Over the course of 6 months the eruption plume dispersed across the Northern Hemisphere. Figure 1 shows the maximum SO₂ column amount retrieved during the 6-month period and illustrates that SO₂ from the Holuhraun eruption was dispersed over large parts of the Northern Hemisphere, including poleward of the Arctic circle. For the majority of the time, the plume circulated around the pole and the northern regions (see animation in S1), overpassing Scandinavia, eastern Europe, Russia, Greenland and Canada several times. The plume overpassed Europe on multiple occasions, most often northern Europe (Schmidt et al., 2015; Ialongo et al., 2015; Zerefos et al., 2017; Steensen et al., 2016; Twigg et al., 2016; Boichu et al., 2016), but also Italy (22 October 2014) and Spain, and reached as far south as Morocco, Algeria (on 5 and 6 November 2014 respectively) and the Greece–Macedonia–Albania region (5 and 6 January 2015).

Figure 2 Time series of SO₂ mass as a function of day from 1 September 2014. The plot shows the IASI mass (obtained from the retrieval) with error bars in black and the fitting of the flux retrieval (Sect. 5) with the cyan line. Red circles and orange line show the masses corrected for cloud presence and their fitting of the retrieval.

IASI orbits are grouped into 12 h intervals in order to have two maps each day of IASI-retrieved SO₂ amount and altitude. These maps are gridded onto 0.125^∘ latitude and longitude boxes following Carboni et al. (2016). The SO₂ mass time series is obtained by summing the mass values of the regularly gridded map for each 12 h period. The same procedure has been used for the errors, this means that the sums of errors of grid boxes are considered as errors on the total masses (this could be an error overestimation but we cannot consider the usual errors in quadrature due to the possible presence of systematic error, e.g. the errors are not independent). The time series of SO₂ mass, together with the errors, are presented in Fig. 2. The largest SO₂ mass is found in September 2014 (up to 0.25 Tg), when the eruption was most powerful. The SO₂ mass decreases during October 2014 (with some peak values around 0.1 Tg), then increases around end of November and beginning of December 2014 (up to 0.15 Tg). The SO₂ mass steadily declines during January and February 2015 as the eruption comes to an end. There is no detection of a SO₂ plume attached to the vent in the second half of February (and the SO₂ mass for this period, reaching a value up to 0.01 Tg, is from a non-Icelandic source).

The SO₂ mass present between two altitude levels was estimated using the method by Carboni et al. (2016) to produce the vertical distribution of SO₂. In this study the vertical distribution of SO₂ was estimated every 12 h from 0 to 10 km with a vertical resolution of 0.5 km for all latitudes north of 30^∘ N. Both the young emitted plume, as well as the mature plume that had been transported around in the Northern Hemisphere for a few days, are included in the distribution. The time series of the SO₂ vertical distribution for the Holuhraun eruption is shown in Fig. 3. The centre of mass of the plume closest to the vent can be used as a rough estimate for the injection height. Figure 3 shows the time series of two data sets: (i) the vertical distribution and (ii) the altitude of the centre of mass of the SO₂ values within 500 km of the vent. The altitude of the centre of mass is less than 4 km for the majority (96 %) of the measurements.

Figure 3 SO₂ vertical distribution in kilometres above sea level. The colours represents the mass of SO₂ and dark red represents values greater than 0.03 Tg. Every column of the plot is generated from an IASI map (one every 12 h; Fig. 1 and Supplement files). The black diamonds show the altitude of the centre of mass (CM) computed with the IASI pixels within 500 km from the vent, and the red and black dots show the DOAS measurements and other observations from Pfeffer et al. (2018).

4 Quantifying satellite-retrieval underestimation

Part of the SO₂ plume can be missed by this IASI scheme, and the derived SO₂ masses should be considered as a minimum. The presence of low thermal contrast can prevent the detection of SO₂. While cloud below the plume should not be a problem for this scheme, the presence of cloud above the plume will smooth the spectral signature of SO₂, causing underestimation of the retrieved SO₂ amount (Carboni et al., 2012). At its most extreme this effect can completely erase the SO₂ signal (for a cloud layer with zero transmittance).

We estimate the percent of SO₂ mass missing due to cloud above the plume as a function of cloud optical depth and altitude above the SO₂ plume using the same simulations with a default atmosphere used in Fig. 6 of Carboni et al. (2012). Using the European Space Agency cloud Climate Change Initiative data set (Stengel et al., 2017, ESA cloud CCI) of Advanced Very High Resolution Radiometer – AVHRR (same platform as IASI, i.e. same local overpassing time) L3 monthly mean statistic, we computed the following:

• Monthly mean histograms of frequency of cloud optical depth at 550 nm, τ, was averaged over the globe. τ is not present in the cloud L3 database for locations without daylight (e.g. visible channels) and most of the Icelandic plume in the winter months is without daylight. As a consequence, here we are assuming the global histogram of frequency of τ is valid over the plume region.

• Monthly mean histograms of frequency of cloud altitude are averaged over the plume region (30 to 90^∘ N). Cloud altitude is available for all locations and during winter months.

We consider the measured mass M_meas to be the difference between the true mass M_corr and the missing one M_miss:

$$\begin{array}{}\text{(1)}& {\displaystyle }{M}_{\mathrm{corr}}-M_{\mathrm{miss}}=M_{\mathrm{meas}}.\end{array}$$

We can rewrite this as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{corr}}\left(\mathrm{1}-{\displaystyle \frac{M_{\mathrm{miss}}}{M_{\mathrm{corr}}}}\right)=M_{\mathrm{meas}}\text{(3)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{corr}}=M_{\mathrm{meas}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\frac{M_{\mathrm{miss}}}{M_{\mathrm{corr}}}}}\right).\end{array}$$

The term within the bracket on the right is the correction factor. We compute the correction factor, C, for every month of the eruption as a function of altitude and apply it to the vertical distribution data set.

$$\begin{array}{}\text{(4)}& {\displaystyle }{M}_{\mathrm{corr}}\left(h\right)=M_{\mathrm{meas}}\left(h\right)\cdot C\left(h\right),\end{array}$$

with

$$\begin{array}{}\text{(5)}& {\displaystyle }C\left(h\right)={\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-Z(h\left)\right)}},\end{array}$$

where Z(h) is the fraction of SO₂ “missed” in the measurements due to cloud above the plume. Z(h) is estimated as the product of the probability of having cloud above altitude h, F(h) times the attenuation due to cloud, A.

$$\begin{array}{}\text{(6)}& {\displaystyle }Z\left(h\right)=F\left(h\right)\cdot A\end{array}$$

The probability of having cloud above h has been estimated from the ESA cloud CCI data for the region considered for the volcanic plume (latitude >30^∘ N) as the number of cloud retrievals above altitude h divided by number of observations.

Attenuation due to cloud (A) is the sum of the frequency of having a cloud with a cloud optical depth f(τ) times the attenuation due to a cloud with the same optical depth a(τ).

$$\begin{array}{}\text{(7)}& {\displaystyle }A=\sum _{\mathit{\tau }=\mathrm{0}}^{n}f\left(\mathit{\tau }\right)a\left(\mathit{\tau }\right)\end{array}$$

f(τ) has been estimated using the monthly mean histogram of frequency of cloud optical depth, estimated over the globe. a(τ) has been estimated by running the SO₂ retrieval using, as IASI measurements, simulated spectra with water cloud above the plume, using the default atmosphere and different optical depths τ. For optical depths bigger than 10 the attenuation is 1 (cloud is opaque and completely masks the SO₂ signal). Figure 4 shows the correction factor together with the SO₂ vertical distribution obtained from the IASI retrieval and SO₂ vertical distribution corrected for underestimation due to cloud cover.

Figure 4 SO₂ vertical distribution in kilometres above sea level. The colours represent the mass of SO₂, and dark red represents values higher than the colour bar. Every column of the plot is generated from an IASI map (one every 12 h). Panel (a) shows the data obtained from the IASI maps. Panel (b) is (a) times the correction factor (to include SO₂ that statistically has been missed by the IASI measurements due to cloud above the plume). Panel (c) is the correction factor C(h).

In the following section the emission fluxes have been estimated with both the original SO₂ masses (obtained from IASI retrieval scheme explained in Sect. 2) and the masses corrected by cloud cover (this section).

We estimate the missing SO₂ mass due to thermal contrast by comparison with OMI SO₂ (Theys et al., 2015) for the months of September and October 2014, as the OMI data set does not cover all of the time period of the eruption due to the lack of solar radiance during the winter. We compare, with visual inspection, the daily maps of IASI and OMI and identify the parts of the plume missing from IASI detection (and consequently missing in the IASI retrieval). Table 1 presents the list of all areas identified and the SO₂ estimate from OMI (BIRA-IASB algorithm).

Table 1 Areas of plume missing from IASI detection scheme and OMI SO₂ mass estimate in the area, assuming two different altitudes.

The total mass of SO₂ missed due to thermal contrast is estimated to be a few percent of the emission estimate by IASI. In particular, the missing plume for the first 2 months, estimated with OMI 0–1 km, has a total mass of 0.08 Tg of SO₂ that corresponds to 3 % percent of the emission estimate by IASI (see Sect. 5 for emission estimate). Cloud could also be a problem for OMI, for such low-altitude SO₂, and OMI values could be off by a factor of 2–3, which will change the underestimation to 6 %–9 % (instead of 3 %).

5 Daily SO₂ fluxes

The time series of SO₂ fluxes and the coefficient of an average exponential decay of SO₂ (the state vector x) are retrieved from the time series of SO₂ mass (the measurement vector y) using the optimal estimation scheme of Rodgers (2000). The solution (x, the vector of parameters that we want to retrieve) gives the most probable x given the measurements y and the a priori knowledge x_a. The state vector x is found by minimising the cost function:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathit{\chi }}^{\mathrm{2}}=[\mathit{y}-\mathbf{F}(\mathit{x}\left){]}^T{\mathbf{S}}_{\mathrm{e}}^{-\mathrm{1}}\right[\mathit{y}-\mathbf{F}\left(\mathit{x}\right)]+[\mathit{x}-{\mathit{x}}_{\mathrm{a}}{]}^T{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}[\mathit{x}-{\mathit{x}}_{\mathrm{a}}],\end{array}$$

where F is the forward model (that simulates the measurement given the state vector), x_a is the a priori value of the state vector and S_e and S_a are the measurement and a priori error covariance matrices. S_e and S_a are diagonal matrices with the variances (square of errors) of y and x_a respectively as their diagonal elements. The a priori values used were 0.2±0.2 Tg day⁻¹ for flux and 2±2 day for the e-folding time: these values reflect the variability of lifetime found in the literature for low tropospheric emission (Carn et al., 2016; Lee et al., 2011).

To retrieve the time series of SO₂ fluxes the state vector x was defined as follows. The first element of the state vector is the average SO₂ e-folding time (λ) for the period analysed; the following x_i elements are the average emission flux f_i between the two IASI estimates of SO₂ mass at time t_i−1 and t_i:

$$\begin{array}{}\text{(9)}& {\displaystyle }\mathit{x}=(\mathit{\lambda },f_{\mathrm{1}},f_{\mathrm{2}},\mathrm{\dots },f_n).\end{array}$$

To define the forward model F we consider that the SO₂ mass m measured by the satellite is a function of time and that the mass decays proportionally to the mass itself, plus a source term of flux f is added. These terms give a first-order differential equation:

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\mathrm{d}m}{\mathrm{d}t}}=-km+f.\end{array}$$

Assuming a constant flux f over the time interval Δt between two consecutive mass estimates m_i and m_i−1, the solution becomes

$$\begin{array}{}\text{(11)}& {\displaystyle }{m}_i=m_{i-\mathrm{1}}{e}^{-\frac{\mathrm{1}}{\mathit{\lambda }}\mathrm{\Delta }t}+f\mathit{\lambda }(\mathrm{1}-e^{-\frac{\mathrm{1}}{\mathit{\lambda }}\mathrm{\Delta }t}),\end{array}$$

where λ (with $\mathit{\lambda }=\mathrm{1}/k$ ) is the average e-folding time. Equation (11) is the forward model F(x).

Figure 5 Time series of emission fluxes of SO₂. To facilitate the comparison with the flux estimates from Schmidt et al. (2015), (a) presents values (for the first month only) of the retrieved SO₂ flux time series from IASI measured mass in black with grey error bars, the retrieved flux time series from corrected IASI mass in red with orange error bars and DOAS ground measurements from Pfeffer et al. (2018) in blue (blue bars show the errors, dotted bars show the maximum and minimum values measured that day). Panel (b) shows the 6-month time series with the same colour coding.

We used optimal estimation (OE) to fit the measured time series of mass loading with a forward model that reproduces the time series as a function of emission fluxes and SO₂ e-folding time. This forward model (Eq. 11) is the inverse of Eq. (6) in Theys et al. (2013), section “Delta-M method”. Both equations derive from the solution of the same differential equation (Eq. 10). In Theys et al. (2013) the fluxes are obtained assuming the e-folding time; here, the fluxes and the averaged e-folding time and their uncertainties are estimated simultaneously based on OE.

Figure 2 shows the fitting of the measurements y with the forward model F, while Fig. 5 shows the retrieved fluxes with errors.

The flux time series follows a similar pattern to the SO₂ mass time series, with the highest values in September 2014. IASI determined fluxes are provided here with the ground-based measurement fluxes provided in parenthesis (described in Sect. 5). The September daily average was 0.06 Tg day⁻¹ (0.09) and the September daily maximum was 0.26 Tg day⁻¹ on the 20 September (0.19). The October average flux was 0.023 Tg day⁻¹ (0.083), with peak values of 0.15 and 0.12 Tg day⁻¹ on the 11 and 19 October (0.10). The November average flux was 0.029 Tg day⁻¹ (0.069). The flux increases in the second part of November to a maximum of 0.13 Tg day⁻¹ on the 23 November (0.13). The estimates for December, January and February show decreasing fluxes with monthly averages of 0.016, 0.006 and 0.005 Tg day⁻¹ respectively (0.026, 0.028, 0.016). The monthly averages are lower than those measured by the ground-based measurements, while the maximum daily averages for each month are generally higher. The UV ground-based measurements for the dark months of December, January and February are sparse, with only 10 measurements over these 3 months. There was only 1 day with measurements at the beginning of December, then 6 days with measurements in the second half of January and 3 days with measurements in February. The extrapolated flux from the ground-based measurements through December to the first half of January is consistent with the error bars from the IASI estimates. The differences in monthly means between the ground-based measurements and the IASI flux estimates in the sunny months are explainable by low values with large error. The fluxes calculated for September 2014 are consistent with Schmidt et al. (2015) (e.g. up to 0.12 Tg day⁻¹ during early September, 0.02–0.6 Tg day⁻¹ between 6 and 22 September, 0.06–0.12 Tg day⁻¹ until the end of September).

The total mass of SO₂ emitted by the eruption is obtained by summing all the fluxes, f_i, output by the retrieval and multiplying by the corresponding Δt_i. The error associated with the total mass emitted is obtained by adding in quadrature the errors δf_i multiplied by the time interval Δt_i. The maximum value of total mass emitted is obtained by summing all the fluxes plus the uncertainty and the minimum value is obtained by summing all the fluxes less the uncertainty (negative values are set to zero). This gives a total mass of emitted SO₂ of 4.4±0.8 Tg with a maximum of 11.6 Tg and a minimum of 0.4 Tg. The retrieved averaged e-folding time is 2.4±0.6 days.

Note that, in using an average SO₂ e-folding time for the entire eruption period, any variation in e-folding time will be interpreted as an inverse variation in the estimated flux; i.e when the “real” e-folding time is higher than the retrieved one, the flux will be overestimated and vice versa. The flux uncertainties include the errors in flux due to the variation in e-folding time. The SO₂ lifetime can vary significantly on a daily basis mainly as a function of water vapour and solar irradiation. Also, note that any e-folding time shorter than the retrieved one can fit the measurements and give higher fluxes. Given these caveats, the value of the e-folding time (2.4±0.6 days) is consistent with the mean lifetime of 2.0 days estimated for September 2014 (Schmidt et al., 2015). Their estimate was based on minimising the difference between the SO₂ amount from the NAME dispersion model and the IASI and OMI satellite measurements.

Figure 5 shows that higher values (peaks) of SO₂ flux often alternate with lower values (below 0.02 Tg day⁻¹), even during periods that have been identified as generally characterised by higher fluxes. This intermittent flux behaviour has important implications in terms of the estimate of total SO₂ emitted.

Figure 6 Time series of the SO₂ column amount (DU) measurements at Sodankylä as a function of day from the 1 September 2014. Black symbols are the Brewer measurements with error bars, red symbols are the mean and standard deviation of all the IASI measurements within 200 km, and blue lines represent the maximum and minimum of the IASI measurements.

Gauthier et al. (2016) used thermal infrared (TIR) data from SEVIRI on board the geostationary satellite Meteosat Second Generation (MSG) to retrieve an SO₂ mass time series from 1 September 2014 to 25 November 2014. Their daily retrieved mass values are lower compared to the IASI values reported here due to the smaller geographic area considered and possibly due to the lower sensitivity of SEVIRI. Nevertheless they estimate a total SO₂ emitted mass of 8.9±0.3 Tg for the period September 2014 to November 2014, which is a factor of 2 higher than IASI.

Our understanding is that the Gauthier et al. (2016) scheme produced valid measurements when the SO₂ loading increased with respect to the previous measurement. This leaves data gaps when the SO₂ measured remains constant or decreases. The resulting data gaps were filled by linearly interpolating between two valid measurements. This has a tendency to bias flux estimates in favour of increasing SO₂ loading. The data set of Gauthier et al. (2016) contains several days with no valid flux estimates and, for these data gaps, the interpolation of valid measurements into the gaps could account for their discrepancies with our data set. The ground-based measurements reported in Pfeffer et al. (2018) suggest an intermediate value of 7.3 Tg over this time interval.

6 Comparison with ground-based and near-source measurements

The conditions of the Holuhraun eruption are significantly different to eruptions investigated in previous studies using IASI (Carboni et al, 2016), because the plume from Holuhraun was confined to altitudes between the surface and 6 km at high-northern latitudes and because the eruption took place during the autumn and winter months. As a result there is less radiance and low (or negative) thermal contrast between the surface and the first layer of the atmosphere. These conditions lead to lower sensitivity for measurements in both the UV and TIR spectral range. Nevertheless, a cross-comparison with other available measurements is informative when assessing our results. The following comparisons are an addition to previous comparisons done between UV satellite measurements and a dispersion model (Schmidt et al., 2015) and balloon measurements (Vignelles et al., 2016).

First the IASI data set is compared with ground-based Brewer measurements of the SO₂ column amount of the mature plume over Finland. Successively, the plume height is compared with near-source measurements in Iceland using ground- and aircraft-based visual observations, web camera and NicAIR II infrared images, triangulation of scanning DOAS instruments and the location of SO₂ peaks measured by DOAS traverses as reported in Pfeffer et al. (2018).

The Brewer ground measurements (Ialongo et al., 2015) were made at Sodankylä (67.42^∘ N, 26.59^∘ E). The SO₂ column amounts are routinely obtained from the direct solar irradiances at wavelengths of 306.3, 316.8 and 320.1 nm by using the same Brewer algorithm as for the ozone retrieval (Kerr et al., 1988). The method is based on the Lambert–Beer law, which describes the attenuation of the direct solar irradiance reaching the Earth's surface at certain wavelengths due to the atmospheric constituents. In order to avoid the effects of stray light at short wavelengths, the measurements corresponding to large air mass values (after 14:20 UT) are not included. The SO₂ column amounts in Sodankylä are typically close to zero with an estimated detection limit of about 1 DU. Ialongo et al. (2015) compared the SO₂ column amount values from the Brewer and OMI measurements in Sodankylä during September 2014. The differences between OMI and Brewer retrievals were usually smaller than 2 DU.

The comparison here is performed by averaging all the IASI pixels that pass quality control within 200 km of the ground measurements. The Brewer instrument measures in the UV and thus only in daylight conditions. This means that only the first month of the eruption can be considered. Figure 6 shows the time series of the SO₂ column amount from both ground measurements and IASI as a function of time. All the plume episodes (with SO₂ amount larger than 2 DU) are consistent between the two data sets with the exception of 15 and 19 September when the plume only passes over the northern part of the 200 km circle in the IASI data and does not pass over the ground measurement station. This means that, while IASI presents high SO₂ measurements, elevated values are not observed in the Brewer measurements.

There are a few days of low (less the 2 DU) SO₂ reported by only one of the two instruments (IASI or Brewer), meaning that the detection limit of one of the instruments is not exceeded. This is consistent with the IASI minimum error for low-altitude plumes (2 DU for plumes below 2 km,  Carboni et al., 2016) and with the Brewer minimum error (2 DU,  Sellitto et al., 2017).

The ground-based measurements of the eruption cloud-top height were collected from multiple techniques including ground- and aircraft-based observations, web camera, ScanDOAS, MobileDOAS and NICAIR II IR camera (Pfeffer et al., 2018). The ScanDOAS and MobileDOAS approaches and IASI retrievals, provide the height of the centre of mass of the plume, while the other techniques provide the height of the top of the plume. Figure 3 presents the ground-based and IASI altitudes together. In general, the IASI and ground-based altitudes agree that the altitudes varied mainly between 1 and 3 km; however they do not agree particularly well on any specific day.

The time series of the ground-based (Pfeffer et al., 2018) and IASI flux estimates are presented in Fig. 5. Within error, they generally agree. There are a few significant differences between the two data sets: two values in October–November 2014 and some values at the end of January and February. The ground-based measurements in November alternate between very high and very low values. On 5 November 2014, the ground-based value is significantly higher than the IASI flux estimate for that day. Pfeffer et al. (2018) suggest the high values in November could be due to degassing from a continually replenishing lava lake contributing to the total gas in addition to the degassing from the magma being erupted. The plume altitude was less than 2 km on this day and under these conditions IASI values can be underestimates. The total mass emitted can be estimated using the ground-based measurements as the integral below the red line in Fig. 2 (i.e. interpolating flux values). Even if the majority of fluxes are consistent with each other within the error estimate, the total mass calculated by IASI (4.4 Tg) and IMO (9.6 Tg) differs by a factor of 2. The 5 November discrepancy contributes significantly to this difference.

7 Conclusions

The first satellite-based SO₂ flux data set of the full 2014–2015 Holuhraun eruption has been estimated using the IASI instruments. The data set provides a flux estimate every 12 h for the entire eruption. Thermal infrared spectrometers are the only satellite instruments that could follow the SO₂ plume around the Arctic in the absence of solar irradiation during the winter months of the eruption. The low altitude of the SO₂ plume and cold underlying surface reduces IASI sensitivity to SO₂; however the results compare reasonably well with ground-based near and distal measurements. The observations show that the Holuhraun plume passed over large parts of the Northern Hemisphere during the eruption. The time series of the vertical distribution of SO₂ showed a low-altitude plume confined mainly within 0–6 km. The time series of SO₂ masses showed a maximum of 0.25 Tg of atmospheric loading in September 2014. A new optimal estimation scheme was developed to calculate daily SO₂ fluxes and e-folding time based on satellite-retrieved atmospheric SO₂ burdens. Application of the method gave estimates of SO₂ flux of up to 200 kt day⁻¹. The minimum total mass of SO₂ was calculated to be 4.4±0.8 Tg and the average SO₂ e-folding time was found to be 2.4±0.6 days. By comparison with the OMI data set during the summer months of the eruption we estimate that the SO₂ masses missed due to low thermal contrast are of the order of a few percent (3 %) of the total emission. We estimated the SO₂ mass missed due to cloud above the SO₂ plume that is masking the signal, using AVHRR cloud CCI data set monthly mean.

Results show a correction factor increasing with decreasing altitude from one (no underestimation of SO₂ masses) at the top of the cloud around 9 km up to a factor of 2 (we measure 50 % of the true mass) for plumes between 0 and 1 km. Applying this correction resulted in a total mass, emitted during the 6 months of the eruption, of 6.7±0.4 Tg and little change in the average e-folding time (2.5±0.7). The IASI flux data reported here are representative of every 12 h and with no data gaps. Other emission estimates include the interpolation (or extrapolation) of fluxes where there are no measured data; this could produce discrepancies between estimates made using different methods.