Enhanced characterization of SO₂ plume height and column density using the second UV spectral band of TROPOMI

Volcanic emissions of sulfur dioxide (SO₂) affect the environment, climate, and society. Their detection and quantification rely extensively on remote sensing techniques, which are used to track SO₂ and monitor volcanic activity worldwide. In particular, nadir-viewing satellites measuring total SO₂ vertical column densities (VCDs) have provided valuable insights into volcanic emissions for decades. However, the determination of the SO₂ layer height (LH) is more challenging. In this study, we present an improved SO₂ LH (and VCD) retrieval algorithm, applicable to the second UV spectral band (BD2) of the TROPOspheric Monitoring Instrument (TROPOMI). This band exhibits a stronger SO₂ absorption than the third band (BD3) that is traditionally used for SO₂ retrievals. To assess the impact of various spectral, atmospheric, and observation conditions, we conducted sensitivity analyses from a set of synthetic spectra representative of TROPOMI measurements using the Look-Up Table COvariance-Based Retrieval Algorithm (LUT-COBRA). Our results demonstrate that BD2 retrievals result in more accurate estimates of the SO₂ heights and columns, particularly in the upper troposphere and lower stratosphere (UTLS), with LH errors reduced by at least a factor of 2. The algorithm was applied to real TROPOMI observations of volcanic eruptions and degassing episodes, and compared to BD3 retrievals. BD2 shows an improved sensitivity, with less noise, and a detection limit as low as 2 DU, surpassing the current operational TROPOMI SO₂ product by an order of magnitude. Furthermore, our plume height estimates align closely with independent measurements from the Infrared Atmospheric Sounding Interferometer (IASI) and Microwave Limb Sounder (MLS), confirming the reliability of the approach.

1 Introduction

Sulfur dioxide (SO₂) is an important trace gas emitted by volcanoes. Once released into the atmosphere, it undergoes oxidation, leading to the formation of sulfate aerosols that can significantly impact aviation (Schmidt et al., 2014; Kristiansen et al., 2024), human health (Carlsen et al., 2021; Heaviside et al., 2021), ecosystems (Tortini et al., 2017; Weiser et al., 2022), climate (Robock, 2000; Haywood and Boucher, 2000; Solomon et al., 2011; Santer et al., 2014; Marshall et al., 2025), and atmospheric chemistry (Stevenson et al., 2003; Lamotte et al., 2021). Among the volcanic ejecta (i.e., rock fragments, ash, water vapour, carbon dioxide, and trace gases), SO₂ is one of the most abundant. It is readily detectable using ultraviolet (UV) and thermal infrared (TIR) remote sensing techniques (Kern et al., 2013; Theys et al., 2013).

Ground-based instruments can be used to study volcanic emissions (Galle et al., 2010; Scollo et al., 2014), but their limited coverage restricts the analysis to small eruptions and leaves many volcanoes worldwide unmonitored. Space nadir-viewing sensors represent an excellent alternative, and have yielded key information on the total SO₂ vertical column density (VCD) for more than 40 years (Carn et al., 2016). However, the retrieval of the layer height (LH), defined as the altitude at which the plume concentration peaks, is a more recent advancement, although essential for several applications. Indeed, information on the SO₂ height can serve as a proxy for volcanic ash (Thomas and Prata, 2011), and tracking SO₂ clouds is therefore crucial for air traffic safety, as SO₂ and colocated aerosols can cause long-term damage to aircraft engines (Carn et al., 2009; Sears et al., 2013). Beyond operational concerns, the plume altitude is a key parameter in atmospheric and climate models, as it helps to further evaluate the transport, chemical transformation, and effects of volcanic gas emissions on air quality (Schmidt et al., 2015; Milford et al., 2023) and radiative forcing (Fuglestvedt et al., 2024; Zhuo et al., 2024; Axebrink et al., 2025). From a volcanological perspective, the knowledge of SO₂ LH also provides deeper insights into eruption dynamics, including eruption type, emission rate, intensity, and underlying processes (Mastin et al., 2009). Moreover, because satellite SO₂ retrievals are sensitive to the assumed vertical distribution of sulfur dioxide, reliably determining the plume height is crucial to quantify the estimates of volcanic SO₂ emissions (Yang et al., 2010).

In the thermal infrared, using respectively the Infrared Atmospheric Sounding Interferometer (IASI; Clerbaux et al., 2009) and the Cross-track Infrared Sounder (CrIS; Zavyalov et al., 2011), Carboni et al. (2012), Clarisse et al. (2014), and Hyman and Pavolonis (2020) have retrieved SO₂ heights with an excellent sensitivity in the upper troposphere and lower stratosphere (UTLS), even for SO₂ VCDs as low as 1 DU (Dobson Unit; 1 DU = 2.69×10¹⁶ molec. cm⁻²). Nonetheless, at altitudes below 5 km, the sensitivity to SO₂ considerably decreases due to a strong absorption by water vapour in the infrared and an overall reduced thermal contrast.

In comparison, SO₂ LH retrievals in the ultraviolet are less developed. The pioneering UV fitting algorithms, based on full radiative transfer simulations, have been applied to measurements from the Ozone Monitoring Instrument (OMI; Levelt et al., 2006) by Yang et al. (2010) and the Global Ozone Monitoring Experiment-2 (GOME-2; Munro et al., 2016) by Nowlan et al. (2011). Such retrieval schemes are time-consuming, and the reported results lack accuracy and precision. Another approach, the Full-Physics Inverse Learning Machine (FP_ILM) algorithm, has been developed and tested on GOME-2 by Efremenko et al. (2017), the TROPOspheric Monitoring Instrument (TROPOMI) onboard the Sentinel-5 Precursor (S-5P) platform (Veefkind et al., 2012) by Hedelt et al. (2019), and OMI by Fedkin et al. (2021). This method, based on machine learning schemes, is computationally efficient and is now implemented in the TROPOMI operational processing system. However, the FP_ILM algorithm is only sensitive to SO₂ height for VCDs larger than 20 DU, limiting its application to relatively large volcanic events.

Recently, the retrieval of the SO₂ slant column density (SCD) from TROPOMI has been greatly improved through the development of the COvariance-Based Retrieval Algorithm (COBRA; Theys et al., 2021), which notably results in a strong reduction of the retrieval noise. This advancement has naturally led to the development of an extension for the plume height retrievals by Theys et al. (2022), which combines COBRA with an iterative look-up table approach (LUT-COBRA) to reconstruct the SO₂ signal. Interestingly, this algorithm can derive SO₂ LHs for VCDs as low as about 5 DU.

Here, our main objective is to further develop LUT-COBRA to get a fast, accurate, precise, and more sensitive retrieval algorithm, especially for low SO₂ columns and SO₂ plumes in the UTLS. To achieve this, we exploit the second UV spectral band (BD2) of TROPOMI, that covers wavelengths below 310 nm, where SO₂ exhibits strong absorption features, and which offers better spectral performance than the third UV band (BD3) conventionally used for SO₂ retrievals (Hedelt et al., 2019; Theys et al., 2022). As a matter of fact, short UV wavelengths have been commonly used for ozone (O₃) profile retrievals (Veefkind et al., 2021), e.g., by Mettig et al. (2021), who exploited TROPOMI measurements between 270 and 329 nm and demonstrated high-quality retrievals. Since a single molecule of SO₂ absorbs nearly twice as strongly as a molecule of O₃ in the short UV (McPeters, 1993), this strategy has the potential to provide a better sensitivity than the BD3 retrievals of Theys et al. (2022). As a first step, we generated a set of SO₂ spectra under typical TROPOMI measurement conditions to assess the sensitivity of LUT-COBRA in BD2, and study the impact of different spectral, atmospheric, and observation conditions on the retrieval quality. From this synthetic framework, we designed an optimal look-up table and further refined our algorithm to analyze real TROPOMI BD2 measurements.

