A satellite chronology of plumes from the April 2021 eruption of La Soufrière, St Vincent

Satellite instruments play a valuable role in detecting, monitoring and characterising emissions of ash and gas into the atmosphere during volcanic eruptions. This study uses two satellite instruments, the Infrared Atmospheric Sounding Interferometer (IASI) and the Advanced Baseline Imager (ABI), to examine the plumes of ash and sulfur dioxide (SO₂) from the April 2021 eruption of La Soufrière, St Vincent. The frequent ABI data have been used to construct a 14 d chronology of a series of explosive events at La Soufrière, which is then complemented by measurements of SO₂ from IASI, which is able to track the plume as it is transported around the globe. A minimum of 35 eruptive events were identified using true, false and brightness temperature difference maps produced with the ABI data. The high temporal resolution images were used to identify the approximate start and end times, as well as the duration and characteristics of each event. From this analysis, four distinct phases within the 14 d eruption have been defined, each consisting of multiple explosive events with similar characteristics: (1) an initial explosive event, (2) a sustained event lasting over 9 h, (3) a pulsatory phase with 25 explosive events in a 65.3 h period and (4) a waning sequence of explosive events. It is likely that the multiple explosive events during the April 2021 eruption contributed to the highly complex plume structure that can be seen in the IASI measurements of the SO₂ column amounts and heights. The bulk of the SO₂ from the first three phases of the eruption was transported eastwards, which based on the wind direction at the volcano implies that the SO₂ was largely in the upper troposphere. Some of the SO₂ was carried to the south and west of the volcano, suggesting a smaller emission of the gas into the stratosphere, there being a shift in wind direction around the height of the tropopause. The retrieved SO₂ heights show that the plume had multiple layers but was largely concentrated between 13 and 19 km, with the majority of the SO₂ being located in the upper troposphere and around the height of the tropopause, with some emission into the stratosphere. An average e-folding time of 6.07±4.74 d was computed based on the IASI SO₂ results: similar to other tropical eruptions of this magnitude and height. The SO₂ was trackable for several weeks after the eruption and is shown to have circulated the globe, with parts of it reaching as far as 45^∘ S and 45^∘ N. Using the IASI SO₂ measurements, a time series of the total SO₂ mass loading was produced, with this peaking on 13 April (descending orbits) at 0.31±0.09 Tg. Converting these mass values to a temporally varying SO₂ flux demonstrated that the greatest emission occurred on 10 April with that measurement incorporating SO₂ from the second phase of the eruption (sustained emission) and the beginning of the pulsatory phase. The SO₂ flux is then shown to fall during the later stages of the eruption: suggesting a reduction in eruptive energy, something also reflected in ash height estimates obtained with the ABI instrument. A total SO₂ emission of 0.63±0.5 Tg of SO₂ has been derived, although due to limitations associated with the retrieval, particularly in the first few days after the eruption began, this, the retrieved column amounts and the total SO₂ mass on each day should be considered minimum estimates. There are a number of similarities between the 1979 and 2021 eruptions at La Soufrière, with both eruptions consisting of a series of explosive events with varied heights and including some emission into the stratosphere. These similarities highlight the importance of in-depth investigations into eruptions and the valuable contribution of satellite data for this purpose; as these studies aid in learning about a volcano's behaviour, which may allow for better preparation for future eruptive activity.

1 Introduction

La Soufrière (13.33^∘ N, 61.18^∘ W, summit elevation of 1220 m), a volcano on the island of St Vincent, entered a phase of explosive activity on 9 April 2021 after having been in a lower level state of eruption, including the slow extrusion of a lava dome, since late December 2020 (Global Volcanism Program, 2021a; Joseph et al., 2022). The volcano is part of the eastern Caribbean volcanic arc and has erupted on at least five occasions since the 18th century in both explosive (1718, 1812, 1902–1903, 1979) and lava-dome-forming eruptions (1971–1972, 1979) (Robertson, 1995; Pyle, 2017).

Eruptive activity since 1970 included the non-explosive extrusion of a basaltic andesite lava dome into the flooded crater from November 1971–January 1972 and a violent series of explosions that began in the flooded summit crater on 13 April 1979 (Aspinall et al., 1973; Shepherd et al., 1979). The 1979 explosions were very well documented at the time and persisted for 2 weeks. They were followed by about 6 months of lava dome extrusion across the crater floor, which by then had been infilled by pyroclastic ejecta (Brazier et al., 1982; Fiske and Sigurdsson, 1982; Shepherd and Sigurdsson, 1982). The explosive eruption of La Soufrière in April 1979 was one of the first eruptions to have occurred during the “satellite era” (Fiske and Sigurdsson, 1982). Observations from the infra-red radiometer on the Synchronous Meteorological Satellite 1 (SMS-1) were used to measure the growth of the volcanic ash plumes (Krueger, 1982), the Stratospheric Aerosol and Gas Experiment (SAGE) measured stratospheric aerosol from the eruption (McCormick et al., 1982), and SO₂ from the eruption was observed by the Total Ozone Mapping Spectrometer (TOMS; Carn et al., 2003, 2016).

After the 1979 eruption, La Soufrière showed no detectable signs of activity or unrest until late 2020. An increase in seismicity was noted by seismologists at The University of the West Indies Seismic Research Centre (UWI-SRC) in November and December 2020 (Joseph et al., 2022). An effusive eruption began on 27 December 2020 with the emplacement of a new lava dome, which grew over the following months. On 8 April 2021, the alert level was raised to the highest level and an evacuation of the highest risk communities was ordered (Global Volcanism Program, 2021b). The first of a number of explosive eruptive events began on 9 April 2021 at 08:41 LT (12:41 UTC) (Joseph et al., 2022). Multiple explosive events occurred over the following 2 weeks, with the last explosive event occurring on 22 April 2021. Activity during this eruptive period led to the closure of local airports, and ashfall affected much of St Vincent, as well as neighbouring islands, including the Grenadines, Barbados and Saint Lucia (Global Volcanism Program, 2021b). Details of the eruption, its impacts and the crisis management can be found within a special issue on the eruption published by the Geological Society (Robertson et al., 2024).

Technological developments in the 4 decades since the 1979 eruption mean that numerous aspects of the 2021 eruption were observed with multiple satellite instruments, including thermal anomalies; dome growth; lightning; and the evolution of the SO₂, ash, and sulfate plumes (Global Volcanism Program, 2021a, b; Smart and Sales, 2021; Babu et al., 2022; Thompson et al., 2022; Yue et al., 2022; Horváth et al., 2022; Dualeh et al., 2023; Bruckert et al., 2023; Camejo-Harry et al., 2024; Esse et al., 2024). Measurements made by satellite instruments allow the detection of SO₂ and ash and quantification of the plume properties, which is essential for assessing the potential hazard to aviation (Prata and Tupper, 2009; Thomas and Watson, 2010; Lechner et al., 2017) and providing more sophisticated information on eruption source parameters (Aubry et al., 2021). Of particular value is the ability of satellite instruments to track the evolution of these plumes as they are transported away from the source. Geostationary instruments with a high temporal resolution (e.g. up to 30 s) are extremely valuable for identifying hazardous plumes and can also help to characterise eruptive events (e.g. Gupta et al., 2022; Prata et al., 2022).

