Aerosol Composition and Extinction of the 2022 Hunga Plume Using CALIOP

We use the CALIOP (Cloud-Aerosol Lidar with Orthogonal Polarization) instrument to determine the microphysical properties of the stratospheric aerosol plume after the Hunga eruption in 2022, the largest so far after the Pinatubo in 1991. In the early stages, low depolarization (<2 %) is found everywhere except in patches of high depolarization (up to 35 %) detected within the plumes of sulfur compounds up to 3 d after the eruption. As standard CALIOP L2 products are not operational in the case of the Hunga aerosol plume, we implement an iterative method of successive approximations to retrieve extinction profiles, by estimating the aerosol optical depth (AOD) and then the lidar ratio (LR). The AOD of the plume at 532 nm is between 0.5 and 1.25 on the first four days, then decreases rapidly and stabilizes at 0.047 ± 0.011 for March 2022. The LR is initially between 60 and 80 sr, consistent with the early growth of sulfate aerosol particles, and then decreases to 48 ± 6 sr between late January and late March 2022. Results are compared and validated with the solar occultation instrument SAGE III (Stratospheric Aerosol and Gas Experiment) on board the International Space Station (ISS) and Mie calculations. A comparison with limb-viewing instruments highlights significant quantitative disagreements in extinction and AOD estimates, which we attribute, in part, to the unusual size distribution of the aerosols within the Hunga plume.

1 Introduction

After becoming active on 20 December 2021, the Hunga submarine volcano (20.57° S, 175.38° W) produced a spectacular explosive phase, with a Volcanic Explosivity Index (VEI) of ∼6, close to Mount Pinatubo eruption in 1991 (Poli and Shapiro, 2022), on 15 January 2022 by 04:10 UTC (Poli and Shapiro, 2022; Gupta et al., 2022; Yuen et al., 2022). This eruption injected material up to an altitude of ∼58 km by 04:30 UTC, with most of the plume settled between 26 and 34 km altitude afterwards (Carr et al., 2022; Khaykin et al., 2022; Podglajen et al., 2022). In the first hours, volcanic ash, ice and sulfur dioxide were seen in the stratosphere by geostationary satellite instruments (Himawari-8, GOES-17) before giving way to a plume where sulfate aerosols (SA) and water vapor coexisted for a few weeks (Sellitto et al., 2022; Carn et al., 2022; Legras et al., 2022). This was due to the rapid conversion of sulfur dioxide, a process enhanced by the abundant presence of water vapor (Millán et al., 2022; Zhu et al., 2022; Sellitto et al., 2024). The injected water vapor was estimated at 146 ± 5 Tg, which means an instantaneous 10 % increase of the global stratospheric content (Millán et al., 2022; Khaykin et al., 2022). As a result, the small, spherical, non-depolarizing liquid particles detected by spaceborne and in-situ observations (Legras et al., 2022; Kloss et al., 2022; Baron et al., 2023) in the days following the eruption formed a long-lived SA plume that circled the globe multiple times, spreading slowly in latitudes between the southern polar region and 20° N (Taha et al., 2022; Legras et al., 2022; Khaykin et al., 2022; Schoeberl et al., 2023b; Sellitto et al., 2024) and still detected at the end of 2023 (Knepp et al., 2024). The particle size distribution (PSD) of the Hunga plume exhibited an effective radius close to 0.4 µm (Khaykin et al., 2022; Baron et al., 2023; Duchamp et al., 2023; Asher et al., 2023; Boichu et al., 2023; Knepp et al., 2024), a mode width σ∼1.25 (Duchamp et al., 2023; Knepp et al., 2024) and a total density number of 3–5 cm⁻³ (Duchamp et al., 2023). This mode established quickly, within two weeks, and persisted many months. It differs from the background but also from other recent stratospheric eruptions of medium intensity (Wrana et al., 2023) and established much faster than after the Pinatubo eruption (Boichu et al., 2023). Considering these PSD parameters and a H₂SO₄ weight proportion of 70 %, Duchamp et al. (2023) estimated the total mass of H₂SO₄ in SA at 0.66 ± 0.1 Tg for the Hunga. This corresponds to an initial SO₂ emission of 0.44 Tg in the stratosphere, in line with modest early 0.4–0.5 Tg estimates (Millán et al., 2022; Carn et al., 2022; Asher et al., 2023) compared to the 18 ± 4 Tg of Pinatubo (Guo et al., 2004) and even to the 1.5 Tg of Raikoke (de Leeuw et al., 2021), with a VEI of 4 (Firstov et al., 2020). However, recent studies suggested that the injected SO₂ mass burden of the Hunga eruption could have been larger, with 1.0–1.3 Tg SO₂ column estimates using satellite data (Sellitto et al., 2024), geological processes (Wu et al., 2025) and modeling (Bruckert et al., 2025). The monomodal distribution assumption used for the retrieval of aerosol PSD with satellite data could lead, in the event of a second peak in the actual PSD, to a significant underestimation of the total number density and therefore of the total mass (Wrana et al., 2025). Sellitto et al. (2022) showed that water vapor warming impact on the top of atmosphere radiative forcing dominated over aerosol cooling effect, leading to a – very unusual for a volcanic eruption – net warming of the climate system two weeks after the eruption. In the longer term, other studies estimated a dominant aerosol effect slightly decreasing surface temperatures in the Southern Hemisphere in 2022–2023 (Schoeberl et al., 2023a; Gupta et al., 2025). The Hunga aerosol layer is likely 1.5–10 times more effective at producing a net radiative cooling than those resulting from the El Chichón and Pinatubo eruptions (Sellitto et al., 2025). The aim of this study is to take advantage of the properties of the CALIOP lidar to better characterize the composition of the young plume, and to implement a method for reconstructing the AOD, LR and extinction profiles over the first three months following the eruption, before comparing and validating them with others products.

2 Data and Methods

2.1 CALIOP and Optical Properties Retrieval

The Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) is a two-wavelength lidar on board the CALIPSO satellite (Vaughan et al., 2004; Winker et al., 2010). In this work, we use the L1 v4.51 532 nm total attenuated backscatter averaged in 40 km bins in the horizontal, filtering all profiles which have been contaminated by a low-energy shot. The native vertical resolution is preserved. The molecular backscatter – used to compute the 532 nm attenuated scattering ratio (ASR) – is calculated following Hostetler et al. (2006). The exceptional case of Hunga, due to the height and intensity of its plume (well separated from the tropopause and with good signal-to-noise ratio for night orbits), allows us to retrieve the optical properties of the plume over the first few weeks after the eruption. First, the AOD is estimated from the ASR at aerosol-free altitudes above and below the plume. Then, we calculate the corresponding first guess LR by dividing the AOD by the vertically integrated attenuated backscatter after subtracting molecular backscatter. We use an iterative method of successive approximations to remove the aerosol optical attenuation from above the aerosol layer and retrieve the aerosol extinction profiles with the final LR. The details of the algorithm are given in the next section (Sect. 2.2). LR depends on some microphysical properties of the observed aerosols such as size, shape, and chemical composition (Fernald et al., 1972). It usually varies from about 5 to 100 sr, with low values generally associated with large particles and higher values suggesting the presence of small and/or absorbing particles, making it useful for aerosol type classification (Müller et al., 2007; Kim et al., 2018). We present daily results based both on individual orbits (Fig. 4) and on daily averages (Fig. 5c). For each day, the chosen individual orbit, termed the “most significant orbit of the day”, is the one containing the plume with the highest mean ASR above 20 km and over a threshold value of 1.8. Then an average over a selected latitude range is performed inside the plume. On the other side, daily zonal averages are calculated by discarding orbits where a plume is detected but processing is not possible due to a lack of exploitable signal below the plume, that is a range free of aerosols (this case appears progressively from the end of February). Plume free orbits, frequent in the early stage, are counted as zero in the average of the AOD. We limit ourselves to night-time data, as daytime data are more noisy. Measurements in the South Atlantic Anomaly (SAA) are excluded throughout the study, between 90° W and 60° E in the Southern Hemisphere, where the Earth's magnetic field disturbs CALIOP measurements (Noel et al., 2014). Due to solar activity, CALIOP was not operating on 18 January and between 20 and 26 January. This procedure palliates standard CALIOP aerosol L2 products which unfortunately are not suitable for the type of event considered in this study (see Appendix A). Note that throughout the manuscript, references to CALIOP-derived AOD correspond to the plume-layer AOD, i.e., the aerosol optical depth integrated only over the vertical extent of the detected plume. This differs from the usual column-integrated AOD, but this distinction has little impact on the results because the background extinction outside the plume in the stratosphere is very small compared with the plume extinction.

2.2 Detailed Optical Properties Retrieval

The lidar equation, after correction of the range effect, is:

$$\begin{array}{}\text{(1)}& {\mathit{\beta }}^{\prime }(\mathit{\lambda },z)=\mathit{\beta }(\mathit{\lambda },z)\exp \left(-\mathrm{2}\underset{z}{\overset{+\mathrm{\infty }}{\int }}\mathit{\sigma }(\mathit{\lambda },z^{\prime })\mathrm{d}{z}^{\prime }\right)\end{array}$$

where β is the backscatter coefficient, σ is the extinction coefficient and β^′ is the retrieved attenuated backscatter stored in the standard CALIOP L1B data product.

The two coefficients are separated into a molecular, a particle and an ozone component: $\mathit{\beta }={\mathit{\beta }}_{\mathrm{m}}+{\mathit{\beta }}_{\mathrm{p}}$ and $\mathit{\sigma }={\mathit{\sigma }}_{\mathrm{m}}+\mathit{\eta }{\mathit{\sigma }}_{\mathrm{p}}+{\mathit{\sigma }}_{{\mathrm{O}}_{\mathrm{3}}}$ where η is the multi-scatter parameter (the O₃ backscatter is neglected but not the extinction due to absorption). Molecular and particle LR are respectively $S_{\mathrm{m}}={\mathit{\sigma }}_{\mathrm{m}}/{\mathit{\beta }}_{\mathrm{m}}$ and $S_{\mathrm{p}}={\mathit{\sigma }}_{\mathrm{p}}/{\mathit{\beta }}_{\mathrm{p}}$. According to Hostetler et al. (2006, p. 28–29), we can solve the molecular part using ${\mathit{\beta }}_{\mathrm{m}}(\mathit{\lambda },z)=\frac{\mathrm{3}}{\mathrm{8}\mathit{\pi }{k}_{\text{bw}}\left(\mathit{\lambda }\right)}{\mathit{\sigma }}_{\mathrm{m}}(\mathit{\lambda },z)$ and ${\mathit{\sigma }}_{\mathrm{m}}(\mathit{\lambda },z)={\mathcal{N}}_{\mathrm{m}}\left(z\right)Q_{\mathrm{s}}\left(\mathit{\lambda }\right)$ with k_bw(532 nm)=1.0313 and $Q_{\mathrm{s}}\left(\mathrm{532}\mathrm{nm}\right)=\mathrm{5.167}\times {\mathrm{10}}^{-\mathrm{31}}{\mathrm{m}}^{\mathrm{2}}$, where 𝒩_m is the molecular number density. ${\mathit{\sigma }}_{{\mathrm{O}}_{\mathrm{3}}}(\mathit{\lambda },z)={\mathcal{N}}_{{\mathrm{O}}_{\mathrm{3}}}\left(z\right)Q_{{\mathrm{O}}_{\mathrm{3}}}(\mathit{\lambda },T)$, where ${\mathcal{N}}_{{\mathrm{O}}_{\mathrm{3}}}$ is the ozone number density and $Q_{{\mathrm{O}}_{\mathrm{3}}}$ is the ozone cross-section, dependent on temperature T, measured by Serdyuchenko et al. (2014). The attenuated scattering ratio is defined as $\text{SR}={\mathit{\beta }}^{\prime }/{\mathit{\beta }}_{\mathrm{m}}$ at 532 nm. For a pure molecular profile, it is therefore:

$${\text{SR}}_{\mathrm{m}}\left(z\right)=\exp \left(-\mathrm{2}\underset{z}{\overset{+\mathrm{\infty }}{\int }}\left({\mathit{\sigma }}_{\mathrm{m}}\right(\mathit{\lambda },z^{\prime })+{\mathit{\sigma }}_{{\mathrm{O}}_{\mathrm{3}}}(\mathit{\lambda },z^{\prime }\left)\right)\mathrm{d}{z}^{\prime }\right)$$

which is always lower than 1.

We now define Γ as the ratio of SR to SR_m:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle }\mathrm{\Gamma }\left(z\right)=& {\displaystyle \frac{\text{SR}\left(z\right)}{{\text{SR}}_{\mathrm{m}}\left(z\right)}}\\ {\displaystyle }=& {\displaystyle }\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{p}}\left(z\right)}{{\mathit{\beta }}_{\mathrm{m}}\left(z\right)}}\right)\exp \left(-\mathrm{2}{\int }_z^{+\mathrm{\infty }}\mathit{\eta }{\mathit{\sigma }}_{\mathrm{p}}\left(z^{\prime }\right)\mathrm{d}{z}^{\prime }\right)\end{array}\end{array}$$