This paper is organized as follows. Section 2 presents the TROPOMI instrument and the methodology to improve SO₂ retrievals. Section 3 provides the theoretical basis of LUT-COBRA, followed by an overview of the sensitivity tests. In Sect. 4, we describe the developments made to our algorithm to process actual TROPOMI BD2 measurements and discuss the main results for various examples of volcanic emissions, comparing our retrievals with available correlative data, such as TROPOMI in BD3, IASI, and the Microwave Limb Sounder (MLS; Waters et al., 2006). In conclusion, Sect. 5 summarizes the key points of our study and outlines possible future steps.

2 TROPOMI

TROPOMI is a hyperspectral nadir-viewing sensor onboard the Sentinel-5 Precursor satellite, launched in 2017 as part of a European Space Agency (ESA) and Copernicus program. This spectrometer performs daily measurements of the solar radiation reflected by the Earth surface and backscattered by the atmosphere. It employs 8 spectral bands covering the ultraviolet, visible (VIS), near-infrared (NIR), and shortwave infrared (SWIR) spectral regions to retrieve information on tropospheric and stratospheric trace gas constituents, such as HCHO (De Smedt et al., 2018), CO (Borsdorff et al., 2018), NO₂ (van Geffen et al., 2020), CH₄ (Lorente et al., 2021), O₃ (Spurr et al., 2018), and SO₂ (Theys et al., 2017, 2019), as well as on aerosol (Torres et al., 2020) and cloud (Loyola et al., 2018) properties from local to global scales. The S-5P platform operates on a polar orbit at an altitude of approximately 800 km, crossing the equator at 13:30 LT (local time) and acquiring data over a 2600 km wide swath. In this work, we compare the UV spectral bands 2 and 3, for which the swath is divided into 448 and 450 across-track positions (or “rows”), respectively. This configuration ensures a high spatial resolution, with a ground pixel size at nadir of 3.5×5.5 km² (across-track × along-track). Additional details about the TROPOMI design and specifications can be found in Veefkind et al. (2012).

While SO₂ is traditionally retrieved between 310 and 340 nm (Hedelt et al., 2019; Theys et al., 2021, 2022), here we leverage wavelengths below 310 nm to benefit from the strong SO₂ differential absorption features (see Fig. 1a). For TROPOMI, this means exploiting the second band instead of the third one. BD2 spans wavelengths from 300 to 326 nm with a spectral sampling of 0.065 nm, whereas BD3 covers longer wavelengths, between 310 and 405 nm with a coarser sampling step of 0.2 nm. BD2 also has a finer spectral resolution of 0.5 nm, compared to the 0.55 nm of BD3, according to Veefkind et al. (2012). These characteristics suggest that BD2 provides more information per wavelength than BD3. However, a notable limitation of this second band, as highlighted in Fig. 1a, is the strong ozone absorption at shorter UV wavelengths. In addition, Fig. 1b shows that the signal-to-noise ratio also decreases below 310 nm, resulting in higher noise levels. Note that between 310 and 320 nm, the BD2 SNR resampled to the coarser BD3 spectral sampling appears higher than the native BD3 SNR.

Figure 1 (a) Absorption cross sections (ACS) of SO₂ (in orange) from Bogumil et al. (2003) and O₃ (in black) from Serdyuchenko et al. (2014), for a temperature (T) of 293 K. For illustration purposes, the data are convolved using a Gaussian function with a standard deviation of 0.5 nm. (b) Signal-to-noise ratio (SNR) estimated from TROPOMI Level-1 data for the second (in blue) and third (in red) UV spectral bands. The BD2 SNR scaled to the BD3 spectral sampling (in green) is also shown.

As a test case, we retrieved the SO₂ SCD from both TROPOMI bands using COBRA (Theys et al., 2021). The results are illustrated in Fig. 2a and b for the third and second bands, respectively. These measurements, taken on 18 January 2022, highlight the volcanic degassing of Popocatepetl and Wolf (see Table 4) in the western tropics. The SCDs from both bands are in excellent agreement, with a standard deviation of 0.37 DU for the difference between BD3 and BD2, considering only pixels where both exceed 0.5 DU. Moreover, we find the retrieval noise to be significantly reduced in BD2 compared to BD3, with SCD standard deviations of 0.13 and 0.30 DU, respectively, between 10–11° N, a SO₂-free region. This is promising and suggests better performance in the second band, despite the greater contribution of ozone (see Fig. 1a) and the lower SNR (see Fig. 1b) around 305 nm. The improvement is mainly attributed to the stronger SO₂ absorption, combined with an excellent spectral sampling in BD2, approximately three times better than in BD3. Encouraged by these results, we applied a similar strategy to retrieve the SO₂ plume height and vertical column density. The next section describes the retrieval approach and provides an overview of the synthetic sensitivity tests performed on the algorithm to evaluate its performance.

Figure 2 SO₂ slant column density retrieved from TROPOMI UV measurements performed on 18 January 2022, using BD3 (310–326 nm) in (a) and BD2 (305–320 nm) in (b). Black triangles indicate the locations of the Popocatepetl and Wolf volcanoes (see Table 4).

3 SO₂ retrievals and sensitivity analyses

3.1 LUT-COBRA description

Our TROPOMI BD2 retrieval algorithm relies on the Look-Up Table COvariance-Based Retrieval Algorithm (LUT-COBRA), previously developed by Theys et al. (2022) for BD3. Here, we review the theoretical foundations of this approach. Comprehensive discussions on the joint retrieval of SO₂ plume heights and column densities can be found in Yang et al. (2010).

Let y denote the total slant optical depth (SOD) at the top of the atmosphere (TOA), defined as the negative logarithmic ratio between the measured wavelength-calibrated radiance I and irradiance I₀, that is $y=-\log \left(\frac{I}{I_{\mathrm{0}}}\right)$. A general formulation is

$$\begin{array}{}\text{(1)}& y=y_{{\mathrm{SO}}_{\mathrm{2}}}+y_{\mathrm{bckg}}+\mathit{ϵ},\end{array}$$

where $y_{{\mathrm{SO}}_{\mathrm{2}}}$ describes the SO₂ slant optical depth, y_bckg the background contribution (e.g., from O₃ absorption), and ϵ the measurement error. The total SOD can be approximated by the following Taylor expansion,

$$\begin{array}{}\text{(2)}& y-y_{{\mathrm{SO}}_{\mathrm{2}},i}\approx \mathrm{\Delta }\mathrm{VCD}{\displaystyle \frac{\partial y_{{\mathrm{SO}}_{\mathrm{2}},i}}{\partial \mathrm{VCD}}}+\mathrm{\Delta }\mathrm{LH}{\displaystyle \frac{\partial y_{{\mathrm{SO}}_{\mathrm{2}},i}}{\partial \mathrm{LH}}}+y_{\mathrm{bckg}}+\mathit{ϵ},\end{array}$$

with $y_{{\mathrm{SO}}_{\mathrm{2}},i}=y_{{\mathrm{SO}}_{\mathrm{2}}}({\mathrm{VCD}}_i,{\mathrm{LH}}_i)$ the SOD at the linearization point i, [ΔVCD, ΔLH] the increments in column density and plume height, and $\left[\frac{\partial y_{{\mathrm{SO}}_{\mathrm{2}},i}}{\partial \mathrm{VCD}},\frac{\partial y_{{\mathrm{SO}}_{\mathrm{2}},i}}{\partial \mathrm{LH}}\right]$ the associated partial derivatives (Jacobians).

Equation (2) can be solved with LUT-COBRA, an iterative method that models the SO₂ spectrum through a look-up table containing precomputed SO₂ SODs and Jacobians for a wide range of scenarios. Basically, the algorithm starts by extracting these data for conditions close to the measurements, and taking into account a priori values for the LH and VCD. Each iteration then refines the estimates of SO₂ height and column density, which are subsequently used to update the optical depths and Jacobians. This process is repeated until convergence to the solution ${\widehat{x}}_i=[{\mathrm{VCD}}_i,{\mathrm{LH}}_i]$, given by