This paper uses data from the Infrared Atmospheric Sounding Interferometer (IASI) on the three MetOp satellites and the Advanced Baseline Imager (ABI) on the Geostationary Operational Environmental Satellite (GOES) to analyse the plumes of ash and SO₂ from the 2021 eruption of La Soufrière. Further details on the instruments and methods used to retrieve information about the plumes are given in Sect. 2. The eruption sequence and plume characteristics are discussed in Sect. 3. These instruments provide a complementary view of the eruption, with the high temporal resolution of ABI allowing the eruption sequence to be evaluated, while using IASI it possible to study the dispersion of the plumes as they travel across the globe.

2 Methods

2.1 The Infrared Atmospheric Sounding Interferometer

2.1.1 Instrument

The Infrared Atmospheric Sounding Interferometer (IASI) is a Fourier transform spectrometer on-board three meteorological satellite instruments: MetOp-A, -B and -C launched in 2006, 2012 and 2018, respectively, with data from all three instruments used in this study. The instrument's field of view consists of four circular pixels, each with a 12 km diameter (at nadir) within a 50 km by 50 km square at nadir (Clerbaux et al., 2009). The instruments measure across a wide spectral range within the infrared part of the electromagnetic spectrum between 645 and 2760 cm⁻¹ (3.6–15.5 µm) with a high spectral sampling of 0.25 cm⁻¹ and apodised spectral resolution of 0.5 cm⁻¹ (Blumstein et al., 2004). More information on the IASI instrument can be found in Clerbaux et al. (2009). Within the instrument's spectral range there is sensitivity to volcanic ash (v-shaped absorption feature between 750 and 1250 cm⁻¹; Clarisse et al., 2010a) and SO₂ (three absorption features ν₁, ν₃ and ν₁+ν₃, centred at 8.7, 7.3 and 4 µm, respectively). A number of methods have been developed to extract information about SO₂ from the IASI spectra (e.g. Clarisse et al., 2008; Walker et al., 2011, 2012; Clarisse et al., 2012; Carboni et al., 2012; Clarisse et al., 2014). Retrieval techniques have also been developed for other volcanic gases and aerosols including H₂S (Clarisse et al., 2011), CO (Martínez-Alonso et al., 2012), sulfate (Guermazi et al., 2021) and ash particles (Clarisse et al., 2010a, b; Maes et al., 2016; Ventress et al., 2016; Taylor et al., 2019). Each IASI instrument obtains near-global coverage twice a day, and as they are infrared measurements there is no break in coverage at night and during high-latitude winters. This coverage and sensitivity to gases and aerosols associated with volcanic eruptions makes the IASI instruments well suited for studying the evolution of volcanic plumes. This study uses retrieval schemes developed by Walker et al. (2011, 2012) and Carboni et al. (2012) for quantifying SO₂ emissions from La Soufrière and Sears et al. (2013) for the detection of volcanic ash.

2.1.2 Retrievals

In this study, two methods have been used to analyse SO₂ plumes from the La Soufrière eruption. The first method is a linear retrieval, which is applied in this case to detect pixels that contain elevated quantities of columnar SO₂; see Walker et al. (2011, 2012). The second method is an optimal estimation retrieval scheme, which has been applied to the flagged pixels to quantify the column amount, height, and the effective radiating temperature, and the errors associated with each of these. In broad terms, this method works by comparing the IASI-measured spectra against SO₂ spectra simulated by the fast-radiative transfer model RTTOV (version 9; Saunders et al., 1999); see Carboni et al. (2012, 2016, 2019) for more details. The following changes have been made to the retrieval setup.

• A higher a priori height has been used. In most previous applications, an a priori height of 400 hPa (∼ 7.6 km) has been used. This has been changed to 150 hPa (∼ 14.3 km; upper troposphere) to reflect the higher injection height of the La Soufrière eruption as reported in Global Volcanism Program (2021b). The a priori variance has been kept at 500 hPa.

• In previous applications of the retrieval, the 1σ plume thickness was constrained and set to 100 hPa. This is more appropriate for a plume in the lower troposphere (∼ 1.9 km at 400 hPa), but due to the exponential change in pressure with height in the atmosphere, this thickness is inappropriate in the upper troposphere and lower stratosphere (UTLS) (∼ 4.8 km at 150 hPa). In this application the a priori plume thickness has been set to 30 hPa (∼ 1.2 km at 150 hPa). This is the minimum thickness that can be set given the spacing of the pressure levels in the RTTOV forward model.

The results have been divided into descending (satellite travelling N to S; ∼ 09:30 local overpass time at the Equator) and ascending (satellite travelling S to N; ∼ 21:30 local overpass time at the Equator). This approach can lead to artefacts in the results at the point where the data cross into the next day. This is particularly notable in the ascending orbits at La Soufrière, with the date change occurring around the location of the volcano. To minimise this impact, the descending and ascending results have been offset from each other: the descending results are a 24 h composite of the descending nodes of orbits starting on each date, while the ascending results are compiled from 24 h of ascending nodes of orbits starting from midday on each date and up until midday of the following day. In this way the artefact is moved to the other side of the globe, minimising its impact. In this paper, the ascending results are referred to by the start date.

To estimate the total mass of SO₂ the retrieved column amounts are first gridded to a regular grid of 0.125 by 0.125^∘(roughly 13.5 by 13.91 km at 13^∘ N, 61^∘ W): the retrieved values are interpolated to fill gaps created by the IASI field of view and where pixels fail the quality control. Using the area of each grid box, the column amounts are converted to a mass and the total mass on each day is then computed by summing all the gridded masses within a defined region (−45 to 45^∘ N and −180 and 180^∘ E). The total mass errors are computed in the same way, which may lead to an overestimation of the mass error, but to sum the errors in quadrature could lead to an underestimation due to the systematic errors. Note that this overestimation in the errors is carried forward into the flux and total emission errors. A vertical distribution of the SO₂ mass can be obtained by combining the mass grid with a height grid generated from the retrieved heights.

Error analysis by Carboni et al. (2012) explored the effects of cloud and volcanic ash on the iterative SO₂ retrieval, demonstrating that ash with an optical depth of 1 (at 550 nm) can significantly affect the output of the retrieval. An ash optical depth of 2 was shown to cause a 50 % underestimation in the SO₂ column amount, and an ash optical depth of 5 was shown to mask the SO₂ signal completely. To investigate the effect of ash on the La Soufrière results, a linear ash retrieval as described in Sears et al. (2013) has been applied to the IASI data for the first few days after the first eruptive event to detect pixels containing volcanic ash that may affect the retrieved SO₂ values. This retrieval is run at three pressure levels (400, 600 and 800 hPa) to obtain three estimates of the ash optical depth. Pixels are then flagged as containing volcanic ash if any one of these ash optical depths exceeds a threshold. In this study, the results for the 400 hPa (∼ 7.6 km) level are used as it is the closest level to the retrieved heights. The linear ash retrieval makes a number of assumptions, which means that the retrieved ash optical depths are not necessarily accurate, and it is being used here to simply indicate the possibility of ash affecting the SO₂ retrieval output. Note that there are spectral similarities between volcanic ash, desert dust and desert surfaces that can lead to false ash flags in desert regions such as the Sahara (e.g. Prata et al., 2001; Simpson et al., 2003; Park et al., 2014). To avoid this, the ash linear retrieval has been run for a smaller region (−5 to 25^∘ N, −68 to −20^∘ E) than the main analysis.