In practice, +∞ is the top of the vertical domain of CALIOP detection, that is 40 km, above which attenuation is assumed to be negligible.

We consider now an aerosol layer which is isolated between a bottom altitude z_B and a top altitude z_T, that is there is no aerosol above z_T and there is a layer without aerosols below z_B. The AOD of this layer is ${\int }_{z_{\mathrm{B}}}^{z_{\mathrm{T}}}\mathit{\eta }{\mathit{\sigma }}_{\mathrm{p}}\left(z^{\prime }\right)\mathrm{d}{z}^{\prime }$. As β_p(z)=0 for z≤z_B and z≥z_T, the AOD can be obtained as:

$$\begin{array}{}\text{(3)}& \text{AOD}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\ln \left({\displaystyle \frac{\mathrm{\Gamma }\left(z_{\mathrm{B}}^{-}\right)}{\mathrm{\Gamma }\left(z_{\mathrm{T}}^{+}\right)}}\right)\end{array}$$

where $\mathrm{\Gamma }\left({z_{\mathrm{B}}}^{-}\right)$ and $\mathrm{\Gamma }\left({z_{\mathrm{T}}}^{+}\right)$ stand, respectively, for a layer free of aerosol below and above the aerosol layer. In practice, for an aerosol layer with no other aerosols above, $\mathrm{\Gamma }\left(z_{\mathrm{T}}^{+}\right)$ can be set to 1. For z≤z_B, we require a layer of Γ fluctuating around a minimum which is thicker than 1 km to decide that it can be considered as free of aerosols and usable to measure the AOD. Of course, this is liable to artifacts especially for small AOD and therefore our estimate in such cases must be taken as a lower bound.

It is also necessary to correct Γ for the fact that the background is not totally free of aerosols and for calibration issues. Unlike Martinsson et al. (2022) who use a constant correction factor but for lower altitudes, we estimate it from the attenuated scattering ratio observed in 2021 in a state of CALIOP laser similar to 2022 and for a weakly perturbed stratosphere (besides the remaining aerosols of the 2019–2020 Australian wildfires). For each Hunga case, the 2021 attenuated scattering ratio is averaged at the same latitude and time of the year between 17 km and the altitude of the bottom of the plume. This average is used as a reference for the AOD calculation from the same average made in 2022 under the plume introducing a corrective factor. This factor usually larger than one divides $\mathrm{\Gamma }\left(z_{\mathrm{B}}^{-}\right)$. As the background aerosols are possibly replaced and then missing within the Hunga plume, their effect might not need to be compared and therefore we show results obtained with and without the correction (Fig. 5).

Once the AOD is known, and assuming that the particle lidar ratio S_p is constant over the aerosol layer (which means uniform microphysical properties), we have:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\displaystyle }C\left(z\right)& {\displaystyle }={\displaystyle \frac{\mathrm{\Gamma }\left(z\right)}{\mathrm{\Gamma }\left(z_{\mathrm{T}}\right)}}\\ {\displaystyle }& {\displaystyle }=\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{p}}\left(z\right)}{S_{\mathrm{p}}{\mathit{\beta }}_{\mathrm{m}}\left(z\right)}}\right)\exp \left(-\mathrm{2}{\int }_z^{z_{\mathrm{T}}}\mathit{\eta }{\mathit{\sigma }}_{\mathrm{p}}\left(z^{\prime }\right)\mathrm{d}{z}^{\prime }\right)\end{array}\end{array}$$

where C(z) is known from the observations and the molecular properties.

For a given S_p, this equation can be solved from the top. Since S_p is not known, a determination of S_p and σ_p can be done iteratively using:

$$\begin{array}{}\text{(5)}& S_{\mathrm{p}}{\mathit{\beta }}_{\mathrm{m}}\left(z\right)\left(C\right(z)\exp \left(\mathrm{2}\underset{z}{\overset{z_{\mathrm{T}}}{\int }}\mathit{\eta }{\mathit{\sigma }}_{\mathrm{p}}\left(z^{\prime }\right)\mathrm{d}{z}^{\prime }\right)-\mathrm{1})={\mathit{\sigma }}_{\mathrm{p}}\left(z\right)\end{array}$$

The first iteration neglects aerosol attenuation for a first guess of S_p:

$$\begin{array}{}\text{(6)}& \mathit{\eta }{S}_{\mathrm{p}}\underset{z_{\mathrm{B}}}{\overset{z_{\mathrm{T}}}{\int }}{\mathit{\beta }}_{\mathrm{m}}\left(z\right)\left(C\right(z)-\mathrm{1})\mathrm{d}z=-\text{AOD}\end{array}$$

This integral can be calculated using the trapezoidal rule over a set of discretized levels. The following iterations calculate σ_p(z) downward from z_T using Eq. (5). If the vertical discretization is defined, counting from the top, with a step Δz and if z_T and z_B are, respectively, the levels i and j with j>i, the first part of the iteration calculates the extinction and the attenuation factor d by initializing them as:

$$\begin{array}{rl}{\displaystyle }d\left(i\right)& {\displaystyle }=\mathrm{1}\\ {\displaystyle }{\mathit{\sigma }}_{\mathrm{p}}\left(i\right)& {\displaystyle }=S_{\mathrm{p}}{\mathit{\beta }}_{\mathrm{m}}\left(i\right)\left(C\right(i)-\mathrm{1})\end{array}$$

and then for k in $[i+\mathrm{1},j]$:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\displaystyle }d\left(k\right)& {\displaystyle }=d(k-\mathrm{1})\exp \left(\mathrm{2}\mathit{\eta }{\mathit{\sigma }}_{\mathrm{p}}(k-\mathrm{1})\mathrm{\Delta }z\right)\\ {\displaystyle }{\mathit{\sigma }}_{\mathrm{p}}\left(k\right)& {\displaystyle }=S_{\mathrm{p}}{\mathit{\beta }}_{\mathrm{m}}\left(k\right)\left(C\right(k\left)d\right(k)-\mathrm{1})\end{array}\end{array}$$

The second part of the iteration recalculates S_p as:

$$\begin{array}{}\text{(8)}& S_{\mathrm{p}}^{*}={\displaystyle \frac{\text{AOD}}{\mathit{\eta }}}{\left[\sum _{k=i}^j{\mathit{\beta }}_{\mathrm{m}}\left(k\right)(\mathrm{1}-C(k\left)d\right(k\left)\right)\mathrm{\Delta }z\right]}^{-\mathrm{1}}\end{array}$$