$$\begin{array}{}\text{(3)}& {\widehat{x}}_{i+\mathrm{1}}={\widehat{x}}_i+{\left({\mathbf{k}}_i^T{\mathbf{S}}^{-\mathrm{1}}{\mathbf{k}}_i\right)}^{-\mathrm{1}}{\mathbf{k}}_i^T{\mathbf{S}}^{-\mathrm{1}}\left(y-y_{{\mathrm{SO}}_{\mathrm{2}},i}-\stackrel{\mathrm{‾}}{y}\right),\end{array}$$

where k_i is a two-column matrix of SO₂ Jacobians, and all the contributions other than SO₂ are treated through a generalized error covariance matrix S and its associated mean spectrum $\stackrel{\mathrm{‾}}{y}$. Indeed, LUT-COBRA does not directly fit the background, but instead, it considers an ensemble of N measured spectra $y_l^{\prime }$ ($l\in [\mathrm{1},N]$) that are unaffected by SO₂. These spectra are characterized by a mean optical depth $\stackrel{\mathrm{‾}}{y}$, so that,

$$\begin{array}{}\text{(4)}& \mathbf{S}={\displaystyle \frac{\mathrm{1}}{N-\mathrm{1}}}\sum _{l=\mathrm{1}}^N\left(y_l^{\prime }-\stackrel{\mathrm{‾}}{y}\right){\left(y_l^{\prime }-\stackrel{\mathrm{‾}}{y}\right)}^T,\end{array}$$

is used to statistically represent the combined error term y_bckg+ϵ of Eq. (1). Moreover, the retrieval errors, calculated as the square root of the diagonal elements of the solution covariance matrix, i.e.,

$$\begin{array}{}\text{(5)}& {\widehat{x}}_{\mathrm{err},i}=\sqrt{{\left({\mathbf{k}}_i^T{\mathbf{S}}^{-\mathrm{1}}{\mathbf{k}}_i\right)}^{-\mathrm{1}}},\end{array}$$

can be used to describe the quality of the solutions found by the algorithm, as demonstrated below.

3.2 Simulation conditions and algorithm settings

To evaluate the performance of LUT-COBRA across different wavelength ranges (i.e., spectral bands), a set of synthetic spectra was generated using the Linearized Discrete Ordinate Radiative Transfer (LIDORT) scalar model (Spurr and Christi, 2019). Although LIDORT neglects polarization effects, these typically alter TOA radiances by less than 5 % (Lerot et al., 2014; Escribano et al., 2019), mainly in a spectrally smoothed way, suggesting only a minor influence on the retrievals. Simulations were performed under typical TROPOMI measurement conditions between 300 and 330 nm, with a spectral sampling of 0.065 nm for BD2 and 0.2 nm for BD3. The observation conditions were arbitrarily fixed to a solar zenith angle (SZA) of 10°, a viewing zenith angle (VZA) of 0°, a relative azimuth angle (RAA) of 0°, a surface height of 0 km, and a surface albedo of 5 %. It should be noted that clouds, aerosols, and rotational Raman scattering were not included in this study. While these factors may introduce uncertainties of a few percent in the radiative transfer (Yang et al., 2010; Fedkin et al., 2021), and up to 5 km bias on the SO₂ LH for highly absorbing aerosols collocated with SO₂, our aim was to maintain a simple and controlled synthetic framework. Furthermore, we used the US Standard (US STD) atmospheric profiles of ozone, temperature, pressure, and air density (Anderson, 1986). For sulfur dioxide, a Gaussian function with a standard deviation of 0.5 km was applied to various combinations of SO₂ heights and columns. Each profile was actually centered on a specific layer height, ranging from 1 to 25 km in 1 km increments, and from 25 to 45 km in 5 km increments, and was scaled to the following VCDs: 1, 2, 5, 10, 15, 20, 25, 30, 40, 50, 75, 100, 125, 175, 250, and 300 DU. For high SO₂ columns, a coarser sampling of VCD nodes was used to limit the total number of simulations. However, sensitivity tests indicate that interpolation errors in the SODs and Jacobians are negligible under these conditions, with differences of less than 0.5 % in the radiative transfer, even at shorter UV wavelengths. These profiles were interpolated on an altitude grid from 0 to 50 km, with steps of 0.5 km up to 28 km, then 1 km. The temperature-dependent absorption cross sections of ozone by Serdyuchenko et al. (2014) and sulfur dioxide by Bogumil et al. (2003) were used as input to the simulations, after convolution with a Gaussian function characterized by a standard deviation of 0.47 nm to approximate the TROPOMI Instrument Spectral Response Function (ISRF) in the UV. It is important to stress that the above description stands for our “baseline” set of synthetic spectra. In addition, to investigate the role of influencing quantities, this reference set is complemented by many other synthetic spectra corresponding to perturbed conditions in terms of atmospheric profiles (e.g., temperature, pressure and air density, O₃ and SO₂), absorption cross-section datasets, albedos, and surface heights. This will be further addressed in Sect. 3.3 and Appendix A.