Figure 1 Ash optical depths at 550 nm from the IASI linear ash retrieval for 10 April 2021. This retrieval was run assuming a height of 400 hPa. The volcano's location is indicated by a red triangle. The results are a composite of multiple orbits which have been divided into descending (∼ 09:30 local overpass time at the Equator) and ascending (∼ 21:30 local overpass time at the Equator). The ascending results are a composite of orbits from 12:00 UTC on 10 April to 12:00 UTC on 11 April. For further details see Sect. 2.1.2.

Table 1 Percentage of IASI SO₂ pixels on 9 to 12 April 2021 that may be affected by volcanic ash. This has been calculated for the region −5 to 25^∘ N, −68 to −20^∘ E. AOD refers to ash optical depth. Note that the results are split into descending and ascending nodes, with the ascending results consisting of orbits beginning after 12:00 UTC on the date indicated and ending at 12:00 UTC of the following day. See Sect. 2.1.2 for further information.

Figure 1 shows examples of an optically thick ash plume travelling east from the volcano on 10 April 2021. Analysis of the percentage ash detections (Table 1) show that the SO₂ column amounts may be strongly affected by the presence of ash, and so the column amount and mass values presented in this paper should be considered minimum estimates, especially in the first few days after the eruption began (before the ash falls out). Future work should consider the simultaneous retrieval of ash and SO₂ so as to reduce the impact of ash on the SO₂ iterative retrieval. Additionally, SO₂ may be missed where it falls below the retrieval detection limit.

2.1.3 Estimating SO₂ e-folding time and flux

The decrease in the total mass of SO₂ (m) in the atmosphere with time (t) can be described as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}m}{\mathrm{d}t}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}m,\end{array}$$

where τ is the average SO₂ e-folding time. The e-folding time refers to the lifetime of SO₂ in the atmosphere and incorporates the loss of SO₂ due to oxidation and deposition and where the SO₂ amount falls below the detection limit of the instrument being used. Fitting Eq. (1) to the total SO₂ atmospheric burden obtained with a satellite instrument can give a simple estimate of the average e-folding time over a given period. For an eruption like La Soufrière, where there are multiple emission events, the total SO₂ mass loading is a function of both the average e-folding time and variable SO₂ flux (f):

$$\begin{array}{}\text{(2)}& m_i=m_{i-\mathrm{1}}{e}^{-\frac{\mathrm{1}}{\mathit{\tau }}\mathrm{\Delta }t}+f\mathit{\tau }(\mathrm{1}-e^{-\frac{\mathrm{1}}{\mathit{\tau }}\mathrm{\Delta }t}),\end{array}$$

where i is the time step and Δt is the time interval between measurements.

Carboni et al. (2019) uses Eq. (2) within an optimal estimation approach to estimate both the SO₂ flux at each time step and an average e-folding time for the entire eruptive period. This approach has been applied here to the total masses obtained for La Soufrière. The average e-folding estimate can be strongly influenced by the a priori value, and thus for this study an independent estimate of the e-folding time was computed by fitting Eq. (1) to the IASI SO₂ masses computed for the 23 to 30 April 2021 (the period after the last explosive event occurred). This is shown in Fig. 2 and produced an e-folding estimate of 5.47 d. This value was used as the a priori in the optimal estimation approach to generate a time-varying emission flux for the 14 d emission period (9 to 22 April 2021). The optimal estimation scheme produced an average e-folding time of 6.07±4.74 d. On days with no explosive events, as identified with the ABI data (see Sect. 2.2.2), the a priori SO₂ flux and uncertainty are set to 0.

Figure 2 IASI SO₂ masses for 23–30 April 2021, which is the period after the last eruptive event at La Soufrière. A simple and independent estimate of SO₂ e-folding time has been computed by fitting Eq. (1) to the IASI total SO₂ mass estimates. This has been used as an independent a priori in the optimal estimation method for obtaining the flux and average e-folding time for the eruptive period (9–22 April).

To compute the total mass of SO₂ erupted over the studied period, the fluxes were summed and then multiplied by the time difference (Δt) between the descending and ascending orbits (assumed to be 12 h). To estimate the total erupted mass error, the individual flux errors are summed in quadrature (i.e. a quadratic mean of the errors), and then multiplied by Δt. An upper and lower bound of the total emission estimate is obtained by summing the SO₂ flux plus or minus the errors at each time step (excluding negative values), respectively, and multiplying by Δt.

2.2 The Advanced Baseline Imager (ABI)

2.2.1 Instrument

The Advanced Baseline Imager (ABI) is on-board three of the Geostationary Operational Environmental Satellite platforms: GOES-16 (or GOES-East), GOES-17 (or GOES-West) and GOES-18, launched in 2016, 2018 and 2022, respectively. The GOES-16 instrument is used in this study. Recent images produced with this instrument can be accessed at https://www.star.nesdis.noaa.gov/goes/fulldisk.php?sat=G16 (last access: 5 April 2023). The area imaged by GOES-16 (positioned at 75.2^∘ W) covers most of the continental United States, Eastern Canada, Central and South America and the Caribbean, extending across the Eastern Pacific and a large part of the North Atlantic Ocean. The ABI instrument has 16 channels spanning the 0.45 to 13.7 µm spectral range, with the spatial resolution for the different channels varying from 0.5 to 2 km at nadir (Schmit et al., 2005). Having seven channels between 7.3 and 13.3 µm means that the instrument has sensitivity to ash and SO₂, and the instrument has previously been used for the detection of ash and SO₂ and quantification of ash properties (e.g. Pavolonis et al., 2020). Previous versions of the ABI instruments have also been used for this purpose (e.g. Yu et al., 2002; Ackerman et al., 2008). During the eruption of La Soufrière, the ABI sensor recorded a new full disc image every 10 min. In addition, on GOES-16 there are two moveable mesoscale regions, covering an area of 1000×1000 km (at subsatellite point), which can provide data every minute (Schmit et al., 2017). These are moved to provide higher temporal coverage for events such as severe weather, hurricanes and forest fires (Schmit et al., 2017). During the La Soufrière eruption, one mesoscale region was moved over the volcano. This started at 09:00 UTC on 10 April (missing the earlier eruptive events) and ended at 06:00 UTC on 16 April 2021 (before the last three eruptive events). In this study, the ABI data have been used to identify the start and end times of each eruptive event and for determining the height of the ash cloud. When available, the 1 min mesoscale data have been used to identify the start and end times of the eruptive events. The 10 min resolution data have been used for the remaining days. The 10 min data have been used for the height analysis. While the ABI instrument has sensitivity to ash and SO₂, it has not been used to quantify the amount of either in this study. Instead, the ABI data have been used to document the sequence of explosive events produced by the 2021 La Soufrière eruption.

2.2.2 Analysis of ash

The Satpy Python package has been used to load true- and false-colour products from the ABI data for the eruption period. In both cases, the RGB components of this were then normalised and used to produce maps showing the volcanic plumes from this eruption (9–22 April). The false-colour images have been constructed by assigning the 12.3–10.3 µm, 11.2–8.4 and 10.3 µm channels to red, green and blue, respectively (for further information on the ash RGB product, see Fuell, 2018). Table 2 shows how to interpret the images. An example of these plots can be seen in Fig. 3. In addition, the 11.2–12.3 and 10.3–11.2 µm brightness temperature differences (BTDs) have been calculated. The 11.2–12.3 µm BTD combination has been widely used to distinguish volcanic ash in satellite data and uses the positive transmission gradient between 10 and 12 µm caused by volcanic ash (Prata, 1989a, b). The 10.3 µm channel is less affected by water vapour (Lindsey et al., 2012), and thus the 10.3–11.2 µm BTD may be more appropriate in tropical atmospheres. A threshold can be used to flag pixels containing ash. The BTD approach can be limited for a number of reasons; for example, false detections can occur in regions with high surface emissivities (e.g. deserts) and strong temperature inversions or due to desert dust, while ash clouds may not be identified if they are too optically thin or thick or if there are significant quantities of ice or water within the plume (Rose et al., 1995; Simpson et al., 2000; Prata et al., 2001). There is also no generally applicable threshold that can be set for ash detection, with the most appropriate value varying with region and time.