where the sum is performed using the trapezoidal rule.

The new estimate of S_p is then obtained by combining the previous estimate with Eq. (8) as:

$$\begin{array}{}\text{(9)}& S_{\mathrm{p}}\leftarrow {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(S_{\mathrm{p}}+S_{\mathrm{p}}^{*}\right)\end{array}$$

The iteration is made until S_p converges. The convergence is greatly accelerated by Eq. (9) and is in practice obtained in less than 8 steps for all investigated cases.

In this study, we neglect multiple scattering by setting η=1.

Figure 1 shows an example (Cloud C1 from 17 January) of the procedure described above from the attenuated scattering ratio (Fig. 1a–c; computed with the attenuated backscatter stored in the L1B files) to the aerosol extinction σ_p (Fig. 1g). The standard deviation from the mean (Fig. 1j–l) is chosen to calculate the error bars for the individual cases. For this specific case, we calculate an AOD of 1.24 ± 0.13 and then a LR of 70.9 ± 2.5 sr.

Figure 1 (a) CALIOP 532 nm attenuated scattering ratio SR from orbit 2022-01-17T17-24-22-ZN on 17 January 2022. (b) Same as (a) but averaged between 18 and 14° S. (c) Horizontal section of (a) at 17.5 km altitude. (d) Blue curve: same as (b) (SR, no molecular attenuation). Orange curve: (b) taking into account molecular attenuation including O₃ transmission (Γ, Eq. 2). (e) Same as (d) with a reduced scale to better see variations in the vicinity of 1. (f) Orange curve: same as (e). Blue dotted line: corresponds to $\mathrm{\Gamma }\left(z_{\mathrm{B}}^{-}\right)$. (g) Aerosol extinction σ_p retrieved together with S_p. (h) One-way attenuated aerosol extinction ${\mathit{\sigma }}_{\mathrm{p}}^{\prime }$ which is ${\mathit{\sigma }}_{\mathrm{p}}/\sqrt{d}$. (i) Blue curve: measured CALIOP attenuated scattering ratio. Orange curve: reconstructed attenuated scattering ratio from retrieved σ_p, S_p and d. Green curve: reconstructed non attenuated scattering ratio from retrieved σ_p and S_p. (j) Mean ($\mathrm{\Gamma }\left(z_{\mathrm{B}}^{-}\right)$, blue dotted line), median (red dotted line) and standard deviation (blue area) of Γ calculated below the aerosol layer. (k) and (l) are respectively the corresponding aerosol extinction σ_p and LR and AOD from Γ values (j).

2.3 SAGE III/ISS

The Stratospheric Aerosol and Gas Experiment III (SAGE III) instrument, on board the International Space Station (ISS), has been providing measurements of solar occultations since June 2017 (Cisewski et al., 2014). The instrument provides aerosol extinctions at nine wavelengths from 384 to 1543 nm in 0.5 km vertical steps between 0 and 45 km altitude. The instrument observes about 15 sunrises and 15 sunsets per day with a latitudinal range which varies depending on the period of year. The aerosol extinctions are retrieved as residuals of a spectral multilinear fit for O₃ and NO₂ but do not require any size distribution assumptions unlike instruments with limb-scatter geometry like OSIRIS (Bourassa et al., 2007) and OMPS-LP (Loughman et al., 2018). Due to its low geographical sampling, SAGE III only began to observe the Hunga SA plume from February 2022, with several profiles on only 3 d (7, 15 and 16 February), whereas the plume was detected for most of the days in March 2022. For this study, we use the version 5.3 of the SAGE III/ISS level 2 solar aerosol product.

2.4 OMPS-LP

The Ozone Mapping and Profiler Suite Limb Profiler (OMPS-LP) instrument, on board the Suomi-NPP satellite, provides aerosol extinction profiles from measurements of the scattered solar radiation in 1 km vertical steps. In this study, we use the 745 nm band from versions v2.1 and v2.5 of the aerosol extinction retrieval algorithm developed by NASA (Taha et al., 2021) and two other OMPS-LP products: the 745 nm band of the v1.3 aerosol extinction product developed by the University of Saskatchewan (USASK) (Bourassa et al., 2023) and the 869 nm band of the v2.1 aerosol extinction coefficient product developed by the Institute of Environmental Physics (IUP) after a conversion to 745 nm (Rozanov et al., 2024). Of the three slits separated by 250 km at the tangent point designed to observe the atmosphere (Jaross et al., 2014), we use the middle one as it has better straylight performance and pointing knowledge. The extinction is averaged daily over all orbits of that day outside the SAA zone and after a horizontal interpolation to a standard latitude grid of 1.1° resolution that corresponds to the mean resolution of OMPS-LP in the considered range of latitudes.

2.5 OSIRIS

The Optical Spectrograph and InfraRed Imaging System (OSIRIS), on board the Odin satellite, measures vertical profiles of spectrally dispersed, limb scattered sunlight from the upper troposphere to the lower mesosphere (Llewellyn et al., 2004). In this work, we use v7.3 of the level 2 Odin/OSIRIS stratospheric aerosol extinction coefficient profile at 750 nm (Rieger et al., 2019) using a multiwavelength retrieval that improves the accuracy of the extinction product by reducing sensitivity to the unknown PSD in the inversion. This product provides aerosol extinction coefficient in 1 km vertical steps between 0.5 and 44.5 km altitude.

2.6 Mie Calculations

We convert 532 nm CALIOP AOD results to 756 nm and obtain the theoretical LR as a function of the PSD (see Sect. 3.2) by calculating extinction and backscatter coefficient using miepython (see Data availability), a python code implementing the Mie theory according to Wiscombe (1979). As input, we use fixed real part refractive indices from the GEISA database (Armante et al., 2016), for a temperature of 215 K and a H₂SO₄ weight proportion w_s of 70 % (Tabazadeh et al., 1997; Duchamp et al., 2023) considering that the rest of the liquid droplet is water (Biermann et al., 2000), that is $\mathit{\{}{n}_{\mathrm{532}},n_{\mathrm{756}}\mathit{\}}=\mathit{\{}\mathrm{1.439},\mathrm{1.438}\mathit{\}}$. As the SA have a very low absorption in the shortwave spectral range (Palmer and Williams, 1975), we fix the imaginary part of the refractive index to 10⁻⁶. Unimodal lognormal size distribution is assumed for the Hunga SA plume following the results of Duchamp et al. (2023). From this article (Fig. 2), we derive median radius r_m=0.35 µm and mode width σ=1.25, producing a conversion factor from 532 nm results to 756 nm of 0.815.

3 Results

3.1 Plume Composition

The Hunga eruption injected volcanic ash, SO₂ and water vapor into the stratosphere (Millán et al., 2022; Legras et al., 2022). The ash/ice component rapidly dissipated, likely due to sedimentation (Sellitto et al., 2022; Khaykin et al., 2022). A remaining ash or ice plume was observed by CALIOP on 15 January between 33 and 39 km altitude with a high depolarization ratio (over 30 %, see Fig. 2a), then it was observed again on 16 January by CALIOP and by a ground-based lidar over Reunion island on 20 January (Baron et al., 2023). After that date, it could be followed only by OMPS-LP (Taha et al., 2022). The main plume of sulfur compounds was revealed by geostationary imagers (Himawari-8 and MSG-1, see Appendix B) and IASI observations within a day of the main eruption suggesting two separate aerosol clouds named C1 and C2 (Sellitto et al., 2022; Legras et al., 2022), moving westward, which were also observed by CALIOP (Fig. 2b) and by UV satellite instruments (Carn et al., 2022). These observations suggested that the conversion of SO₂ into sulfates was very fast for the two clouds, especially for the moistier C1 (Sellitto et al., 2022). Both clouds exhibited a high scattering ratio with very low depolarization in their bulk indicating the predominance of small spherical particles. Nevertheless, two small patches were observed, in the lower part of C2, of higher depolarization close to 35 % with a large color ratio (CR), which suggested remaining ice crystals or ash.

Figure 2 6 CALIOP night orbits on 15 January (a), 16 (b), 17 (c, d) and 19 (e, f). For each panel, from left to right, are represented the 532 nm ASR, the 532 nm depolarization ratio (DR in %, orthogonal channel/parallel) and the 1064 nm/532 nm color ratio (CR, 1064 nm backscatter coefficient is not available above 30 km). The corresponding CALIOP ground tracks on RGB-Ash maps from Himawari-8 and MSG-1 are provided in Appendix B.

Two successive orbits from 17 January (Fig. 2c and d) revealed, in turn, the two clouds exhibiting maximum depolarization values of 4 %–5 % within their main structures, which could have corresponded to dissipating ice or falling ash. Some other pockets of depolarization were observed in other aerosol clouds below 18 km altitude (not shown). The higher depolarization spots also exhibited a larger CR compared with other plume sections, consistent with the dissipation of solid particles (Prata et al., 2017). At this stage, the aerosol index from UV satellite instruments did not reveal any ash in the two clouds (Taha et al., 2022).

The two clouds were seen again on 19 January (Fig. 2e and f) with depolarization everywhere below 2 %, before the interruption of CALIOP operations. During this interruption, the two clouds elongated under the zonal shear (Legras et al., 2022) and were no longer easily distinguishable after 27 January when CALIOP resumed operation. All the subsequent orbits showed very low depolarization and CR, consistent with a plume containing only small spherical particles.

We note an increase in the CR towards the bottom of plume, by a factor 2 to 4 (Fig. 2c and d). Although this feature could reflect an actual change in the particle size, it is quantitatively consistent with differential attenuation between the two wavelengths keeping a constant size distribution and we interpret it as such.

Figure 3 shows that the apparent CR calculated from attenuated backscatter increases with the partial optical depth from the top of the plume and so may vary significantly across the plume when the AOD (that is the total optical depth) is large, even with uniform optical properties. For instance, this can be seen for cloud C1 on 17 January, with theoretical CR values of about 0.4 for small optical depths (top of the cloud), increasing with optical depth (toward the bottom of the cloud) and reaching values of 1.2 and 1.75 (green and orange PSDs, respectively) at an optical depth of 1.24 (CALIOP AOD value for this case) in agreement with the observed CR values provided by CALIOP (Fig. 2d).

Figure 3 Theoretical CR (1064$/$532 nm) calculated with a Mie code (see Sect. 2.6) as a function of the partial optical depth from the top of the cloud for possible PSD assumptions in the early plume (r_eff=0.2 µm with σ=1.6 in blue, r_eff=0.3 µm with σ=1.4 in orange and r_eff=0.4 µm with σ=1.3 in green).

3.2 Optical Properties

Figure 4 shows the temporal evolution of the 532 nm AOD (a) and LR (b) along selected orbits over the first few weeks of plume dispersion in the stratosphere. Early AOD values were:

1. 16 January: ${\text{AOD}}_{\text{C1}}=\mathrm{0.51}\pm \mathrm{0.07}$; ${\text{AOD}}_{\text{C2}}=\mathrm{0.59}\pm \mathrm{0.08}$

2. 17 January: ${\text{AOD}}_{\text{C1}}=\mathrm{1.24}\pm \mathrm{0.13}$; ${\text{AOD}}_{\text{C2}}=\mathrm{1.01}\pm \mathrm{0.12}$

3. 19 January: ${\text{AOD}}_{\text{C1}}=\mathrm{0.62}\pm \mathrm{0.08}$; ${\text{AOD}}_{\text{C2}}=\mathrm{0.89}\pm \mathrm{0.09}$

Using a similar method for 16 January, Sellitto et al. (2022) obtained values of ≈0.55 for C1 and ≈0.60 for C2. On 19 January, a peak in AOD, reaching ≈1 at 500 nm and lasting about 12 h, was recorded by ground-based photometric observations from the AERONET station of Learmonth, western Australia (Boichu et al., 2023). CALIOP measurements then were interrupted but a 532 nm ground-based lidar located on Reunion Island (21.08° S, 55.38° E) provided Hunga plume AOD maximum values of 0.84 and 0.55 respectively on 21 and 22 January (Baron et al., 2023) also detected in AERONET stations measurements on Reunion Island (Boichu et al., 2023). The values of the ground measurements were consistent with the CALIOP values recorded before (0.89 and 0.62 for 19 January) and after the interruption (0.15 and 0.24, respectively, for 27 and 28 January). AOD values rapidly decreased during the first few weeks after the eruption. In March 2022, the AOD curve tended to stabilize at a monthly value of 0.047 ± 0.011 (error calculations with independent data, see Appendix C).

The LR calculations fluctuated mainly between 40 and 60 sr between 27 January and 25 March, with a slightly downward trend over time, resulting in a mean value of 48 ± 6 sr over this period. LR values, like AOD values, were highest during the four days following the eruption, when the plume was compact and localized, with values between 60 and 80 sr for C1 and C2. The Reunion Island lidar measured early LR within the Hunga peak aerosol plume between 68 and 75 sr at 532 nm between 21 and 23 January (Kloss et al., 2022; Baron et al., 2023).

Some properties of the Hunga plume were comparable to those produced by other volcanic eruptions. The two-month average of the LR was similar to the 48 sr value at 532 nm and 16 km altitude observed for the Nabro SA plume in June 2011 (Sawamura et al., 2012) and to the unconstrained assumed value 50 sr used in CALIOP L2 retrieval. One hypothesis for the high early LR values is the presence of residual ash (see Sect. 3.1), as ash particles tend to increase the LR (Vaughan et al., 2021). For instance, Prata et al. (2017) reported mean LR values of 69 ± 13 sr from CALIOP observations at a mean plume-top altitude of 12.5 km for the Puyehue–Cordón Caulle eruption in June 2011, while Lopes et al. (2019) retrieved LR values of 76 ± 27 sr between 18 and 19.3 km altitude for the Calbuco eruption in April 2015. Both eruptions were characterized by substantial ash loads. However, in the present case, the low depolarization values indicate that ash was unlikely to be present in significant amounts in the bulk of the plume, even during the early period. We therefore tentatively attribute the high early LR values primarily to microphysical differences in the SA, namely a smaller effective radius while maintaining a narrower size distribution, consistent with the initial particle growth reported in recent studies based on both measurements and modeling (Asher et al., 2023; Boichu et al., 2023; Kahn et al., 2024; Bruckert et al., 2025), and supported by Mie calculations (Fig. 6).

Figure 4 (a) Evolution of the AOD retrieved from CALIOP 532 nm attenuated backscatter between 16 January and 25 March 2022. Blue and orange circles correspond respectively to C1 and C2 early clouds. When we can no longer easily distinguish these two plumes, from 27 January onwards, we use the green circles to represent the case of the most significant orbit of the day over all the night orbits of the day outside the SAA zone (see Sect. 2.1). (b) Same as (a) for LR. The error bars come from standard deviation to the mean calculations (see Fig. 1j).

In order to test our method, we compare the results with extinction profile data from passive spaceborne sensors at coinciding measurement points. Figure 5a and b show collocated individual profiles of CALIOP, SAGE III and OMPS-LP on 15 February near the equator and 20 March near 20° S. Such coincidences are rare events over the period due to the scarce spatiotemporal sampling of SAGE III profiles (see Sect. 2.3). In both cases, the very good agreement between CALIOP and SAGE III extinction profiles is striking in terms of profile shape and values. Data from the three OMPS-LP products agree regarding the altitude of the extinction peak, but significantly underestimate its magnitude compared to direct extinction measurements from SAGE III, which we take as a reference (see Sect. 2.3). The discrepancy arises from the limited sensitivity of OMPS-LP in the early stages, as its long limb path prevents it from seeing into the plume effectively and the retrieval dependency on particle size assumptions (Bourassa et al., 2023). The gap becomes smaller with time as the plume gets diluted, as seen by comparing 20 March with 15 February. The IUP v2.1 and NASA v2.5 exhibit higher peak values than the other versions, bringing them closer in magnitude to SAGE III and CALIOP. Below the peak, the larger extinction in NASA v2.1 is due to an artifact originating from a convergence problem in the algorithm inducing a larger AOD when integrated. This issue was addressed in NASA v2.5 by substantially increasing the number of iterations (Taha et al., 2025), and was reduced in v2.5, although not completely eliminated. These results are consistently observed in other cases of coincidences (see Appendix D).

In Fig. 5c, we compare the evolution of AOD daily zonal means with the same products of Fig. 5a and b, adding the observations from the OSIRIS instrument (see Sect. 2.5). The figure also shows the SAGE III average over March 2022. The SAA longitudes are excluded for all instruments but SAGE III. Until early March, all curves show the same variations with two wide gaps corresponding to two periods where the bulk of the plume was located in the SAA region. CALIOP values on 16 January are lower than those on 17 and 19 January due to sampling (also visible in Fig. 4). Over January and early February, the 756 nm CALIOP converted curve is above the others with the OMPS-LP NASA v2.5 providing the closest match – except on 3 and 4 February, when enhanced values were detected in all limb datasets. In the second half of February, the converted CALIOP curve and the OMPS-LP curves agree within the uncertainties which range between 15 % and 50 % for CALIOP (see Fig. C1) and, with both OMPS-LP NASA versions, above all other curves. After 1 March, when the dispersion removes the SAA dependency, CALIOP AOD falls below the OMPS-LP NASA data which remain nearly at a constant level – with v2.5 exhibiting, however, larger variability than v2.1 – and match SAGE III. This may be due to the difficulty to measure AOD on many CALIOP orbits as the plume gets diluted, weakening the average signal, while the higher sensitivity of OMPS-LP still produces reliable results. However, the CALIOP curve falls within the other OMPS-LP products and, when the AOD is large (panel b), individual profiles still match SAGE III. Based on these results, we fit an exponential-decay curve to CALIOP before March and SAGE III data in March, suggesting an overall decay of the AOD with an e-folding time of 19.3 d. The use of a correction factor (see Sect. 2.2) for CALIOP data does not significantly change the results (orange vs. ochre curves). The OMPS-LP values of AOD, and differences between the three products, are consistent with existing literature (Sellitto et al., 2022; Taha et al., 2022; Khaykin et al., 2022; Bourassa et al., 2023) and their differences with SAGE III are on average smaller than with individual profiles for which the low sensitivity bounds the values. The AOD curve from OSIRIS is far below the others with values ∼0.01. The substantial low bias in OSIRIS aerosol extinction during the Hunga time period has been previously noted by Rozanov et al. (2024). It is believed to be caused by PSD assumptions in the retrieval that are amplified relative to OMPS-LP by the higher solar zenith angles under which OSIRIS operates. The unusual aerosol size distribution of the Hunga plume leads to the significant quantitative disagreements in AOD estimates from scattered light measurements due to the different size distribution assumptions, as well as the sensitivity of the instrument.

Figure 5 (a) Collocated profiles from CALIOP, SAGE III and OMPS-LP on 15 February. CALIOP profile is averaged between 4 and 2° S at 178.5° E and 14:43 UTC. SAGE III profile is measured at 3.2° S, 178.4° E and 18:17 UTC. OMPS-LP profile is averaged between 4 and 2° S at 171.7° W and 00:54 UTC. (b) Collocated profiles from CALIOP, SAGE III and OMPS-LP on 20 March. CALIOP profile is averaged between 21.3 and 19.3° S at 87.2° E and 20:40 UTC. SAGE III profile is measured at 20.3° S, 88.1° E and 00:15 UTC (21 March). OMPS-LP profile is averaged between 21.3 and 19.3° S at 95.5° E and 07:16 UTC. CALIOP night orbits and OMPS-LP day orbits will always have a time lag of about half a day. To compensate for this, we choose a different spatial band between CALIOP and OMPS-LP, taking into account the ERA5 angular velocity using panels (d) and (f) from Fig. 1 in Legras et al. (2022). Shading corresponds to the two OMPS-LP profiles that are averaged over the 2° latitude band. (c) Comparison of daily averaged AOD retrieval excluding profiles in the SAA longitudes from CALIOP (see Sect. 2.1) with several products from passive sensors measuring extinction (SAGE III, OMPS-LP and OSIRIS). Calculations are made between 17.5 and 31.5 km altitude and averaged between 20 and 10° S. The ochre curve represents CALIOP AOD without correction (see Sect. 2.2). The OMPS-LP products contain a large amount of data per day, whereas SAGE III and OSIRIS provide much less data due to their measurement geometry. This is why OSIRIS is missing a few days, and the SAGE III values (red, accounting all longitudes, and green, excluding SAA longitudes) for March are monthly averages. The black boxes represent time periods when the main part of the plume is in the SAA zone. For all the three plots, CALIOP 756 nm is converted from CALIOP 532 nm using Mie calculations (see Sect. 2.6). Data from OMPS-LP NASA are shown after applying filtering (see Appendix E).

Figure 6 Theoretical LR calculated using a Mie code (see Sect. 2.6) as a function of the effective radius r_eff and the mode width σ of a lognormal distribution. The black cross represents the mean values from SAGE III of r_eff and σ – respectively 0.40 µm and 1.29 – inside the plume in March 2022 between 30 and 10° S taken from Fig. 2 in Duchamp et al. (2023).

Assuming the Hunga plume consists of small spherical particles, the observed LR is liable to be compared with Mie theoretical calculations. Figure 6 shows the numerically calculated 532 nm LR as a function of PSD parameters for a lognormal distribution. LR is very sensitive to PSD variations for r_eff ranging from ≈0.1 to 0.5 µm and for low σ values. The black cross, representing the average r_eff and σ within the plume observed by SAGE III for March 2022 in the tropics, corresponds to a theoretical LR of ≈52 sr, well within the range of values determined by our CALIOP estimate, enhancing the overall confidence in our results.

4 Conclusions

The CALIOP lidar monitored the formation and evolution of Hunga aerosol plumes in the stratosphere from the day of the eruption to the end of its mission in mid-2023, fuelling numerous studies (Sellitto et al., 2022; Legras et al., 2022; Khaykin et al., 2022; Bourassa et al., 2023; Boichu et al., 2023; Bruckert et al., 2025). The quality and high resolution of its night-time measurements, with a good signal-to-noise ratio, provides several insights into the composition of the plume. First, the detection of highly depolarized pockets in the early plumes suggests the presence of remaining ice crystals or ash. The optical properties of aerosols for well-separated plumes in the first few days after the eruption, as well as for the uniform, dispersed plume in the following weeks, could then be estimated. The AOD of the plume is between 0.5 and 1.25 on the first four days, then, values decrease and stabilize at 0.047 ± 0.011 for March 2022. The LR values initially range between 60 and 80 sr, consistent with the early growth of sulfate aerosol particles as reported in both measurements and modeling studies and supported by Mie calculations, before decreasing to 48 ± 6 sr over the 27 January–25 March period. The data are compared and validated on individual profiles with direct extinction measurements by solar occultations from SAGE III when it detects the plume, and with Mie theoretical calculations. CALIOP values are higher than those of the various OMPS-LP products. In terms of daily averaged AOD, the NASA data show the closest match during the early, concentrated stage of the plume (until March 2022). For individual extinction profiles, both IUP and NASA v2.5 products provide slightly better agreement than the other OMPS-LP products.

Appendix A: NASA Standard CALIOP L2 Product

Figure A1 3 CALIOP night orbits on 17 January (a, C2), 19 (b, C1) and 27 (c). For each panel, from left to right, are represented the 532 nm extinction coefficient, the 532 nm LR (extinction/backscatter), whether the LR is constrained (1) or not (0) to retrieve the extinction and the 532 nm stratospheric AOD (SAOD).

NASA has developed standard L2 level products derived from the attenuated backscatter and offers orbital files of type “05kmAPro” (https://doi.org/10.5067/CALIOP/CALIPSO/CAL_LID_
_L2_05kmAPro-Standard-V4-51) and “05kmALay” (https://doi.org/10.5067/CALIOP/CALIPSO/CAL_LID_L2_
05kmALay-Standard-V4-51) containing variables we retrieve as part of our study such as the 532 nm extinction coefficient, the 532 nm stratospheric AOD (SAOD), the initial and final 532 nm LR and associated quality flags, among others (Young et al., 2008).

Figure A1 presents NASA L2 data for 3 orbits containing the Hunga SA plume on 17 (a, C2), 19 (b, C1) and 27 January (c). Most of the time, the LR is not constrained and its value does not change from that set at the outset, which is 50 sr. When the value is constrained, the algorithm seems to produce abnormally low LR values under 50 sr and sometimes close to 0 sr. In comparison, the values we obtain in the plume are respectively 72.5 ± 3.5, 62.6 ± 4.6 and 54.7 ± 12.4 sr for orbits from panel (a–c). This can lead to important differences for the retrieval of extinction and SAOD. For instance, mean NASA SAOD is around 0.4 (maximum around 0.7) while we find 1.01 in averaging inside the plume for orbit from panel (a).

Appendix B: CALIOP Ground Tracks on Himawari-8 and MSG-1 RGB-Ash Maps

Figure B1 exhibits the location of the CALIOP tracks shown in Fig. 2 over the aerosol plume of the Hunga as detected from geostationary imagers. The sulfur plume is visible as a greenish cloud. On Fig. B1a, the sulfur plume is just emerging from the collapsing umbrella west of the Hunga and the high altitude aerosol cloud sampled by CALIOP in Fig. 2a, which has traveled faster, cannot be seen because it is too thin. One day later, Fig. B1b shows the well separated C1 and C2 plumes, the first being the highest one in Fig. 2b. On 17 January, the plumes have separated due to wind shear and they are crossed separately on two orbits in Fig. B1c and d. The infrared images are somewhat polluted by the strong emission of Australia ground. On 19 January, the plumes which have now elongated over the Indian ocean, but can still be distinguished, are again crossed successively by two orbit tracks in Fig. B1e and f.

Figure B1 Composite RGB image from infrared channels of Himawari-8 and MSG-1 following the Eumetsat Ash Recipe (Legras et al., 2022). The two satellites are patched at 90° E. Each panel (a–f) corresponds to the same panel of Fig. 2 and the displayed CALIOP orbit track is shown in red. The timing of the CALIOP section and the RGB image differs by less than 10 min. The location of the Hunga volcano is indicated on panels (a), (b) by a red square.

Appendix C: Error Calculation of an Averaged Variable

To calculate the error bars for an averaged variable, on a daily basis or over several weeks, we use the following formula, which applies when the data are independent:

$$\begin{array}{}\text{(C1)}& {\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{X}}={\displaystyle \frac{\mathrm{1}}{\sqrt{N}}}\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\mathit{\sigma }}_i^{\mathrm{2}}}\end{array}$$