With this configuration in mind, we processed the simulated spectra for each condition. The outputs comprise the SO₂-free radiance, $I_{{\mathrm{O}}_{\mathrm{3}}}$, which only includes O₃ absorption, and the combined radiance, $I_{{\mathrm{SO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{3}}}$, which includes absorption from both SO₂ and O₃. These spectra were then used to compute the SO₂ SOD as,

$$\begin{array}{}\text{(6)}& y_{{\mathrm{SO}}_{\mathrm{2}}}=-\ln \left({\displaystyle \frac{I_{{\mathrm{SO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{3}}}}{I_{{\mathrm{O}}_{\mathrm{3}}}}}\right).\end{array}$$

From the different combinations of plume heights and column densities, the LH and VCD Jacobians in Eq. (2) were derived by finite difference. As an example, the dependence of the SO₂ SOD on the SO₂ altitude and column density is presented on Fig. 3 for various SO₂ profiles. Overall, the results are in line with the findings of Yang et al. (2010) and Nowlan et al. (2011). In particular, the sensitivity to the plume height increases with higher SO₂ abundances, as the signal becomes more pronounced. The LH Jacobians are also more intense for lower-altitude plumes. This can be attributed to the strong vertical gradient of the air mass factors in the lower atmosphere (see Fig. B1 in Appendix B), where the measurement sensitivity increases rapidly with the altitude. As a result, small changes in plume height induce significant variations in the SO₂ signal, leading to larger LH Jacobians in that region. In contrast, the VCD Jacobians tend to be larger for higher-altitude plumes, where the AMFs, and thus the sensitivity to the SO₂ column, are higher. However, their amplitude decreases with an increasing SO₂ column, due to the non-linear nature of SO₂ absorption, which leads to a saturation of the signal. Overall, for both Jacobians, the results indicate an improved performance below 310 nm, consistent with the discussions of Sect. 2. This is especially evident for high plume heights (above 15 km), as noted by Yang et al. (2010).

Figure 3 SO₂ LH (a) and VCD (b) Jacobians for different column densities and plume heights.

In a next step, the calculated SO₂ SOD and Jacobians were incorporated into a small look-up table. Before the retrievals, noise levels, representative of BD2 or BD3, could be added to the synthetic spectra to be fitted, as

$$\begin{array}{}\text{(7)}& y^{*}=y+{\displaystyle \frac{y_n}{\mathrm{SNR}}},\end{array}$$

where y_n is a spectrum of random real numbers sampled from a univariate Gaussian distribution with mean 0 and variance 1. In case such noise was introduced, we considered a diagonal covariance matrix, with ${\mathbf{S}}_{ii}=\left(\frac{\mathrm{1}}{{\mathrm{SNR}}_i^{\mathrm{2}}}\right)$ for the ith wavelength, consistent with the data shown in Fig. 1b. In the absence of noise, however, we used a unit covariance matrix for the retrievals. Note that for the synthetic retrievals, we assumed the background contribution to be perfectly known, which means that SO₂ spectra, instead of total spectra, were fitted (see Eq. A1). In practice, the algorithm is considered to have converged when the variation in plume height is less than 0.25 km and the change in VCD is below 5 % between successive iterations. To ensure realistic solutions and computational efficiency, the number of iterations is limited to 10. Results of these synthetic tests are presented below.

3.3 Synthetic retrievals and error analysis

In this section, we first evaluate the performance of LUT-COBRA for both UV bands of TROPOMI, accounting for their respective spectral characteristics in terms of sampling, resolution, and SNR (see Sect. 2). Then, we assess the systematic errors on the retrievals by perturbing input data of the radiative transfer. These synthetic analyses were carried out to understand how the algorithm responds to forward model and measurement uncertainties. Spectra in BD2 and BD3 were computed from different input SO₂ VCDs and LHs, detailed in Table 1.

Table 1 Overview of the SO₂ VCDs and LHs as well as the corresponding a priori values provided as input to LUT-COBRA.

The a priori inputs were arbitrarily selected, but our tests have revealed that this choice has a negligible impact on the retrieval quality. Figure 4 summarizes the resulting biases for the LH and VCD, in both bands, averaged over 100 noisy spectra. It can be seen that the LH and VCD biases are globally reduced in BD2, especially for lower SO₂ heights and column densities. They are typically within the convergence criteria (0.25 km for the LH and 5 % for the VCD) of the algorithm. We also observe that the standard deviations of the retrievals are, as expected, greatly reduced compared to BD3. Figure 4 thus demonstrates a better accuracy and indicates that, in principle, there is no fundamental limitation in retrieving the SO₂ plume height and column density from the TROPOMI BD2.

Figure 4 (a) Absolute difference between the retrieved and expected SO₂ LHs as a function of the input parameters. (b) Relative difference between the retrieved and expected SO₂ VCDs, expressed as a percentage, and shown as a function of the input parameters. In both cases, retrievals are performed using synthetic spectra with added noise, in BD2 (305–320 nm, blue) and BD3 (310–326 nm, red). Circles and squares indicate the corresponding SO₂ VCDs or LHs. Error bars represent the standard deviation associated with each retrieval. Note that the data are slightly offset along the x axis for clarity.

A key question is whether BD2 can enhance the retrieval precision compared to BD3. To investigate this, we estimated the SO₂ LH retrieval error, as defined in Eq. (5), for both bands. The results, given in Fig. 5, correspond to a prescribed SO₂ column of 2 DU, with plume heights going from the lower troposphere to the upper stratosphere. For BD2, retrieval errors are shown for different fitting ranges, varying the lower wavelength while keeping the upper limit fixed at 320 nm. In contrast, BD3 is evaluated using a single spectral range (310–326 nm), similar to Theys et al. (2022). From Fig. 5, it can be concluded that the use of BD2 results in a remarkable improvement in precision compared to BD3. When the same lower wavelength limit (310 nm) is applied, BD2 already shows slightly lower LH errors in the UTLS relative to BD3, confirming the positive impact of the finer spectral sampling. Figure 5 further demonstrates that extending the BD2 spectral range to shorter wavelengths significantly enhances the retrieval performance, particularly for high-altitude plumes (i.e., above 25 km). For a lower limit of 305 nm, the retrieval error is reduced by approximately a factor of 2. This is consistent with the behaviour observed for the air mass factors (see Fig. B1). At high altitudes, AMFs are greater for wavelengths larger than 310 nm, but they vary less, whereas the information on the plume height primarily comes from this weak altitude dependence and, indirectly, from the temperature dependence of the SO₂ absorption cross sections. However, for shorter wavelengths (or higher SO₂ VCDs), the height dependence of the sensitivity functions is increasing. The algorithm is therefore able to extract information more easily from the spectrum, thereby reducing retrieval errors. On the other hand, including shorter wavelengths in the fit may introduce significant systematic errors in the retrieval. In the remainder of this section, we evaluate the potential sources of uncertainty (e.g., on ozone) and estimate the associated errors.

Figure 5 BD2 (in colours) retrieval errors of SO₂ LH for different lower limits of the fitting spectral range, with the upper limit fixed at 320 nm. Results are compared with those from BD3 (in black) using a 310–326 nm fitting window. All errors are shown for a SO₂ VCD of 2 DU.

In this context, the sensitivity of LUT-COBRA in BD2 was tested by inverting noise-free spectra simulated with a single deviation from the baseline conditions of our look-up table (see Sect. 3.2). These perturbations were applied to atmospheric parameters, including temperature, air density (pressure), O₃ and SO₂ profiles, as well as to the SO₂ absorption cross sections, albedo, and surface height. In each case, variations of different magnitudes were applied, as indicated in Table 2.

Bogumil et al. (2003)(Birk and Wagner, 2018)

Table 2 Overview of all the sensitivity tests of LUT-COBRA. Analyses performed without noise in TROPOMI BD2, over the following fitting spectral range (sampling): 305–320 nm (0.065 nm). For the baseline conditions, US STD denotes the US Standard data from Anderson (1986). The profile notation stands for North (NH) or South (SH) Hemisphere, at mid- (±30–60°) or polar (±60–90°) latitudes during the Winter–Spring (WS) or Summer–Fall (SF) seasons. These data are drawn from the climatology of Lamsal et al. (2004).

For the temperature, we used the climatology of Lamsal et al. (2004), selecting profiles representative of different seasons in the northern hemisphere. Note that we also analyzed a case in the tropics (not reported in Table 2). These distributions, originally provided for different ozone column densities, were interpolated to a fixed VCD of 345.7 DU (matching the standard value) so that only the shape of the profile differs from our baseline conditions. The effect of air density and pressure, in turn, was evaluated by adjusting the air concentration in each layer of the US Standard atmosphere (Anderson, 1986), as defined in our LUT. Regarding the ozone contribution, we examined two sources of uncertainty: the shape of its profile and the corresponding VCD. Given the strong O₃ absorption in BD2, a broader range of scenarios was explored. For the vertical distribution, we considered four profiles (Lamsal et al., 2004), spanning latitudes from 30° S to 90° N (tropics and southern hemisphere profiles are not listed in Table 2) and all seasons. As with the temperature analysis, these profiles were scaled to the standard O₃ VCD to test the effect of a change in profile shape independently. Afterwards, to assess the impact of the VCD, we tested values ranging from 310 to 380 DU (only two cases are provided in Table 2), while keeping the same standard profile shape. For SO₂, we investigated the impact of different vertical distributions by fitting spectra corresponding to rectangular, Lorentzian, and Gaussian profiles of varying standard deviations (σ_G), for the LHs and VCDs given in Table 1. The influence of the SO₂ absorption cross sections was also studied, by comparing the SEOM-IAS (Birk and Wagner, 2018) dataset with that from Bogumil et al. (2003). Finally, we looked at uncertainties related to observation conditions, in particular the albedo and surface height.

All the above sensitivity analyses were performed within the 305–320 nm spectral range, selected as an optimal balance, minimizing LH retrieval errors (see Fig. 5) and reducing systematic uncertainties related to the ozone absorption. Figure 6 summarizes the mean systematic errors on the retrieved LHs and VCDs, assuming extreme and more realistic (see Table 2) perturbations in each case, for a relatively low column density (5 DU) and UTLS altitude (13.5 km).

Figure 6 BD2 LUT-COBRA retrieved biases for the LH (a) and VCD (b), considering a SO₂ height of 13.5 km and column density of 5 DU as solutions. In particular, the results are shown for perturbations on the temperature (T prof.) and air density/pressure (Air prof.) profiles, O₃ vertical distribution (O₃ prof.) and its VCD, SO₂ profile type (SO₂ prof.) and its standard deviation (SO₂ std dev.), SO₂ absorption cross section (SO₂ ACS), albedo and surface height (z_surf).

Overall, the total systematic errors are estimated by combining the individual effects of each source of uncertainty in quadrature, yielding 1.91 km for the LH and 2.58 % for the VCD. Since some contributions are partially correlated (e.g., temperature and air density profiles, or the O₃ profile and its VCD), these values should be interpreted as indicative and approximate estimates of the total error. The temperature, pressure and air density, as well as the sulfur dioxide profiles appear to have a negligible impact on the retrieval accuracy. Similarly, the albedo and surface height have little effect on the quality of the results, especially since these two observation conditions are generally well constrained in practice. In contrast, the choice of SO₂ absorption cross-section dataset has a more important influence on the retrieved solutions, particularly for the plume height. A comparison of the two datasets (see Fig. C1) shows that their temperature dependence is similar, but the dataset from Birk and Wagner (2018) assumes a stronger SO₂ absorption than that of Bogumil et al. (2003), especially below 310 nm. As a result, SO₂ air mass factors, as defined by Eq. (B1), decrease at the wavelengths with the highest absorption. Since the SO₂ SOD depends on the product of the AMFs and ACS, the fitted signal is “smoothed”, ultimately leading to systematically higher SO₂ LH. Importantly, we observe from Fig. 6, that the ozone profile, particularly its vertical shape, plays a key role on the retrievals, as expected given its strong absorption at shorter UV wavelengths. Nevertheless, the obtained LH and VCD biases remain relatively modest. Note that we conducted further tests using alternative baseline conditions (i.e., with various look-up tables) for the temperature, O₃ profile, SO₂ ACS, and albedo. These complementary analyses are detailed in Table A1 and lead to consistent conclusions. The results are therefore promising, but they highlight the need to mitigate the influence of ozone. Our strategy to better account for its strong contribution is described in Sect. 4.1.

4 Application of LUT-COBRA to TROPOMI measurements

4.1 Algorithm developments

Following the synthetic analyses, we developed a more comprehensive look-up table and refined LUT-COBRA for application to actual TROPOMI BD2 observations. As explained in Sect. 3.1, a LUT is a set of pre-calculated SO₂ spectra for various conditions, which allows efficient retrievals by avoiding the computational cost of online radiative transfer simulations. In this updated LUT, we considered a broader range of scenarios, with particular emphasis on the role of ozone, which emerged as a key factor affecting the retrieval quality (see Sect. 3.3). As a first step, we performed a statistical analysis of atmospheric profiles, including temperature, air density, and ozone. Our objective was to build a compact set of representative O₃ profiles parametrized based on the total ozone column only. This was motivated by the need to limit the LUT size while adequately representing the ozone profile variability in the retrievals. The idea of using the O₃ VCD to inform about the shape of the O₃ profile is not new (Labow et al., 2015 and references therein). It has been extensively used for satellite total ozone retrievals (e.g., by Lerot et al., 2014) and, in our application, has the potential to reduce the uncertainties around the O₃ profile (see Fig. 6), given that the total ozone column is retrieved from TROPOMI with an excellent accuracy (Garane et al., 2019). To this aim, we used the climatology of Lamsal et al. (2004) that provides vertical distributions of temperature and ozone density for several O₃ VCDs, latitudes, and seasons. From these data, we derived the following mean atmospheric profiles,

$$\begin{array}{}\text{(8)}& {\overline{P}}_i={\displaystyle \frac{\sum _{j,k}{\mathit{\alpha }}_{i,j,k}{P}_{i,j,k}}{\sum _{j,k}{\mathit{\alpha }}_{i,j,k}}},\end{array}$$

with $P_{i,j,k}$ the vertical distribution (of temperature, air density or ozone) for a VCD i, latitude j and season k, and ${\mathit{\alpha }}_{i,j,k}$ a corresponding weighting factor which is calculated from a frequency distribution based on the OMI Level-3 ozone product (Bhartia, 2012). As an example, Fig. 7a shows the mean ozone density profile for 385 DU, resulting from different contributions, with their respective weighting factors. This highlights the regions and seasons where such a high concentration of ozone is most likely to be found. Similarly, Fig. 7b presents the calculated mean ozone distributions for various column densities. They thus correspond to the most probable profiles for such O₃ VCDs, and it can be seen that the O₃ peak shifts towards lower altitudes as its concentration increases, providing valuable information about the profile shape that helps reduce the associated uncertainties in the SO₂ retrievals.

Figure 7 (a) O₃ density profiles for various latitudes and seasons (in colours), considering an O₃ VCD of 385 DU. The profile notation stands for North (NH) or South (SH) Hemisphere, at mid- (±30–60°) or polar (±60–90°) latitudes during the Winter–Spring (WS) or Summer–Fall (SF) seasons. The weighting factor is indicated in parentheses for each distribution. The mean profile (in black) is also shown. (b) Mean O₃ profiles generated for different O₃ VCDs.

Accordingly, we calculated mean O₃ profiles for different O₃ column densities, summarized in Table 3. For consistency, mean air density and temperature profiles were computed using the same approach. These mean profiles were then used as input for radiative transfer simulations, covering a wide range of conditions (see Table 3). In total, approximately 60×10⁶ total radiance spectra and 180×10³ O₃ radiance spectra were produced using a high-performance computer (HPC). Note that we considered the same absorption cross-section datasets as previously (see Sect. 3.2) for ozone (Serdyuchenko et al., 2014) and sulfur dioxide (Bogumil et al., 2003). However, these spectral data were pre-convolved here with a 0.05 nm box-car function, sampled at 0.065 nm and interpolated to the mean temperatures. Afterwards, the simulated radiances were exploited to calculate the SO₂ optical depths and Jacobians. These data were convolved with the TROPOMI IRSF (varying as a function of the detector row) and corrected for the solar-I₀ effect (Aliwell et al., 2022), before being incorporated into a comprehensive look-up table to analyze TROPOMI BD2 measurements.

Table 3 Physical parameters that define the SO₂ slant optical depth look-up table.

The retrieval of SO₂ LH and VCD in BD2 is fundamentally the same as in BD3. The core of the algorithm is explained in more detail by Theys et al. (2022) and here we only describe the new developments applied to LUT-COBRA.

First, one needs to determine the parameters related to the reflection (i.e., albedo and surface height), geometry (i.e., SZA, VZA, RAA), and O₃ amount corresponding to a measured spectrum. These pieces of information are essential for the retrievals. Indeed, during the process, a sub-LUT depending only on the SO₂ LH and VCD, and representative of the measurements, is extracted by linear interpolation. The observation angles can usually be derived from the TROPOMI Level-1 files. For both the albedo and surface height, we assumed a Lambertian equivalent reflector (LER), similar to Theys et al. (2022), and characterized by effective parameters. This model provides a simplified representation of the complex mechanisms of reflection by the surface, clouds, and aerosols. In BD3, the effective height is computed from the surface and cloud altitudes weighted by the cloud fraction, based on data from the TROPOMI Level-2 product. The effective albedo, on the other hand, is calculated from the reflectance at 340 nm, a spectral region free from any absorption by sulfur dioxide or ozone to ensure a correct computation. Regarding the O₃ input, the corresponding operational offline TROPOMI Level-2 product (Spurr et al., 2018) is utilized. In addition, as outlined in Sect. 3.1, a priori values on the sulfur dioxide plume height and column density have to be provided to LUT-COBRA. The SO₂ LH a priori is set to 7 km, or LER height + 2 km if the LER height is higher than 5 km, and the prior VCD is defined from the SOD and AMF found in the operational TROPOMI SO₂ column product, assuming an altitude of 7 km. However, all these input data are not available in BD2 and must therefore be read from BD3. It is important to remind that the two bands do not share the same number of across-track positions (as explained in Sect. 2), and there is also a shift in the flight direction (along-track positions). We performed this BD3-to-BD2 transition by directly interpolating the BD3 data into BD2 using the coordinates (latitudes and longitudes) provided in TROPOMI Level-1 files. This approach proved to be a good approximation, while simplifying our script and optimizing the computation time.

Once all BD2-interpolated data are obtained, the covariance matrix needed for the retrieval (Eq. 4) is calculated for each row. A set of SO₂-free radiance spectra measured arbitrary for 300 pixels along the flight direction is selected for a given orbit. The SO₂ VCDs and corresponding retrieval errors from the TROPOMI Level-2 product are used to filter the pixels, discarding those with SO₂ VCDs > 2.5 × VCD errors. Note that to ensure a robust construction of the covariance matrix, at least 100 spectra must be included. However, we observed that the covariance matrix, as described in Eq. (3), sometimes becomes ill-conditioned in BD2, resulting in unphysical values in the inverse matrix for some wavelengths, and thereby compromising the retrievals. To address this issue, we performed an eigen-decomposition of the covariance matrix, and refined the inverse of the covariance matrix by discarding eigenvalues of S that were too small, as

$$\begin{array}{}\text{(9)}& \mathbf{S}=\sum _{i=\mathrm{1}}^n{\mathit{\lambda }}_i{v}_i{v}_i^T\to {\mathbf{S}}^{-\mathrm{1}}=\sum _{i=\mathrm{1}}^k{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_i}}{v}_i{v}_i^T,\end{array}$$

where λ_i and v_i corresponds to the eigenvalue and eigenvector of S, n is the number of wavelengths, k denotes the threshold beyond which S⁻¹ presents erroneous values. This correction, described in more detail by Clarisse et al. (2023), ensures that the inverse matrix is calculated correctly and avoids inconsistent retrieval errors. Tests were performed to identify the threshold minimizing the errors, and a value of 10⁻⁷ on the eigenvalues was found, corresponding approximately to the precision of single-precision floating-point residuals. The procedure was thus applied systematically to all spectra, with typically less than 30 % of eigenvalues discarded to stabilize the inversion, depending on the conditioning of the covariance matrix.

The data were subsequently provided to LUT-COBRA, initiating the fitting process. Note that our analysis is restricted to measurements with SZAs below 65°, in order to avoid conditions of strong ozone absorption that could complicate the retrievals. It should also be mentioned that the algorithm may return values outside the predefined grids (see Table 3) for some iterations. In such cases, the SO₂ height is set to surface height + 1 km or grid maximum − 1 km for the next iteration, depending on whether it falls below or above the limits. Similarly, if the derived SO₂ column density lies outside the VCD grid, the algorithm reverts to the a priori value for the next iteration.

4.2 Results and comparisons with satellite observations

In this section, we analyze the SO₂ LHs and VCDs retrieved from various TROPOMI BD2 observations of volcanic eruptions, as summarized in Table 4. These results are compared with BD3 estimates from LUT-COBRA (Theys et al., 2022) and the operational product (Hedelt et al., 2019), as well as with SO₂ measurements from IASI (Clarisse et al., 2014) and MLS (Pumphrey et al., 2015).

Table 4 Overview of the volcanic emissions analyzed with our TROPOMI BD2 SO₂ plume height retrieval algorithm.

4.2.1 Raikoke (June 2019)

As a first case, we examined the SO₂ plume a few days after the eruption of the Raikoke volcano on 21 June 2019. This explosive event, which was the first activity of the volcano in nearly 95 years, lasted about 24 h, injecting a substantial amount of sulfur dioxide, over 1.5 Tg, up to altitudes of 11 km and beyond (Cai et al., 2022; Vernier et al., 2024). The plume rapidly dispersed across the northern hemisphere, attracting a large scientific interest due to its impact on the environment and atmospheric composition (Grishin et al., 2021; de Leeuw et al., 2021). Such an eruption is particularly well-suited to assess the reliability of our approach, as the large SO₂ abundances facilitate the determination of the plume height.

We first tested different spectral fitting windows and found that the 305–326 nm interval was optimal to capture the strong SO₂ absorption in the short UV and minimize the O₃ (and noise) contribution, thereby enhancing the measurement sensitivity. BD3 LUT-COBRA has a reported detection limit of 5 DU, and the current S-5P operational product is much less sensitive below 20 DU, whereas our algorithm can detect SO₂ plumes down to 2 DU. Owing to this improvement, more detailed spatial and vertical information on the SO₂ plume can be obtained. This wavelength range is therefore adopted for all subsequent cases. Moreover, for the following analyses, only pixels with retrieved LH errors below 2.5 km are considered.

Figure 8 presents the retrieved SO₂ heights and column densities from BD2 and BD3 on 23 June 2019. We can see that the LHs and VCDs of both bands agree very well. The LHs in BD2 seem slightly lower than in BD3 but the plume filaments are better captured in BD2 data, indicating a better sensitivity. To further assess the consistency of our algorithm, we calculated the number of retrieved LHs and VCDs from both bands within bins of 1 km × 1 km and 10 DU × 10 DU, respectively, and performed a linear regression. The results, shown in Fig. 9, demonstrate a good agreement between the two products, which both identify the main characteristics of the SO₂ plume, and support the reliability of LUT-CORA across the two spectral bands.

Figure 8 SO₂ retrievals from TROPOMI UV measurements on 23 June 2019 for the Raikoke eruption (see Table 4). Panels (a) and (b) show the LHs retrieved using BD3 (310–326 nm) and BD2 (305–326 nm), respectively. Panels (c) and (d) display the corresponding VCDs. The volcano location is marked by a black triangle.

Figure 9 SO₂ LHs per 1 km × 1 km bin (a) and VCDs per 10 DU × 10 DU bin (b) retrieved from TROPOMI BD3 (310–326 nm) and BD2 (305–326 nm) measurements on 23 June 2019 for the Raikoke eruption (see Table 4). Each plot includes the 1:1 line (y=x), the correlation coefficient (r), and the linear fit ($y=ax+b$).

4.2.2 Etna (August 2024)

To investigate whether BD2 can enhance the retrieval performance in cases of low SO₂ VCDs in the UTLS, we focused on Mount Etna, Sicily, one of the most active volcanoes in Europe. During the summer of 2024, the volcano experienced a few minor eruptions, releasing SO₂ amounts generally below 10 DU, at altitudes below 10 km. One of the most intense emissions were recorded on 15 August (Sennert, 2024). For that day, the SO₂ columns were relatively small, ranging from 6 to 40 DU, and the plume peak was at about 8 km, making it a relevant case study to assess the retrieval sensitivity under challenging conditions.

Figure 10 presents SO₂ maps derived from TROPOMI BD3 observations using LUT-COBRA and FP_ILM. As recommended by Hedelt et al. (2023), a quality assurance threshold of 0.5 was applied to the operational product, while the suggested LH validity flag, designed for pixels with high SO₂ columns (typically above 15 DU) and retrieval errors exceeding 2 km, was omitted to preserve enough pixels and allow a meaningful comparison with LUT-COBRA. Overall, the two algorithms agree well, displaying similar patterns as well as consistent layer heights and column densities. Also shown in Fig. 10 are the results from BD2, compared to IASI-C measurements (version 4.1) for the descending (AM) orbit. We observe that BD2 values closely match BD3 results for common pixels while providing additional insights, particularly in the region between 17.5–22° E and 37–41° N, where plume features absent in BD3 are detected by BD2. This can be directly attributed to the better detection limit of our algorithm (i.e., VCD as low as 2 DU). In this region, the scatter in BD2 LH is quite strong due to the low SO₂ column amount, reflecting a reduced sensitivity. However, the mean plume height and column remain well constrained over the region, with BD2 estimates of 10.75 km and 3.53 DU, respectively, while IASI-C indicates a mean peak height of 10.45 km and column density of 2.32 DU. Our algorithm is therefore consistent with TIR observations. These results are further supported by retrieval errors that remain below 1.5 km for the layer height and 1 DU for the column density in BD2, which is in line with the synthetic tests presented in Sect. 3.3 and confirms the robustness of our method.

Figure 10 SO₂ retrievals from TROPOMI and IASI-C measurements on 15 August 2024 for the Etna eruption (see Table 4). Panels (a) and (b) show the LHs obtained from TROPOMI BD3 using LUT-COBRA and the operational product, respectively. Panel (c) presents the BD2 results, and panel (d) the observations from the descending (AM) orbit of IASI-C. The corresponding SO₂ VCDs are displayed in panels (e)–(h). The volcano location is marked by a black triangle.

4.2.3 Ruang (April–May 2024)

Then, to further illustrate the added value of BD2, we applied our algorithm to a case involving dispersed plumes. More specifically, we analyzed the volcanic eruptions of Ruang, Indonesia, on 17 and 30 April 2024. In the days following these events, the SO₂ plume spread across the tropical region was observed reaching altitudes of about 20 km, with column densities up to 20 DU (Dodangodage et al., 2025). Under these conditions, the BD3 algorithm struggles to retrieve the SO₂ plume height, likely due to the reduced SO₂ signal. In contrast, BD2 LUT-COBRA, with its improved sensitivity, offers a clearer depiction of the plume dispersion and height.

Figure 11 shows the SO₂ LHs retrieved from TROPOMI BD3 and BD2 measurements, along with the observations from IASI-B, on 19 April and 1 May 2024. For the first eruptive event, the plume is more dispersed, which results in significant noise (i.e., pixels with unrealistically low and high LH values) in BD3. Conversely, the plume from the second eruption is more compact, and although BD3 LUT-COBRA performs slightly better, there is still noticeable noise. BD2, however, provides a greater number of valid pixels with height information and significantly reduces the surrounding noise for both dates. Furthermore, it can be seen that SO₂ heights in BD2 are slightly lower than in BD3 for common pixels, by about 1–2 km. Interestingly, the BD2 results are largely consistent with the observations from IASI-B, which measured peak heights of approximately 17 km on 19 April and 1 May. The comparison with IASI is particularly relevant, as thermal infrared measurements are highly sensitive to plume heights in the mid- to upper troposphere and lower stratosphere. This consistency between our BD2 algorithm and IASI retrievals suggests that our method better represents the overall structure of the plume, especially at higher altitudes, compared to BD3. This is expected, given that the air mass factors are almost constant at such altitude for longer wavelengths, as already explained in Sect. 3.3 and shown in Fig. B1.

Figure 11 SO₂ LHs retrieved from TROPOMI and IASI-B measurements for the 2024 Ruang eruptions (see Table 4). Panels (a) and (b) show the results from TROPOMI on 19 April using BD3 (310–326 nm) and BD2 (305–326 nm), respectively. Panels (c) and (d) present the corresponding TROPOMI LHs on 1 May. Panels (e) and (f) display the LHs retrieved from IASI-B on 19 April during the descending (AM) and ascending (PM) orbits, respectively. Panels (g) and (h) show the corresponding IASI-B LHs on 1 May. The volcano location is marked by a black triangle.

Moreover, we derived the SO₂ mass distribution (i.e., SO₂ columns weighted by the pixel area at each altitude) from IASI and TROPOMI measurements, including both the current operational product and LUT-COBRA results. This is illustrated in Fig. 12. Here, the retrieval method from Hedelt et al. (2019) exhibits significant discrepancies, with plume height values substantially underestimated compared to LUT-COBRA and IASI. As previously mentioned, IASI detects a peak near the tropopause, a feature also captured by our algorithm, whereas BD3 LUT-COBRA estimates that a larger portion of the plume resides at higher altitudes in the stratosphere. We also note that the BD3 results capture a non-negligible SO₂ layer below 2 km but this is unlikely to be true and results from low retrieval sensitivity for the corresponding pixels. Regarding the total SO₂ mass (see Fig. 12), BD3 LUT-COBRA and IASI retrievals provide comparable results, while our BD2 algorithm derives a higher value. This can be explained by the greater number of detected SO₂ pixels in BD2 compared to BD3, as a result of the improved sensitivity of our approach. This leads the algorithm to infer a higher total mass. The discrepancy with IASI mainly arises from the instrument's partial coverage, suggesting that the plume may not be fully captured because of gaps between successive orbits. It is worth noting that when the total SO₂ masses are recalculated, neglecting the SO₂ below approximately 5 km (not sounded by IASI or presumably erroneously retrieved by TROPOMI), we obtain 166.28 kt on 19 April, and 139.86 kt on 1 May, aligning more closely with the IASI estimates.

Figure 12 SO₂ vertical mass profiles on 19 April (a) and 1 May 2024 (b) for the Ruang eruptions (see Table 4), derived from IASI-B measurements, as well as from TROPOMI observations using BD3 and BD2 LUT-COBRA, and the operational product. The total mass corresponding to each dataset is indicated in parentheses.

To further evaluate our retrievals, we analyzed SO₂ measurements from the MLS instrument over the same period. Specifically, we used version 5.1 data and applied the recommended selection criteria (Livesey et al., 2022), along with a minimum threshold on the SO₂ mixing ratio of 25 (unitless). The SO₂ density profiles derived from the ascending (daytime) and descending (nighttime) orbits are shown in Fig. 13, comparing data with a quality flag (q) greater than 0.95 (as recommended by Livesey et al., 2022) and data selected with a less restrictive threshold of 0.6. Note that we used the tropical atmospheric profiles of Lamsal et al. (2004) to convert the results from pressure to altitude. Despite its limited vertical resolution for SO₂ (of approximately 3 km), the MLS observations are globally consistent with both our retrievals and those from IASI (see Fig. 12). On 19 April, MLS indeed detects a LH near 17 km during the night, regardless of the quality flag, and a comparable height during the day when exploiting the lower-quality data. Similarly, on 1 May, a peak is observed around 15 km in all daytime measurements and in nighttime data when using the lower threshold. Although less accurate, these results are nonetheless convincing and demonstrate that our algorithm can derive SO₂ plume heights in good agreement with independent satellite observations.

Figure 13 SO₂ vertical density profiles retrieved from MLS measurements on 19 April (a) and 1 May 2024 (b) for the Ruang eruptions (see Table 4). For each day, comparisons are made between daytime and nighttime observations, for different values of the quality flag (q) used for data selection.

4.2.4 Ulawun (June 2019)

Finally, as an additional demonstration case, we also studied the volcanic eruption of Ulawun on 27 June 2029, which exhibited characteristics similar to the Ruang eruptions in terms of SO₂ column densities and altitudes reached. A comparison of the SO₂ LHs retrieved from TROPOMI in both bands is shown on Fig. 14 and reveals similar results. Our BD2 algorithm provides more detailed vertical information on the plume compared to BD3, showing a better agreement with IASI. The retrieval errors are also significantly lower in BD2, further supporting the enhanced precision and reliability of our approach.

Figure 14 SO₂ retrievals from TROPOMI and IASI-B measurements on 27 June 2019 for the Ulawun eruption (see Table 4). Panels (a) and (b) show the LHs retrieved from BD3 (310–326 nm) and BD2 (305–326 nm), respectively, while panels (c) and (d) display the LHs obtained from IASI-B during the descending (AM) and ascending (PM) orbits, respectively, on the same dates. Panels (e) and (f) show the LH errors of BD3 and BD2. The volcano location is indicated by black triangles.

In Fig. 15a, we show the SO₂ mass profiles derived from TROPOMI and IASI observations, while Fig. 15b presents the corresponding distributions from MLS daytime and nighttime measurements for the Ulawun eruption. The SO₂ vertical distribution retrieved from BD2 is in better agreement with IASI (see Fig. 14) than BD3, whether using LUT-COBRA or the S-5P operational product, with an estimated plume height around 17 km. The differences in total SO₂ mass between LUT-COBRA and FP_ILM are mainly due to the limited sensitivity of the operational product at such altitudes and low concentrations, resulting in only a few SO₂ pixels being detected. The explanation for the discrepancies between TROPOMI and IASI is similar to that discussed for the Ruang case. Despite these differences, it can be seen that the SO₂ profile retrieved from BD2 also aligns with the MLS profiles, showing a peak altitude near 15 km, regardless of the observation time or the applied quality flag.

Figure 15 (a) SO₂ vertical mass profiles on 27 June 2019 for the Ulawun eruption (see Table 4), retrieved from IASI-B measurements, as well as from TROPOMI observations using BD3 and BD2 LUT-COBRA, and the operational product. The total mass corresponding to each dataset is indicated in parentheses. (b) SO₂ vertical density profiles derived from MLS daytime and nighttime observations, using different quality flags (q) for data selection, for the same event.

5 Conclusions

In this study, we significantly improved the sensitivity of TROPOMI SO₂ layer height and vertical column density retrievals. Originally developed for BD3 by Theys et al. (2022), we extended LUT-COBRA to BD2. This adaptation leverages the stronger SO₂ absorption at shorter UV wavelengths and the superior spectral performance of TROPOMI in BD2 (mostly in terms of spectral sampling). We first demonstrated the benefit of BD2 by retrieving SO₂ slant columns using the classical COBRA scheme (Theys et al., 2021). Results show excellent agreement between both spectral bands, and demonstrate better sensitivity to SO₂ in BD2. In a second step, we conducted a wide sensitivity analysis for the SO₂ layer height retrieval, based on simulated spectra, covering various spectral, atmospheric, and observation conditions typical of TROPOMI measurements. We showed that, despite the increased noise and ozone absorption below 310 nm, the TROPOMI BD2 could be reliably used for SO₂ layer height retrievals, with retrieval errors reduced by at least a factor of 2 compared to BD3, the most important gain in sensitivity being in the UTLS. The sensitivity analyses also revealed that uncertainties on many input parameters, like the temperature, pressure and air density, the exact SO₂ profile shape, surface height, and albedo have moderate influence on the retrieval quality. In contrast, we found that uncertainty on the SO₂ absorption cross sections at short UV wavelengths can lead to non-negligible systematic biases on the SO₂ layer height of about 0.7 km. This would deserve more attention in the future. Unsurprisingly, we also found that uncertainties on ozone profiles significantly impact the retrieved layer heights. Overall, the estimated systematic errors, considering extreme and small perturbations of all the investigate parameters, were 1.91 km for the SO₂ LH and 2.58 % for the VCD, which is promising. To mitigate the impact of ozone, we derived mean ozone profiles for different ozone columns using the climatology of Lamsal et al. (2004) and generated an extensive look-up table of spectra covering more than 60 million conditions. Following numerical refinements, our algorithm was applied to real TROPOMI measurements for several volcanic events. The 2019 Raikoke eruption served as an ideal test case due to the high SO₂ abundances released in the UTLS. From this event, we identified an optimized fitting window (305–326 nm) that maximizes the sensitivity to SO₂ absorption while minimizing systematic errors due to ozone profile uncertainties. The detection limit was determined to be 2 DU, a significant improvement over the 5 DU of BD3 LUT-COBRA and 20 DU of the S-5P operational product. Overall, the retrieval results reveal a good consistency between both bands. We also analyzed the Etna eruption on 15 August 2024, a more challenging scenario with lower SO₂ amounts. The BD2 retrievals closely match those from BD3 LUT-COBRA and the operational TROPOMI product, but provide a more detailed spatial representation, due to the improved detection limit. Additional case studies, including the 2024 Ruang and 2019 Ulawun tropical eruptions, confirm that BD2 retrievals are of much better quality that in BD3. Retrievals in BD2 capture more dispersed SO₂ plumes with slightly lower LH estimates than in BD3 but in better agreement with independent IASI and MLS measurements. The SO₂ mass estimates from our algorithm are systematically higher than those from BD3 TROPOMI and IASI, likely due to the increased number of valid pixels (in comparison to BD3) and improved spatial coverage (compared to TIR measurements). The cases of Ruang and Ulawun also show that our SO₂ height retrievals outperform the operational product. The latter apparently fails to reproduce the SO₂ height at the tropical tropopause level. The main limitations likely stem from the intrinsically low vertical sensitivity of BD3 at these altitudes and the decreased sensitivity of the FP_ILM algorithm for lower SO₂ columns.

Future developments will focus on the following aspects:

1. We plan to extend our set of synthetic spectra by including volcanic aerosols (ash or sulphate) and assess their impact on the retrievals. As discussed in Section 3.2, aerosols may introduce non-negligible additional uncertainties in the retrieved SO₂ heights when collocated with the plume. The expansion of our code to analyze more complex volcanic eruptions (like the eruption of Hunga Tonga) and treating aerosols explicitly are a longer-term and challenging objective.

2. We also plan to combine the Optimal Estimation Method (OEM) of Rodgers (2000) with LUT-COBRA to further mitigate uncertainties related to ozone absorption and possibly other parameters. We have performed initial tests on synthetic spectra and found a considerable reduction of the LH and VCD biases.

3. By processing the full time-series of TROPOMI, we plan to study long-term degassing changes at the global scale, by exploiting the combination of sensitive SO₂ height and column retrievals. This would also contribute to improve our understanding of the global sulfur budget, and of volcanic and atmospheric processes. Beyond volcanic sources, we also intend to test our algorithm on anthropogenic SO₂ scenes.

4. It is an objective to extend the validation to more cases. One possibility is to compare the results with dispersion modelling results from the Plume_traj toolkit (Pardini et al., 2017, 2018; Queißer et al., 2019).

Appendix A: Sensitivity analyses of BD2 LUT-COBRA

The synthetic tests were conducted under the assumption that the background contribution was well-known (i.e., $\stackrel{\mathrm{‾}}{y}=y_{\mathrm{bckg}}$). From Eqs. (1) and (3), we thus have

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}y-y_{{\mathrm{SO}}_{\mathrm{2}},i}-\stackrel{\mathrm{‾}}{y}& =y_{{\mathrm{SO}}_{\mathrm{2}}}+y_{\mathrm{bckg}}+\mathit{ϵ}-y_{{\mathrm{SO}}_{\mathrm{2}},i}-\stackrel{\mathrm{‾}}{y}\\ & =y_{{\mathrm{SO}}_{\mathrm{2}}}-y_{{\mathrm{SO}}_{\mathrm{2}},i}+\mathit{ϵ}.\end{array}\end{array}$$