Figure 3 Examples of results from the ABI instrument for the La Soufrière plume at 12:00 UTC on 10 April 2021. (a) True-colour image. (b) False-colour image (12.3–10.3, 11.2–8.4 and 10.3 µm channels assigned to red, green and blue, respectively). (c) Map showing 10.3–11.2 µm BTD. (d) Map showing 11.2–12.3 µm BTD. Animations for each event can be found in Taylor et al. (2023).

Table 2 Outline of how to interpret the ABI false-colour images (see Fuell, 2018 for more information but note that differences in how the data is normalised may lead to variations in colour). An example image is shown in Fig. 3.

Together, the true, false and BTD images have been used to study the evolution of the plumes from La Soufrière. Careful examination of the plots allowed the identification of the approximate start and end times of each eruptive event. It should be noted that the start and end times reported here will be different to those observed on the ground. Here, the start time is the first time that the plume is observed in the ABI data, while the end time incorporates the time taken for the plume to rise and disperse away from the volcano. For this reason, the times reported here are different to those reported in Sparks et al. (2024) using seismic data. Additionally, the ABI analysis is limited by the 10 or 1 min resolution of the instrument and can be complicated by the presence of cloud. The interpretation of the start and end times is a subjective process, with the end time being particularly challenging to establish. In general, the end time was determined when the plume moved away from the volcano. This can be affected by the wind speed and so the end times should be considered approximate. As determining the end time was especially challenging, an end time range has been given for most events. It is also likely that lower-magnitude explosive events or degassing between events may not be identifiable in the ABI data. Note that for the full disc the measurement time over La Soufrière has been computed based on the latitude of the volcano and the measurement start and end times, approximately 243 s from the start time. This has been used to report the approximate times for the full disc results in this study. No such adjustment has been made to the mesoscale results, where the 1 min temporal resolution ensures a higher accuracy.

An estimate of the ash height is obtained by comparing the brightness temperature in the 11.2 µm band (ABI channel 14) for each eruptive event in a 0.1^∘ box around the volcano with the ERA-5 ECMWF temperature profile interpolated to the volcano's location (European Centre for Medium-Range Weather Forecasts, 2021). Here this is referred to as the “brightness temperature (BT) method”. The measured brightness temperature is representative of the temperature at the top of the ash or cloud layer and can thus be used to approximate the height. The 11.2 µm channel is within an atmospheric window that has little sensitivity to water vapour. The minimum (or coldest) value within the 0.1^∘ box around the volcano helps to select the optically thickest part of the plume, and therefore the pixel that is least likely to be affected by radiation from beneath the plume. The 0.1^∘ box size should remove any significant effect due to parallax. In its application here, only one value is returned for each eruptive event (referred to as the “optically thick height method/solution”), which ensures the result from the most optically thick part of the plume is reported but subsequently does not show any variability within the event. This method has been widely used to estimate the ash cloud top height (e.g. Prata and Grant, 2001). There are, however, a few limitations that have been discussed in several papers (e.g. Oppenheimer, 1998; Prata and Grant, 2001; Zakšek et al., 2013) and which are highlighted in Table 3. Despite the various limitations of this method, it is a quick and frequently used approach for estimating ash cloud heights. For this application, an uncertainty of 1 K has been assumed for the instrument measurement encompassing the instrument noise and gaseous absorption above the cloud. Horváth et al. (2022) also applied the BT method to ABI data for La Soufrière for comparison against GOES-17 side view heights and noted a number of limitations with the BT method, including an underestimation of the heights for smaller eruptions due to a warm bias. The results here will likely share these limitations and thus should be treated with caution. However, Horváth et al. (2022) noted that the side view and BT method heights showed better agreement for cold plumes or for plumes spreading around the height of the tropopause.

(Oppenheimer, 1998; Zakšek et al., 2013)(Glaze et al., 1989; Zakšek et al., 2013)Oppenheimer (1998)(Woods and Self, 1992; Holasek et al., 1996a)(Oppenheimer, 1998)(Holasek and Self, 1995)(e.g. Prata et al., 2022)

Table 3 Summary of the sources of error associated with estimating the height of volcanic plumes by comparing the brightness temperature at 11.2 µm (ABI channel 14) with a temperature profile. In this paper this is termed the “brightness temperature (BT) method”. These limitations are discussed in several papers (Oppenheimer, 1998; Prata and Grant, 2001; Zakšek et al., 2013).

In this study, there was commonly a tropospheric and stratospheric solution when using the BT method. Additional information has been obtained by comparing the plume direction and speed with the ERA-5 ECMWF wind profiles. Figure 4a shows the average and standard deviation of the wind direction with height, and the average wind speed with height is shown in Fig. 4b. At heights of less than 5 km, the wind is primarily travelling to the west (e.g. 267.5^∘ at 4.3 km) with wind speeds of less than 10 m s⁻¹. At around 5 km, the wind direction shifts to the east. Between 7 and 17 km the wind direction is fairly consistent (varying between 101 and 126^∘) before shifting back to the west at 19 to 20 km. The easterly winds in the stratosphere are a characteristic associated with a phase of the quasi-biennial oscillation (QBO): alternating strong easterly and westerly zonal winds around the Equator that propagate through the stratosphere to the tropopause (Reed et al., 1961; Baldwin et al., 2001). The wind speed is shown to increase between 5 and 14 km and then fall between 15 and 20 km. Figure 4c and d show the average plume direction and speed. This has been computed by visually identifying the centre of the optically thick plume in the RGB image 1 h or 30 min after the start of each eruptive event (as identified with ABI). For the first 6 d of the eruption, where there is some uncertainty with regard to the height, the plume is shown to broadly be travelling to the east, supporting a tropospheric solution. To take this analysis further, at the location identified as the centre of the optically thick plume, the BT method was used to obtain a second set of height solutions (referred to as the “BT method-centre” solutions). The wind directions at the heights obtained with both BT methods are compared against the horizontal velocity components from the motion of the plume (based on distance travelled in 1 h or 30 min, as above) to assign a degree of confidence to the height results for each eruptive event. Figure 5 shows examples from the first three eruptive events, which were shown to be tropospheric, uncertain and stratospheric, respectively.

There are differences between the heights obtained with the “BT method” applied in this study and the results obtained using the BT method found in Horváth et al. (2022), which may arise from different sampling approaches. This is also the case for the BT heights presented in Sparks et al. (2024) where differences may also arise due to the different meteorological data and channels used.

Figure 4 (a) The average and standard deviation of the ERA5 ECMWF wind direction for different height levels between 9 and 22 April. (b) Average wind direction profile of the ERA5 ECMWF wind direction for different height levels between 9 and 22 April. (c) Plume speed for the different explosive events. (d) Plume direction of travel for the different explosive events.