with ${\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{X}}$ the standard error of the variable $\stackrel{\mathrm{‾}}{X}$, σ_i the standard deviation associated with each element i of $\stackrel{\mathrm{‾}}{X}$ and N the number of elements in $\stackrel{\mathrm{‾}}{X}$.

Figure C1 shows the error bars for the averaged 532 nm CALIOP retrieved AOD. The error bars can be large in the first few days because we do not take into account orbits with no plume (no retrieval possible, so no associated error) in the error calculation so as not to bias the results, whereas we do take them into account as a zero value in the calculation of the mean in order to make the best possible comparison with the other instruments. Between late January and early March, errors are between 15 % and 50 %. After, errors increase again as a result of signal attenuation and noise.

Figure C1 Same as Fig. 5c with error bars for CALIOP 532 nm daily AOD.

Appendix D: Other Coincident Cases

This section presents additional cases of collocations between OMPS-LP, SAGE III, and CALIOP, which further support the conclusions drawn from the cases shown in the main text (Fig. 5a and b).

Figure D1 illustrates cases associated with rare events that allow for good spatial and temporal agreement between the spaceborne instruments, based on profiles acquired on 15 February and 20 March, in the vicinity of the collocations shown in Fig. 5a and b. These collocations are not as stringent than those presented previously; however, this slight spatial and temporal mismatch does not affect the validity of the conclusions, which remain consistent with those presented earlier.