This means that SO₂ spectra, rather than total (SO₂ + O₃) spectra, were processed. In addition to the analyses described in Sect. 3.3, other tests were conducted using alternative look-up table conditions. Table A1 presents the changes made to the baseline conditions, along with the estimated retrieval errors for the same perturbed conditions as in Table 2, assuming a SO₂ height of 13.5 km and column density of 5 DU.

Table A1 Overview of the additional sensitivity tests of LUT-COBRA in BD2, considering different look-up table parameters. The results are given for a SO₂ column density of 5 DU and an altitude of 13.5 km.

Appendix B: Air mass factor (AMF) calculation and altitude-dependence

The SO₂ AMFs can be calculated as

$$\begin{array}{}\text{(B1)}& \mathrm{AMF}={\displaystyle \frac{-\ln \left(\frac{I_{\mathit{\lambda }}^{\mathrm{w}}}{I_{\mathit{\lambda }}^{\mathrm{wo}}}\right)}{{\mathit{\sigma }}_{\mathit{\lambda }}\mathrm{VCD}}},\end{array}$$

where $I_{\mathit{\lambda }}^{\mathrm{w}}$ and $I_{\mathit{\lambda }}^{\mathrm{wo}}$ correspond respectively to the radiance with and without the sulfur dioxide contribution, and σ_λ stands for the SO₂ absorption cross sections. Figure B1 shows SO₂ AMF vertical profiles for various wavelengths and SO₂ VCDs.

Figure B1 SO₂ AMF distributions as a function of the altitude, for different wavelengths and SO₂ VCDs.

Appendix C: SO₂ absorption cross section

Figure C1 makes a comparison between the SO₂ absorption cross sections from Bogumil et al. (2003), used in our look-up tables, and those from Birk and Wagner (2018).

Figure C1 Comparison of the SO₂ absorption cross sections from Bogumil et al. (2003) and Birk and Wagner (2018), for a temperature (T) of 293 K.