Figure 5 Demonstration of the use of the ERA5 ECMWF wind data to determine whether the tropospheric or stratospheric solution is more likely. Examples are from the first three eruptive events (see Table 4). These were classified as tropospheric, inconclusive and stratospheric, respectively.

3 Results

3.1 Eruption sequence

True- and false-colour maps and ash BTDs were produced from the ABI data for the entire eruptive period (9 to 22 April 2021) using the instrument's high temporal resolution. An example from 10 April at 12:00 UTC is shown in Fig. 3. Animations showing the true- and false-colour images and the BTD maps for each eruptive event can be found at https://ora.ox.ac.uk/objects/uuid:ca4e7a94-37c6-4d5f-b94a-1f287e661f8a (Taylor et al., 2023, last access: 2 October 2023). These images have been carefully examined to get start and end times for each of the explosive events identified. The start and end times are reported in Table 4. The heights obtained with the BT method are reported in Table 5.

Using the ABI data, a minimum of 35 explosive events containing volcanic ash were identified. The number identified is limited by the temporal resolution of the instrument and cloud cover that might obscure the plume, and thus could potentially be higher. For example, explosive events that last less than 10 min may not be identified in the 10 min resolution data, while discrete eruptions or pulses that occur within the 1 or 10 min between images may be classified as the same eruption. Times where this is suspected are mentioned in Table 4. Additionally, lower-level activity may not be visible within the ABI data. Figure 6a shows the number of explosive events that started on each day as determined with the ABI data. This is shown to be highest on 10 April and decreases throughout the remainder of the eruptive sequence, with Fig. 6c showing that the repose time between explosive events increases over time. The duration of the eruptive events (Fig. 6b) is also shown to generally increase throughout the 9–22 April eruptive period, as also noted by Joseph et al. (2022).

Table 4 Summary of plumes from La Soufrière seen with the ABI satellite instrument. The start and end times have been determined with the ABI images. These do not necessarily correspond to activity observed on the ground: the start time refers to when the plume is first visible in the ABI data, while the end time incorporates the time taken for the plume to rise and disperse away from the volcano. For a number of reasons discussed in Sect. 2.2.2, the times presented here should be considered approximate. Due to the uncertainty with the end times, in some cases these are presented as a range. The eruption has been split into four phases based on the general character of the events. The phase numbers correspond to (1) the initial explosive event, (2) sustained ash emission, (3) the pulsatory phase and (4) the waning phase.

Table 5 Heights retrieved using brightness temperature (BT) method applied to ABI data.

Figure 6 Summary of ABI results. (a) Number of explosive eruptions on each day. (b) Duration of each explosive event. (c) Repose duration between events. (d) Heights retrieved with optically thick and centre BT methods. (e) Using the ABI data, the 2021 eruption has been divided into the four phases with the timing shown here.

The eruptive events identified from the ABI data have been split into four phases, each with distinct characteristics. The first emission (and first phase of the eruption) is likely to correspond to arrival of gas-rich magma at the surface and loss of the lava dome cap. This is seen with ABI at around 12:54 UTC (08:54 LT) on 9 April and corresponds well with the 08:41 LT (12:41 UTC) eruption start time stated in Joseph et al. (2022). Using the ABI data, this emission appears to be relatively short lived, lasting around 4.3 to 4.8 h, finishing at around 17:14–17:44 UTC (although significant cloud cover makes the end time difficult to distinguish). During this time two plumes are evident: a low-altitude plume can be seen travelling to the west from the volcano, while the main plume travels to the east–northeast, with a height of around 13.11${}_{-\mathrm{0.12}}^{+\mathrm{0.12}}$ km based on the optically thick BT method. This is greater than the 8 km reported by the Belmont Observatory (Global Volcanism Program, 2021b). However, differences between ground-based and satellite measurements can arise for a number of reasons, including the different viewing angles and assumptions, poor visibility on the ground, the plume rising too high for the height to be accurately obtained from the ground, and different definitions of plume height (e.g. Tupper et al., 2004; Tupper and Kinoshita, 2003; Tupper and Wunderman, 2009).

Following a short pause of around 1.5 to 2 h, the second phase of the eruption started at around 19:14 UTC on 9 April. This is characterised by a sustained ash emission event with a duration of around 9.5 to 9.67 h. Height estimates based on the BT method are between 15.32 and 19.20 km. In this case, the plume direction and wind profile do not help determine whether the tropospheric or stratospheric solution is more likely. The structures of the plume suggest that there may have been some strong pulses within this sustained phase (see Table 4), but these cannot be easily resolved from the 10 min ABI images. The observations using the ABI instrument are generally consistent with those made by Joseph et al. (2022), who mention a “sustained” and “pulsing” explosive phase starting at 16:00 UTC on 9 April and finishing at 06:00 UTC on 10 April.

The explosive events identified with the ABI instrument suggest that the eruption entered a third phase at 04:54 UTC on 10 April. This phase appears more pulsatory, with a further 25 discrete ash emission events identified within a 65.3 h period. The use of the 1 min ABI data during this phase greatly helped to identify the event timings: highlighting the considerable advantage of this instrument. Joseph et al. (2022) characterise this phase as one of “ash venting”. The length of these explosive events identified from the ABI data varied from around 16 min to up to 4.65 h. This style of repetitive ash bursts is similar to the “Vulcanian” explosions documented at multiple andesite volcanoes, which is considered to reflect the rapid sealing, pressurisation and failure of magma within the conduit (e.g. Druitt et al., 2002; Watt et al., 2007; Joseph et al., 2022). The first pulse is shown to most likely be in the stratosphere with height solutions of 18.89${}_{-\mathrm{0.19}}^{+\mathrm{0.14}}$ and 19.20${}_{-\mathrm{0.14}}^{+\mathrm{0.16}}$ km for the optically thick and centre BT methods, respectively. The majority of the remaining pulses in this phase have been labelled as most likely tropospheric or uncertain. The height results for each event are shown in Table 5.

The fourth and final phase of the eruption, as defined with the ABI data, began on 13 April at 10:36 UTC and continued to the 22 April. During this phase there were a further eight distinct events, each lasting a few hours, with generally longer repose times than the previous phases and with each event emitting ash into the troposphere (heights between ${\mathrm{2.4}}_{-\mathrm{0.17}}^{+\mathrm{0.18}}$ and 14.42${}_{-\mathrm{0.19}}^{+\mathrm{0.22}}$ km). Figure 6d shows that during this phase the heights for each eruptive event fall, suggesting a decrease in eruptive power. Compared to the earlier phases, the RGB images suggest that the plumes in this waning stage are more optically thin and it is possible, especially for events 29 and 33–35, that the heights obtained are underestimates due to upwelling radiation from beneath the plume. A number of the explosive events in this phase appear to have an explosive eruption which produces a small optically thick cloud, followed by a longer emission of a less optically thick ash cloud, while the last three events consist of small plumes with ash that is optically thin. The final eruptive event emitted ash into the lower troposphere for around 6.3–7.3 h.