Figure D1 Collocated profiles from CALIOP, SAGE III and OMPS-LP. (a) 15 February: CALIOP profile is averaged between 3.4 and 1.4° S at 154.0° E and 16:21 UTC. SAGE III profile is measured at 2.4° S, 155.3° E and 19:50 UTC. OMPS-LP profile is averaged between 3.4 and 1.4° S at 162.7° E and 02:36 UTC. (b) 15 February: CALIOP profile is averaged between 3.0 and 1.0° S at 129.5° E and 18:00 UTC. SAGE III profile is measured at 1.7° S, 132.2° E and 21:23 UTC. OMPS-LP profile is averaged between 3.0 and 1.0° S at 137.2° E and 04:17 UTC. (c) 20 March: CALIOP profile is averaged between 20.5 and 18.5° S at 136.6° E and 17:23 UTC. SAGE III profile is measured at 19.5° S, 134.6° E and 21:09 UTC. OMPS-LP profile is averaged between 20.5 and 18.5° S at 146.0° E and 03:53 UTC.

Appendix E: OMPS-LP NASA Data Filtering

The NASA OMPS-LP product provides both filtered and unfiltered extinction profiles (respectively RetrievedExtCoeff and RetrievedExtCoeff_NOFILT in the data files). The NASA algorithm uses the cloud-detection approach of Chen et al. (2016) to identify cloud-top height and filter out scenes containing normal clouds or polar stratospheric clouds. In version v2.5, this scheme was updated to enable cloud detection in the presence of enhanced aerosol layers (Taha et al., 2025), allowing, when possible, the separate identification of elevated aerosol layers and underlying clouds. This refinement may lead to a slight improvement in extinction measurements in the upper troposphere-lower stratosphere by enabling the removal of clouds located below volcanic aerosol layers. In addition, version v2.5 introduces two supplementary quality-control filters, whereby extinction values are excluded when the retrieved aerosol extinction exceeds 0.1 km⁻¹ or when the retrieval residuals are greater than 50 %, owing to the high likelihood that such retrievals are unreliable.