In contrast to this study, Horváth et al. (2022) count 49 explosive events between 9 and 22 April using the ABI data. The higher number reported by Horváth et al. (2022) is partly because they have divided event 2 (as termed in this paper) into multiple explosive events. Sparks et al. (2024) also document the sequence of eruptive events occurring during the April 2021 eruption. Their analysis is primarily based on seismic data, but they also make note of the first observation times with the ABI instrument. As expected, there are differences between the start times identified with the seismic and ABI instruments, as there is some time before the plume becomes visible with the satellite instrument. On the whole, their ABI observation times (only start time reported) agree well with the ones presented here; however, there are a few differences. Firstly, there are differences in the timings of the events: they used the 10 min full disc rather than the mesoscale data, and other differences may arise from the subjective nature of this exercise. Secondly, this study identifies some events that are not included in Sparks et al. (2024) time series, while Sparks et al. (2024) is able to split event 2 (as referred to here) into multiple events using the seismic data. Finally, Sparks et al. (2024) also identify four phases to the eruption, and while there is some overlap with the phases outlined here, there is disagreement with the timings reflecting the different datasets and metrics used.

3.2 SO₂ dispersion

The linear and iterative SO₂ retrievals have been applied to IASI spectra from 9 to 30 April 2021. The iterative retrieval column amounts for 9 to 26 April are displayed in Fig. . Figure 8 shows examples of the column amount and height outputs from the iterative retrieval for 10 and 11 April. Animations of the iterative retrieval column amount and height outputs for the full time period explored can be found in the Supplement along with animations showing the error outputs. The mapped results have been obtained by gridding the outputs for all three IASI instruments to a regular 0.125 by 0.125 degree grid.

Figure  shows the SO₂ column amounts from the iterative retrieval. In these maps it is possible to see the evolution of the La Soufrière SO₂ plumes between 9 and 26 April. The plume is first apparent as a faintly elevated emission to the east of the volcano in the ascending overpass on 9 April 2021 (incorporates ascending orbits from 12:00 UTC 9 April to 12:00 UTC 10 April). A stronger signal is then visible in the descending orbits on 10 April fanning out to the east of the volcano across the North Atlantic. The general east and southeastward transport of the plume between 9 and 11 April implies that the bulk of SO₂ has been emitted into the troposphere, with Fig. 4a indicating wind directions between ∼ 90–140^∘ dominating in the troposphere between 8 and 17 km. By the ascending orbits on 11 April (incorporates ascending orbits from 12:00 UTC 11 April to 12:00 UTC 12 April), the plume has travelled across the Atlantic and has almost reached the western coast of Africa. The higher values for the column amounts leading back to the source reflect the frequent eruptive events that occurred over this time period. The general eastward transport of the plume is consistent with observations by Babu et al. (2022), who tracked the plume with the Ozone Monitoring Instrument. On 12 April, while the bulk of the plume is still advancing towards the east (Fig. ), a fraction of the plume travels to the south and west of the volcano. The wind directions shown in Fig. 4a imply that for southward or westward transport, either some SO₂ has been emitted into the lower parts of the troposphere (less than 5 km) or more likely some SO₂ has been emitted or lofted into the stratosphere. Over the next few days the plume continues to be transported both to the east and west of the volcano, and by the 21 April the plume has been transported around the circumference of the globe. Note that there is an emission of SO₂ from an eruption at Sangay in Ecuador from 12 April (Global Volcanism Program, 2021d). This combines with the plume from La Soufrière, and the two cannot be easily distinguished from each other from 13 April.

Figure 9b shows the distribution of the SO₂ mass with latitude. The largest concentration of SO₂ (up to 17 April) is between 0 and 20^∘ N, but parts of the plume have reached 45^∘ north and south. This highlights the ability of tropical eruptions to transport SO₂ across both hemispheres, rather than being confined to a narrower latitude band like plumes from eruptions at high latitudes (Schmidt and Black, 2022). Also highlighted in Fig. 9 is a line of elevated SO₂ around 15^∘ S. This is the latitude of the Sabancaya volcano (15.787^∘ S, 71.857^∘ W, summit elevation of 5.96 km) in Peru, which was erupting throughout the studied time period. Plumes from Sabancaya can be distinguished in the iterative retrieval maps in the Supplement.

Figure 7 Continued.

Figure 8 Examples of the IASI SO₂ iterative retrieval gridded column amounts and heights for scenes from 10 and 11 April 2021. The date and orbit direction are indicated on the left-hand side of the figure. Maps of the results for the entire time period studied can be found in the Supplement.

Figure 9 Summary of the IASI SO₂ results generated from iterative retrieval outputs for 9 to 30 April 2021. Each of these has been computed for the −45 to 45^∘ N and −180 and 180^∘ E region. (a) Vertical distribution of the SO₂ mass across the plume with time. The dashed lines show the mean height of the tropopause and 1 standard deviation from that point. (b) Distribution of the SO₂ mass with latitude. The black triangle at 13.33^∘ N indicates the latitude of La Soufrière. The smaller triangle at 15.79^∘ S shows the latitude of Sabancaya in Peru, which also erupted during the studied time period. (c) Time series of the total mass of SO₂ with errors. Using the retrieved SO₂ heights and the tropopause heights, this has been divided into tropospheric and stratospheric SO₂. The height errors have also been used to create two more categories: heights that are stratospheric but given the height error could be tropospheric and heights that are tropospheric but could be stratospheric. (d) The flux of SO₂ with time and the number of eruptive events identified from the ABI data on each day.

3.3 SO₂ mass and flux

A total mass time series, derived from the IASI iterative retrieval output for the −45 to 45^∘ N and −180 and 180^∘ E region, is shown in Fig. 9c. Note that other volcanoes (e.g. Sabancaya in Peru and Sangay in Ecuador) were erupting at this time, and any SO₂ from these eruptions entering the region will affect the total mass, e-folding and flux estimates. However, given that these are smaller emissions than La Soufrière, their impact is negligible and within the reported errors. The IASI-derived SO₂ mass loading peaks at 0.31 ± 0.09 Tg in the descending orbits on 13 April. This is fairly small compared to other eruptions studied by IASI, including Nabro (1.6 Tg), Kasatochi (0.9 Tg), Grímsvötn (0.75 Tg), Copahue (0.72 Tg) and Sarychev (0.6 Tg) in 2011, 2008, 2012 and 2009, respectively (Carboni et al., 2016). It is more similar to the maximum mass loading recorded from Alu Dalafila (0.2 Tg) in November 2008 (Carboni et al., 2016).

The fact that the total SO₂ mass loading for La Soufrière peaks a few days after the 2021 eruption began could be linked to the multiple explosive events which occurred between 9 and 13 April, leading to an increase in the total atmospheric SO₂ loading: SO₂ emitted at a rate greater than the existing SO₂ declines. The Global Volcanism Program (2021c) reported an early satellite estimate of the SO₂ mass of 0.4 Tg on 10 April: larger than the estimate from IASI on the same day 0.24±0.02 Tg, ascending result). Differences between instruments and retrievals are expected due to the varying sensitivities to SO₂, different height assumptions, the differing effects of volcanic ash and aerosols, and different spatial resolutions and observation times. Carn et al. (2016), for example, note that there are often differences between the results from UV and infrared sensors. As mentioned in Sect. 2, the IASI iterative SO₂ retrieval column amounts, and subsequently the total masses, are affected by volcanic ash and instrument sensitivity, and thus the column amounts and total mass values presented here should be considered minimum estimates.

From 13 April the total mass of SO₂ is shown to fall as the SO₂ is removed from the atmosphere through deposition by conversion to sulfate aerosol or dilution below the detection limit of the instrument. The e-folding time used here (see Sect. 2.1.3) describes this loss process. This varies with a number of factors including the latitude of the volcano, the injection height of the plume, meteorological conditions, cloud cover, water vapour, season and the presence of ash (Carn et al., 2016; Zhu et al., 2020, 2022; Schmidt and Black, 2022). Typically, the e-folding time varies from hours to days in the lower troposphere to weeks in the stratosphere. A first estimate of 5.47 d for the SO₂ e-folding time was estimated by fitting Eq. (1) to the IASI total SO₂ masses between 23 and 30 April 2021. Following the application of the Carboni et al. (2019) method to compute the flux and average e-folding time, the average e-folding time was adjusted to 6.07±4.74 d. This is in line with other eruptions including Jebel at Tair (2007) and Merapi (2010), both of which are located in the tropics, emitted plumes between 15 and 18 km, and had e-folding times of between 2 and 4 d (Carn et al., 2016). Given that there was ash emission during the La Soufrière eruption, it is possible that the e-folding time could have been reduced as a result of accelerated oxidation of SO₂ due to reactions on the ash surface, as was seen for the Kelut eruption in a study by Zhu et al. (2020). Some of the SO₂ emitted during the La Soufrière eruption was converted to sulfate aerosol, as is shown in Babu et al. (2022) and Bruckert et al. (2023).

The SO₂ masses from IASI were converted to flux estimates following the method described in Sect. 2.1.3 (note that underestimates in the total masses may affect these results). The daily SO₂ fluxes computed from the IASI results are shown in Fig. 9d along with the number of explosive events which started on each day as determined with ABI data. Both are also reported in Table 6. There is generally a good correspondence between the number of eruptive events and the SO₂ flux. The greatest flux (0.24±0.02 Tg d⁻¹) from the IASI data occurs on 10 April in the descending overpasses which incorporates the second eruptive phase (sustained emission), and the beginning of the pulsatory phase. Following this, the SO₂ flux derived from IASI data is shown to decrease over the third eruptive phase that ends on 12 April. The flux on 12–13 April may incorporate some of the SO₂ emitted by Sangay in Ecuador. The flux remains low for the final phase of the eruption and no emission is visible following the eruption on 22 April. These fluxes agree reasonably well with the flux range reported in Joseph et al. (2022), which varied between 2.76×10⁵ t d⁻¹ (0.276 Tg d⁻¹) on 10 April to 331 t d⁻¹ (0.0003 Tg d⁻¹) on 22 April.

By summing the retrieved IASI SO₂ fluxes and multiplying by the time step between the images (assumed to be 12 h) it is estimated that in total this eruption emitted 0.63±0.5 Tg of SO₂. The error is likely to be overestimated due to the way in which the total mass errors are computed (see Sect. 2.1.2). An alternative way of computing the minimum and maximum total emitted SO₂ is to add and subtract, respectively, the flux and error for each time step and then multiply by the time step rather than sum the errors in quadrature. This gives a minimum and maximum estimates of the total SO₂ mass of 0.24 and 2.43 Tg, respectively.

Table 6 Number of eruptive events identified with ABI and IASI SO₂ fluxes for 9 to 22 April 2021. Note that these estimates may be affected by other volcanoes that erupted during this period.

3.4 SO₂ plume heights

Examples of the iterative retrieval height results are shown in Fig. 8. The heights in the descending orbits on 10 April are variable. The heights in the northern edge of the plume (mostly 10 to 13 km, some less than 10 km) are lower than the rest of the plume. The southern part of the SO₂ plume is higher, primarily between 13 and 18 km, and in some cases exceeds this. Similar structures can be seen in the following few days. There is no obvious gradient to the heights in the upper part of the plume as can be seen in IASI SO₂ height measurements shown in the supplementary material of Koukouli et al. (2022) for 10 and 11 April (using the Clarisse et al., 2014 method). Koukouli et al. (2022) show heights increasing from south to north of the plume, which matches well with the wind directions shown in Fig. 4a. The IASI retrieval used here relies on temperature and water variations in the atmosphere, which do not vary as significantly around the tropopause and may affect the results. Multiple retrieval setups were explored (including varying the retrieval first guess height and varying the plume thickness), but the results were similar in each case. Nevertheless, there is a broad agreement with Koukouli et al. (2022), which reports average heights of 15.7±1.16 km for IASI using Clarisse et al. (2014) retrieval and 14.94±3.87 km from TROPOMI (Hedelt et al., 2019 method) (note that these averages are based on a subset of the plume). Additionally, there is agreement with the 13 to 15 km injection heights obtained by Esse et al. (2024) using backwards trajectory modelling. The IASI SO₂ heights are also generally consistent with ash height estimates reported in Global Volcanism Program (2021b) bulletin report. However, it should be noted that SO₂ and ash are not always co-located spatially, temporally or in height (e.g. Holasek et al., 1996b; Thomas and Prata, 2011; Moxnes et al., 2014; Prata et al., 2017). The IASI results presented here are similar to observations made by the Multi-angle Imaging Spectroradiometer (MISR) onboard the Terra satellite (Yue et al., 2022). Results from this instrument on 10 April showed that there are multiple layers of ash, with two primary layers at 12 and 18 km and some parts of the plume reaching 20 km (Yue et al., 2022). Horváth et al. (2022) analysed 30 of the eruptive events using the ABI side view and found that most of the plumes rose to around the height of the tropopause (16 to 17 km) or entered the lower stratosphere (18 to 20 km). They also reported the heights of the overshooting tops, which went up to 23 km for the largest explosive events. Finally, there is agreement the IASI SO₂ heights here and those from the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) and ICON-ART model that show a volcanic plume layer between 15 and 20 km on 13 April (Bruckert et al., 2023).

The IASI results show that the structures within the plume evolve as the plume stretches across the Atlantic, with the higher values (15 to 17 km) making up the majority of the image. Note that over time there appears to be more speckle in the height results (see animation in the Supplement), which reflects increased uncertainty. The variations in height across the entire SO₂ plume are likely to be linked to the multiple injection heights and pulses which occurred during the earlier eruptive events, as discussed in Sect. 3.1. The lower heights obtained from the ABI data and mentioned in the Global Volcanism Program (2021b) bulletin report, for the later stages of the eruption (e.g. 14–22 April), are not reflected in the IASI results. This is potentially because the retrieval has been setup in a way that is optimised for higher plumes, the fact that these eruptive events emitted less SO₂ (below the detection limit of the retrieval) or due to multiple layers of SO₂.

The evolution of the vertical distribution of the SO₂ mass across the plume is shown in Fig. 9a. This is computed from the gridded column amounts (converted to mass values) and the gridded heights for each day (descending and ascending). On 10 April, SO₂ is shown to be emitted at multiple heights: a smaller emission between 11 and 13 km and a larger quantity between 15 and 18 km. The lower layer largely disappears within 3 d. Between 12 and 18 April the SO₂ mass is shown to be spread over 13 to 19 km. Lines showing the average tropopause height and one standard deviation from this are displayed in Fig. 9a. This was estimated based on the World Meteorological Organisation's (WMO) definition of the tropopause, i.e. that the tropopause occurs at the lowest level where the temperature lapse rate falls beneath 2 K km⁻¹ for at least 2 km. Note that this is computed from the tropopause heights for the pixels containing SO₂ and where the iterative retrieval output passes the quality control: this means that the number of pixels used in this calculation varies over time. Much of the SO₂ appears to be in the upper troposphere and around the height of the tropopause, which supports the observation made by Carboni et al. (2016) that SO₂ from explosive eruptions often ends up at the tropopause. This observation is consistent with the direction of travel of the plume (to the east) and the wind direction at these heights (which has an easterly direction in the upper troposphere and around the tropopause, Fig. 4a). A fairly persistent lower layer of SO₂ can be observed between 3 and 7 km. Some of this can be attributed to emissions at Sabancaya volcano in Peru. From 19 April a layer can also be observed at 14 to 14.5 km. This is around the height of the a priori and so is most likely due to a loss of height information in the IASI spectra as the plume thins.