Figure E1 presents the differences between AOD derived from filtered and non-filtered extinction for both v2.1 and v2.5. In addition to the stronger extinction values and the larger variability already mentioned in the main text, much larger differences are observed between filtered and unfiltered AOD values for v2.5 than for v2.1. In v2.1, the cloud filter primarily identifies isolated cloud structures, whose impact above 17.5 km – the lower altitude limit used in the AOD calculation – is relatively small. In contrast, the v2.5 filtering scheme also flags as clouds some structures that may be located beneath an aerosol plume.

For version v2.5, we illustrate the impact of applying the new filtering scheme on the AOD. The unfiltered data (“raw data”) exhibit very large values compared to the other instruments considered in this study. An initial removal of extinction data exceeding 0.1 km⁻¹ results in a first significant reduction in the AOD values in many cases. Subsequently, the combined application of the other filters (cloud screening and removal of values with residuals greater than 50 %) brings the AOD estimates into overall agreement with the range of values reported by CALIOP and SAGE III (Fig. 5c). The application of such filters leads in a number of cases to a noticeable reduction in the derived AOD by masking extinction values that are likely overestimated, particularly below the Hunga plume. We have only found two orbits where the filters mask likely valid data. This is why we use the filtered extinction for the OMPS-LP NASA AOD calculation. It should be noted that the differences may vary depending on the wavelength considered, as the extinction coefficient generally decreases with increasing wavelength. We decided to present the contribution of the filters separately in this form for v2.5, also due to a practical issue. The input files stipulate a valid range of 0 to 0.1 km⁻¹ for the extinction, so that many client readers (e.g. python libraries) automatically mask data in excess or treat them as missing. The raw data can only be seen after unmasking. Therefore the first step in filtering can be easily missed although it was not intended by the provider. For v2.1, the situation is similar; however, in practice, extinction values never, or only very rarely, exceed 0.1 km⁻¹, so the choice of the data processing approach has no impact on the resulting values of the RetrievedExtCoeff_NOFILT quantity.

Figure E1 Same as Fig. 5c for OMPS-LP NASA v2.1 (purple) and v2.5 (light blue shades) 745 nm daily AOD derived from filtered and non-filtered extinction. The introduction of additional filters in v2.5 allows us to include a so-called semi-filtered curve, in which extinction values exceeding 0.1 km⁻¹ are filtered out.