Figure 9c shows the fraction of SO₂ which has been categorised as either tropospheric or stratospheric. Note that this has been calculated on a pixel-by-pixel basis using the height of the tropopause. Also shown in this plot is the degree of uncertainty with regard to the height of SO₂ relative to the tropopause. Figure 9c highlights that a fraction of the SO₂ mass which while classified as stratospheric could in fact be tropospheric when the height errors are considered. Similarly it shows that the fraction of the mass classified as tropospheric which could be stratospheric. In general it appears that there is a larger amount of SO₂ in the troposphere. For example, in the descending orbits on 13 April, 68.88 % and 31.12 % of the SO₂ was shown to be in the troposphere and stratosphere, respectively. However, based on the error bars, 67.78 % of the tropospheric result could also be classed as stratospheric given the overlap of the error bars with the tropopause height, and 98.65 % of the SO₂ classed as stratospheric could be tropospheric.

4 Comparison with the 1979 eruption

There are a number of similarities between the 1979 and 2021 eruptions at La Soufrière. Shepherd et al. (1979) give a detailed account of the eruptive sequence for the 1979 eruption and report a series of explosive eruptions between 13 and 26 April 1979, with the first explosions occurring from within a water-filled summit crater lake. In the first few hours after the eruption began on 13 April 1979 (09:30 to 15:00 UTC), continuous emissions of steam were punctuated with explosive eruptions producing ash clouds that reached greater than 8 km (Shepherd et al., 1979). This activity was followed by a series of six explosive events with varying durations between 20:05 on 13 April and 15:50 on 14 April. Estimates from satellite data on 13 and 14 April give ash heights of 17 to 18 km for two of these plumes (Shepherd et al., 1979). Following this, between 14 and 17 April, the volcano entered a new phase with only occasional ash emissions rising 1 to 2 km above the vent. Another explosive event occurred on 17 April with ash emissions reaching 18.7 km. Following this, the intensity of activity generally declined, although two more explosive events occurred on 22 and 26 April (Shepherd et al., 1979). By the end of the 1979 explosive sequence, the crater had been infilled with pyroclastic debris, and the summit lake had disappeared. The 2021 explosive eruption sequence shares a number of similarities with the 1979 eruption, even though the early stages of the 2021 eruption were rather different, i.e. 3 months of dome effusion into a “dry” summit crater. Both the 1979 and 2021 eruptions consisted of a series of explosive events and both had a pulsatory phase in the first few days after the start of the eruption. However, it is notable that a larger number of discrete events occurred during the 2021 eruptive period. The later eruptive events in both the 1979 and 2021 eruptions also had a longer repose time.

Another similarity between the 1979 and 2021 eruptive periods is the plume height. Measurements of the 1979 heights were far sparser, with estimates coming from relatively infrequent satellite data, ground- and ship-based observations, and aircraft observations. As mentioned above, Shepherd et al. (1979) gives height estimates for individual eruptive events ranging between 8 and 18.7 km, similar to the heights reported here for the April 2021 eruption. This is similar to this study, which reports a fairly wide range of heights for the 2021 eruption. Krueger (1982) used SMS-1 satellite data to estimate that part of the plume exceeded 18 km, placing it in the stratosphere. Another study using an airborne lidar by Fuller et al. (1982) showed layers of ash between 16 and 19.5 km. These heights for ash are in a similar range to those for SO₂ reported in this study, although again the difficulty of using ash as a proxy for SO₂ is noted. While both Krueger (1982) and Fuller et al. (1982) indicate that parts of the plume reached the stratosphere, it is possible that there was significant variability across the full extent of the plume that was unmeasured, and so it is difficult to quantify how much of the plume was tropospheric and stratospheric for the 1979 eruption.

A significant difference between the two eruptions is the quantity of SO₂ observed by satellite instruments. The 1979 eruption of La Soufrière is one of the first observations of volcanic SO₂ made by the Total Ozone Monitoring Instrument (TOMS), but elevated emissions were only seen on 2 different days (Carn et al., 2003). In total, only ∼ 3 kt of SO₂ was measured with TOMS (Carn et al., 2016), which is substantially lower than the estimates reported here from IASI for the 2021 eruption. Carn et al. (2016) indicates that the TOMS SO₂ estimate is only a fraction of the total emitted with Scaillet et al. (2004), obtaining a total sulfur mass of 0.46 Tg based on petrological methods (∼ 1 Tg SO₂; Carn et al., 2016). Carn et al. (2016) suggests that the majority of sulfur emission may have been in the form of H₂S, with slow oxidation preventing SO₂ detection with TOMS. It is also possible that part of the low SO₂ mass estimate from TOMS could be attributed to the lower resolution of the instrument compared to the satellite instruments observing the 2021 eruption.

5 Conclusions

This study has used satellite data to gain insights into the chronology of events and the various emissions of ash and SO₂ during the April 2021 eruption of La Soufrière. Using data from ABI it is possible to identify the emission of ash, and in this study the instrument's high temporal resolution (up to 1 min) was used to identify and characterise different phases of the April 2021 La Soufrière eruption. Careful examination of true- and false-colour images and BTD maps allowed the identification of the start and end times of 35 explosive events that took place during the eruption; however, this should be considered a minimum estimate. Based on these observations, the following four eruptive phases have been identified: (1) an initial explosive event producing an ash cloud; (2) a strong and sustained ash emission lasting several hours and producing a large optically thick ash cloud; (3) a pulsatory phase with 25 discrete explosive events in a 65.3 h period, varying in duration between 16 min and 4.65 h, each producing optically thick ash clouds; and (4) a waning phase with eight eruptive events with lower plume heights, optically thinner ash clouds, and greater durations between explosive events.

Measurements with retrievals developed for the IASI satellite instrument were able to track the plumes of SO₂ from the 2021 La Soufrière eruption for several weeks after the first eruption occurred. A peak value of 0.31 ± 0.09 Tg was obtained for the total SO₂ mass burden on 13 April. However, it is likely that this is an underestimate due to the effects of volcanic ash on the results and the presence of SO₂ below the detection limit of the instrument. Using the daily IASI SO₂ fluxes it was possible to estimate the total emission of SO₂ to the atmosphere, which was 0.63±0.5 Tg with an upper and lower limit of 2.43 and 0.24 Tg, respectively (considering the full extent of the uncertainties on each day). Again, due to the potential underestimation of the mass, these should be considered minimum estimates. The retrieved heights show interesting structures in the plume that reflect the multiple eruptive pulses; this is consistent with observations on the ground and made by other satellite instruments. The height results suggest that the majority of SO₂ was emitted into the upper troposphere and lower stratosphere, with most of the gas being concentrated between 13 and 19 km.

The 2021 eruption showed some remarkable similarities with the 1979 eruption in terms of both the height and the eruptive sequence. This emphasises the importance of studying these eruptions, utilising as many data sources as possible, in order to better understand the volcano and the character of its eruptions so as to be better prepared for future eruptive events (e.g. Joseph et al., 2022; Barclay et al., 2022).