Radiative impact of the Hunga stratospheric volcanic plume: role of aerosols and water vapor over Réunion Island (21° S, 55° E)

This study attempts to quantify the radiative impact over Réunion Island (21° S, 55° E) in the southern tropical Indian Ocean of the aerosols and water vapor (WV) injected into the stratosphere by the eruption of the Hunga underwater volcano in the South Pacific on 15 January 2022 . Ground-based lidar and satellite passive instruments are used to parameterize a state-of-the-art radiative transfer (RT) model for the first 13 months after the volcano eruption. The descending rate of the aerosol volcanic plume is −8 m d⁻¹. At this rate, aerosols are expected to be present in the stratosphere until the first half of 2025. The overall aerosol and water vapor impact on the Earth's radiation budget for the whole period is negative (cooling, −0.82 ± 0.35 W m⁻²) and dominated by the aerosol impact (∼ 95 %; the remaining ∼ 5 % is due to the water vapor). At the Earth's surface, aerosols are the main drivers and produce a negative (cooling, −1.04 ± 0.36 W m⁻²) radiative impact. Water vapor has hardly any radiative effect at the surface. Between the short-term (months 2 to 4 after the eruption, February–April 2022) and mid-term (months 5 to 14 after the eruption, May 2022–February 2023) periods, the aerosol and water vapor radiative effect at the surface and top of atmosphere (TOA) reduces by 22 % and 25 %, respectively. During the mid-term period, heating $/$ cooling (H / C) rate profiles show a clear vertical difference locally in the stratosphere between the aerosol warming impact (18 to 26 km) and the water vapor cooling (22 to 30 km). The resulting aerosol and water vapor heating $/$ cooling rate profile follows an S-shaped curve with peaks slightly larger for the moist layer (−0.09 K d⁻¹) than for the sulfate layer (+0.06 K d⁻¹).

1 Introduction

More than 1.5 years after the eruption of the Hunga underwater volcano in the South Pacific, the scientific community is still actively investigating the climate impact of the huge amounts of water, steam, and gas that were injected into the atmosphere. The event showed a fast spatiotemporal global dispersion of the stratospheric volcanic matter that circulated the Earth in only 1 week (Khaykin et al., 2022), with small parts of the main aerosol layer dispersed pole-to-pole in 3 months (Taha et al., 2022), first in the form of concentrated patches (Legras et al., 2022). Several features show evidence of a record-breaking atmospheric event. The eruption, equivalent to an energy of 110 Tg of trinitrotoluene (TNT) explosive, is the most powerful volcanic explosion since Krakatau (1883) and Tambora (1815) (Lac et al., 2022). The volcanic plume reached an altitude of 57 km, a coincident estimation resulting from different techniques (Carr et al., 2022; Proud et al., 2022), placing it in the upper stratosphere–lower mesosphere, a record in the satellite era. The mass of water retained in the atmosphere was unprecedented, estimated by Millán et al. (2022) to be 146 Tg injected into the stratosphere (e.g., 37 Tg of water was retained in the stratosphere as a result of cross-tropopause transport after the 1991 Pinatubo eruption; Pitari and Mancini, 2002). In contrast, the mass injection of sulfur dioxide (SO₂) was not that exceptional: ∼ 0.6–0.7 Tg (Carn et al., 2022), which is much smaller than that from previous major eruptions (e.g., 20 Tg for Pinatubo; Bluth et al., 1992). Still, the stratospheric aerosol optical depth (sAOD) has been recorded globally as the largest since the Pinatubo eruption (Taha et al., 2022).

The latter three variables (water, SO₂, and injection height), and co-emitted halogens to a lesser extent, are some of the main factors responsible for the production of volcanic sulfate and for the loss/production of ozone. The initial SO₂ was fully converted into sulfates in less than 2 weeks under the influence of water vapor (WV) (Asher et al., 2023; Legras et al., 2022), whereas volcanic sulfate and water still persist as of today. The fast water vapor injection provided abundant hydroxide (OH), which reacted with SO₂ to form volcanic sulfate at a faster rate than the typical ∼ 30 d (Carn et al., 2016). Higher concentrations of volcanic sulfate led to more rapid coagulation and thus to particles quickly growing in size. In the case of the Hunga volcano, this mechanism is estimated to have halved the SO₂ lifetime and doubled the sAOD (Zhu et al., 2022). This rapid growth and the global persistence of volcanic sulfate aerosols have been demonstrated with AERONET measurements by Boichu et al. (2023), with the occurrence of an unusual “volcanic fine mode” with a peak radius ranging 0.28–0.50 µm. This fine mode was found to be poorly absorbing, although balloon-borne measurements by Kloss et al. (2022) report a moderately absorbing fine mode in the first 10 d after the eruption, indicating small sulfate-coated ash particles. Volcanic sulfate is known to impact ozone depletion by providing additional surface area and suppressing the nitric oxide cycle (Tie and Brasseur, 1995). The transport of volcanic sulfate from the tropics to the Antarctic by the Brewer–Dobson circulation contributed to an increase in ozone concentrations in the middle stratosphere but to a decline in the lower stratosphere at mid- to low latitudes (Lu et al., 2023), while, combined with a cold polar vortex, it contributed to a decreased ozone concentration in the Antarctic (Wang et al., 2023). Because ozone is not primarily emitted during volcanic eruptions, its loss or production by post-eruption reactions is more tedious to estimate (Evan et al., 2023). The effect of the Hunga volcano on stratospheric ozone is still under study.

Water, volcanic sulfate, and the injection height are the main drivers of the impact of the Hunga volcano on atmospheric global circulation (Coy et al., 2022) and climate (Zuo et al., 2022). In particular, the climate forcing will depend on the radiative effect produced by the water vapor longwave (LW) emission and absorption and by the sulfate aerosol longwave and shortwave (SW) near-infrared absorption (Robock, 2000). These interaction mechanisms (emission and absorption) with the longwave and shortwave near-infrared radiation are height-dependent and determine the sign of the differential of energy gained (positive) or lost (negative) in all layers of the atmosphere. Several studies have demonstrated the stratospheric cooling produced by the excess of water vapor injected by the Hunga volcano either locally (Sellitto et al., 2022), zonally (Schoeberl et al., 2022; Vömel et al., 2022; Zhu et al., 2022), or globally (Millán et al., 2022; Schoeberl et al., 2024) at different timescales spanning from instantaneous estimates to 6-month evolutions. As far as volcanic sulfates are concerned, these aerosols usually scatter sunlight back to space, cooling the Earth's surface, and absorb outgoing thermal radiation. Several authors have made the hypothesis that the Hunga volcano eruption could impact climate not through surface cooling due to sulfate aerosols, but rather through surface warming due to the radiative forcing from the excess stratospheric water vapor. To date, assessments of the radiative effect of combined water vapor and aerosols have only been performed for three case studies: during the first 10 d after the eruption by Sellitto et al. (2022), for the first 2 months after the eruption by Zhu et al. (2022), and for the first year (Gupta et al., 2025) and the first 2 years (Schoeberl et al., 2024) after the eruption. Jenkins et al. (2023) evaluated the chances of a temporary global surface temperature anomaly above 1.5 °C over the coming decade caused by the Hunga volcano stratospheric water vapor perturbation.

Here, the impact of water vapor and aerosols on the Earth's radiation budget is estimated over Réunion Island (21° S, 55° E) for the first 13 months after Hunga volcano eruption. Both water vapor and aerosols obtained from ground-based lidar and satellite measurements are used as input in a state-of-the-art radiative transfer (RT) model. The radiative effect is calculated for three scenarios considering aerosols only, water vapor only, and combined aerosols and water vapor.

2 Materials and methods

2.1 The Maïdo instrumentation

The lidar system used in this study, the Li1200 lidar (Dionisi et al., 2015; Vérèmes et al., 2019; Gantois et al., 2024), is located at the Observatoire de Physique de l'Atmosphère à La Réunion (OPAR) Maïdo station (21.079° S, 55.383° E; 2160 m a.s.l.; Baray et al., 2013). The system operates in routine mode at 355 nm, and measurements are made twice a week on Monday and Tuesday nights. In this work, 87 nights of observations were recorded between 19 January 2022 and 15 February 2023. A full description of the system is available in the data paper of Gantois et al. (2024).

The extinction presented in this work is obtained by applying the elastic two-component inversion algorithm (Klett, 1985) using a constant lidar ratio (LR). Several LR values were tested between 30 and 70 sr. The value of 30 sr was fixed for this study. The transmittance method initially used in Baron et al. (2023) for the thick plume observed during the first days after the eruption over Réunion Island was not retained because, with the decreasing aerosol load observed in our longer-term study, the transmittance method would have led to large uncertainties in the LR retrieval at such altitude levels (17–32 km). Although rather unusual for sulfate aerosols, which are more often associated to an LR of 60 sr according to the existing literature (e.g., Lopes et al., 2019), the value of 30 sr is chosen following the results presented in Baron et al. (2023). Indeed the latter found values of LR at 355 nm in the range 29–35 sr with small standard deviations (< 7 sr) by applying the transmittance method during several nights in January 2022. These uncommon lidar ratio values for sulfate aerosols were proved to be stable over time by Duchamp et al. (2023) using Stratospheric Aerosol and Gas Experiment (SAGE-III) observations. The use of 30 sr at 355 nm is also coherent with Mie calculations conducted by Baron et al. (2023) and with the size distribution parameters from Duchamp et al. (2023).

The uncertainty associated to the extinction profiles at 355 nm, and, by extension, to the sAOD at 355 nm, has been calculated considering an uncertainty in the lidar ratio of ±10 sr. This value of ±10 sr corresponds to the largest uncertainty calculated on the lidar ratio at 355 nm by Baron et al. (2023) for the Hunga plume over Réunion Island in January 2022. To compare the lidar-derived and Ozone Mapping and Profiler Suite (OMPS; see Sect. 2.2) sAOD, the wavelength of 745 nm was used (see Sect. 2.2). sAOD at 355 nm was converted into sAOD at 745 nm using a constant Ångström power law, ${\mathrm{AE}}_{\mathrm{355}/\mathrm{745}}$, of −0.14 resulting from our Mie calculation (see Sect. 2.3). The uncertainty associated to the sAOD at 745 nm has been calculated considering both the uncertainty associated to sAOD at 355 nm and the uncertainty associated to ${\mathrm{AE}}_{\mathrm{355}/\mathrm{745}}$, fixed to a constant value of ±1.0 (Baron et al., 2023).

2.2 Satellite and reanalysis data

The Ozone Mapping and Profiler Suite (OMPS) Limb Profiler has been on board the Suomi National Polar Partnership (NPP) since October 2011. Using limb-scattering solar radiation, OMPS provides good-quality aerosols extinction retrievals at several wavelengths: 510, 600, 675, 745, 869, and 997 nm (Taha et al., 2021). As recommended by the latter, we use data product version 2.0 of the aerosol extinction profile at 745 nm to follow the aerosol volcanic plume over Réunion Island from January 2022 to mid-April 2023. These data are provided from 10 to 40 km height on a vertical grid of 1 km. Stratospheric aerosol optical depth calculations are made by integrating the extinction profiles from 17 to 40 km, where 17 km corresponds to the tropopause height over Réunion Island (Bègue et al., 2010). Based on previous studies in the Southern Hemisphere (Bègue et al., 2017; Tidiga et al., 2022), background periods extend from 2012 to February 2014 and from January 2017 to April 2018 to exclude volcanic eruptions (Kelud, Calbuco, Ambae, and Ulawun) and the Australian 2019/2020 biomass burning episode (Black Summer).

Microwave Limb Sounder (MLS) version 5.0 level-3 data are also used to extract the monthly mean water vapor over our site and in the stratosphere during 2021 to serve as a climatological reference. See https://mls.jpl.nasa.gov/data/v5-0_data_quality_document.pdf (last access: 10 January 2025) for more details about this MLS product. The monthly mean of the water vapor in the altitude range of interest in 2021 is 4.5 ppmv. This value sets the climatological reference necessary to parameterize the unperturbed conditions of the water vapor.

The MERRA-2 Stratospheric Composition Reanalysis of Aura MLS (M2-SCREAM) products are used for characterizing the water vapor (WV) and ozone vertical distribution, in particular the 3D 3-hourly GMAO_M2SCREAM_INST3_CHEM product available until 1 April 2023 (https://doi.org/10.5067/7PR3XRD6Q3NQ, GMAO, 2022​​​​​​​). This product, produced at NASA's Global Modeling and Assimilation Office (GMAO), is generated by assimilating MLS and Ozone Monitoring Instrument (OMI) retrievals into the Goddard Earth Observing System (GEOS) Constituent Data Assimilation System (CoDAS) driven by meteorological fields from MERRA-2. Stratospheric water vapor and ozone, among other compounds, are assimilated in M2-SCREAM. Assimilated fields are provided globally at a resolution of 0.5° (latitude) by 0.625° (longitude) from approximately 10 km up to the lower thermosphere. Concretely, the variables of specific humidity (QV, kg kg⁻¹), ozone (O₃, ppmv), mid-layer pressure (PL, Pa), and mid-layer height (H, m) were used. The specific humidity was converted to the actual water vapor pressure and then to the water vapor mixing ratio. All variables were averaged over 4 pixels surrounding the Maïdo coordinates. Assimilation uncertainties for each of the assimilated constituents are calculated from the CoDAS statistical output (Wargan et al., 2023). For the period January 2022 to September 2022, and in the height interval of interest in this study (17–32 km), the uncertainties on the water vapor and ozone are less than 0.2 and 0.13 ppmv (Wargan et al., 2023), respectively.

2.3 The GAME radiative transfer model

2.3.1 Code and parameterization

Radiative fluxes propagating through the atmosphere were calculated with the radiative transfer (RT) model GAME (Dubuisson et al., 1996; Dubuisson, 2004; Dubuisson et al., 2006). For this study, GAME was set up to calculate spectrally integrated upward and downward radiative fluxes in 40 plane and homogeneous layers from 0 to 100 km with a 1 km resolution from 0 to 30 km and a coarser resolution above. The shortwave (SW) spectral range was set from 0.2 to 4.0 µm (wave number resolution of 400 cm⁻¹ from 0.2 to 0.7 µm and 100 cm⁻¹ from 0.7 to 4.0 µm). In the longwave (LW) spectral range, spectral limits were defined between 4.0 and 50.0 µm (115 points at a wave number resolution of 20 cm⁻¹). GAME calculates solar flux values at the boundary of plane and homogenous atmospheric layers by using the discrete ordinates method (Stamnes et al., 1988). Gas (H₂O, CO₂, O₃, N₂O, CO, CH₄, and N₂ are considered) absorption is calculated from the correlated k distribution (Lacis and Oinas, 1991). More details about the computation of the gas transmission functions can be found in Dubuisson (2004) and Sicard et al. (2014). In the longwave spectral range, GAME presents the advantage of the complete representation of the longwave aerosol scattering, in addition to their absorption (Sicard et al., 2014).

For the sake of clarity and comparability with other works, we recall the definition of the direct radiative effect (DRE) of a perturbed vs. unperturbed atmospheric compound on the Earth's radiation budget. At a given height level, L,

$$\begin{array}{}\text{(1)}& \mathrm{DRE}\left(L\right)=\left[F_{\mathrm{p}}^{\downarrow }\left(L\right)-F_{\mathrm{p}}^{\uparrow }\left(L\right)\right]-\left[F_{\mathrm{u}}^{\downarrow }\left(L\right)-F_{\mathrm{u}}^{\uparrow }\left(L\right)\right],\end{array}$$

where F is the radiative flux values for perturbed (p subindex) and unperturbed (u subindex), while the ↓ and ↑ arrows indicate, respectively, the downward and upward flux direction. By that definition, negative (positive) DRE values represent a cooling (warming) effect. The DRE was calculated at two climate-relevant altitude levels: at the top of atmosphere (TOA) and at the bottom of atmosphere (BOA).

As far as GAME parameterization is concerned, temperature and pressure profiles used in both SW and LW simulations are taken from radiosoundings launched from Saint-Denis, the state capital of Réunion Island, 20 km north of Maïdo, every night at 00:00 local time. Aerosols are fully parameterized in GAME by the user in terms of spectrally and layer-mean aerosol optical depth (AOD), single-scattering albedo (SSA), and asymmetry factor (asyF). The layer-mean AOD is distributed vertically according to the profiles of the measured extinction coefficient at 355 nm, whereas SSA and asyF are assumed vertically constant.

The spectral AOD, SSA, and asyF were calculated in the whole spectral range with a Mie code. A monomodal lognormal size distribution was considered, with a geometric median radius of 0.35 µm and a mode width of 1.23. These values, taken from Duchamp et al. (2023), are the retrieved particle size distribution of SAGE-III in the 30–10° S latitude range and corresponding to plume conditions at altitude of maximum extinction averaged between the months of June and August 2022. SAGE-III observations also show that this size distribution at the plume peak height persisted over 17 months with only a small decreasing trend in the size. For the refractive index we used the Gestion et Étude des Informations Spectroscopiques Atmosphériques: Management and Study of Spectroscopic Information (GEISA) spectroscopic database (Jacquinet-Husson et al., 2008). In particular, the refractive index of the binary system H₂SO₄ $/$ H₂O with an H₂SO₄ mixing ratio (in mass, w_t, i.e., the ratio of the H₂SO₄ mass to the total mass of the droplets) of 0.75 and a temperature of 215 K (this temperature corresponds on average to the atmospheric temperature at the height of the volcanic plume) was used. This value was selected in view of the results of Duchamp et al. (2023), who, in their supporting information, show w_t profiles retrieved from zonal average profiles of water vapor (retrieved from MLS, version 5) and temperature (from the ERA5 reanalysis) at latitudes 0, 10, 20, and 30° S for a set of dates between February 2022 and April 2023. At our latitude (20° S), w_t ranges between 0.70 and 0.80 between 22 and 28 km (where the moist layer is located at least until November 2022) for all dates shown. Above 28 km, w_t > 0.80. The real part is defined over the range 0.61–5000.00 µm (wave number resolution of 2 cm⁻¹). The value at 0.61 µm was assumed constant in the range 0.20–0.61 µm. The imaginary part is defined over the range 2.36–23.15 µm (wave number resolution of 0.96 cm⁻¹). The value at 2.36 µm was assumed constant in the range 0.20–2.36 µm. The reader is referred to Biermann et al. (2000) for more details on this dataset. Figure 1 shows the real part (RRI) and imaginary part (IRI) of the refractive index used for w_t = 0.75. The curves for w_t = 0.80 are shown only for comparison. Large spectral variations in the infrared atmospheric window (8–13 µm), which have an important impact on the infrared radiative budget of the atmosphere, are visible. The most astonishing feature of the figure is probably the high absolute values of the IRI, emphasizing the high absorbing properties of sulfate aerosols in the longwave spectral range which induce extremely low SSA (< 0.2) in the whole range. One sees that, between w_t = 0.75 and w_t = 0.80, nearly no difference is visible, so employing w_t = 0.75 at all altitudes (i.e., even above 28 km) should not impact the LW radiative properties. For comparison, the IRI (this study) is 2 to 3 times larger in the atmospheric window than the IRI for mineral dust (Sicard et al., 2014).

Figure 1 Spectral complex refractive index considered for the calculation of the aerosol radiative properties in the shortwave (0.2–4 µm) and longwave (4–50 µm) spectral ranges. See text for details. The infrared atmospheric window (8–13 µm) is indicated by the shaded blue area. The dashed lines are the extrapolation of the dataset used.

The spectral SW surface albedo was interpolated from the four-wavelength AERONET L2.0 annual mean for the year 2022 and assumed constant above 1020 nm. The AERONET site used is “Maido_OPAR”. It is situated at 2160 m a.s.l. and is co-located to the lidar. For the LW broadband surface albedo, we used the Infrared Atmospheric Sounding Interferometer (IASI) December nighttime monthly mean climatology of the surface emissivity (i.e., surface albedo of 1) at 890 cm⁻¹, i.e., 11.24 µm, over the Indian Ocean (Zhou et al., 2013) and set the surface albedo value to 0.01.

To avoid the dependency on solar zenith angle, only daily radiative effects are presented in this work. To do so, each nighttime measurement and parameterization is assumed to be constant for the 24 h of the day considered, and both SW and LW radiative effects are calculated at an hourly time resolution between 00:00 and 23:00 UT. In these calculations, the solar zenith angle is the only parameter that varies. The daily radiative effect is the average of the 24-hourly DRE.

2.3.2 Error estimate

An error budget is performed to quantify the uncertainties on the radiative effect estimations using GAME and caused by the model itself, our parameterization, and the hypothesis made. The GAME model was used in an intercomparison exercise (Halthore et al., 2005) which concluded that it is accurate to within a few watts (< 5 W) for a flux reaching 1000 W m⁻². The impact of this uncertainty on our estimations should be even less, since only daily averaged fluxes are considered. It is thus reasonable to consider an uncertainty in relative terms of 0.5 %.

Two other sources of error are considered: one associated to the lidar ratio selected and another associated to the size distribution selected. New profiles of the extinction inverted using (30 + 10) sr and (30 − 10) sr were used in GAME to quantify the deviation from the nominal (LR = 30 sr) radiative effect estimations. As far as the size distribution is concerned, Duchamp et al. (2023) detected “a small decreasing trend in the size” without quantifying it. We have assumed a decrease in the geometric median radius of −0.01 µm. Thus, a new Mie calculation was performed with a geometric median radius of 0.34 µm, and the resulting radiative properties were used in GAME to quantify the deviation from the nominal (geometric median radius of 0.35 µm) radiative effect estimations. The results from these uncertainties are given in Table 1 in relative terms at BOA and TOA. Logically, the lidar ratio error, which impacts the sAOD error almost proportionally, is by far the largest. We can reasonably consider that the aerosol daily radiative effects are estimated with an uncertainty better than 48 % at TOA and better than 42 % at BOA.

Table 1 Error estimate of the aerosol daily radiative effect.

3 Vertical/temporal evolution of the Hunga volcanic plume over Réunion Island

The historical context of the aerosol load over Réunion Island is shown in Fig. 2 by the temporal evolution of the stratospheric AOD at 745 nm measured by OMPS in the last decade. The background sAOD is measured over the unperturbed years 2012 and 2013. It is (2.6 ± 0.1) × 10⁻³. At each exceptional event the sAOD takes off from this background sAOD, and, since the eruption of Ambae in July 2018, the sAOD over Réunion Island has never turned back to its background value. The sAOD peak produced by the Hunga volcano (0.035) is the highest in the last decade, and it is 4 times higher than the second-highest event (0.009; Calbuco eruption in April 2015). Zonal averages between 30° S and 15° N for the Hunga volcano and 20° S and 90° S for Calbuco showed that the Hunga sAOD was more than double that for the 2015 Calbuco eruption (Taha et al., 2022). The reason why the local and zonal sAOD differences between the Hunga volcano and Calbuco differ lies in the zonal mean stratospheric conditions. In the case of the Hunga volcano, a marked easterly band (Khaykin et al., 2022; Legras et al., 2022) favored direct transport from the Hunga volcano towards Réunion Island (both being at approximately the same latitude). Further back historically, the 40+ year satellite record of monthly sAOD at the scale of the globe (i.e., for the 60° S–60° N latitude band) in Khaykin et al. (2022) shows that only the eruptions of Pinatubo (1991) and El Chichón (1982) exceeded the Hunga one in terms of absolute stratospheric AOD (by a factor of 6 and 3, respectively).

Figure 2 sAOD (17–40 km) at 745 nm from OMPS over Réunion Island. The most important volcanic eruptions (name and date) in the Southern Hemisphere are indicated by vertical red lines. The Australian 2019/2020 biomass burning episode is indicated in orange.

In Fig. 3, the vertical and temporal (January 2022–April 2023) evolution of the Hunga volcanic plume over Réunion Island is analyzed by means of sAOD at 745 nm and profiles of extinction coefficient at 745 nm, water vapor, and ozone. While the monthly OMPS sAOD peak is reached in April and May 2022, the instantaneous lidar sAOD peaks just a few days after the eruption, reaching 0.40 on 21 January. This time difference is an indication of the dispersion time of the volcanic matter injected by the Hunga volcano in the stratosphere in the zonal direction. Other studies confirm that some parts of the volcanic plume dispersed pole to pole in 3 months (Khaykin et al., 2022; Taha et al., 2022). Another indicator of this dispersion in the zonal direction is the standard deviation (calculated as a 15 d rolling standard deviation) associated to OMPS monthly sAOD: once past the first month, it steadily decreases all along the year 2022. The agreement between monthly and instantaneous sAOD, which becomes excellent as of April 2022, is also an indicator of the homogenous dispersion of the volcanic plume in the stratosphere at our latitudes. A decrease in the monthly sAOD is observed after April/May and until November. Then sAOD stabilizes until today (sAOD = 0.012, almost 5 times the background sAOD). From September 2022 on, the lidar sAOD is slightly higher than the OMPS sAOD, although the error bars always overlap one another. It is not clear whether this reflects a systematic difference between both estimations or a limitation of one of the two datasets.

Figure 3 (a) Daily (OMPS) and instantaneous (lidar, nighttime) sAOD (17–40 km) at 745 nm over Réunion Island. Time–height plots of (b) extinction coefficient at 745 nm from OMPS, (c) water vapor, and (d) ozone mixing ratio from M2-SCREAM. The shaded green area in panel (a) is the sum of the standard deviation (calculated as a 15 d rolling standard deviation) and the 10 % relative accuracy of OMPS at 745 nm according to Taha et al. (2021). The red and black lines in panels (b) and (c) report the peak height of the aerosol and water vapor plumes, respectively.

The time–height plot of the extinction (Fig. 3b) clearly shows the height and vertical extension of the volcanic plume which is still present on 15 April 2023 and is located at 18.5–23.5 km height (sAOD = 0.012). The plume peak height has a decreasing tendency since April 2022 at an average steady rate of −244 m per month or ∼ −8 m d⁻¹. Assuming this rate constant in time and a tropopause height on Réunion Island of 17 km (Bègue et al., 2010), the remaining lifetime of the volcanic plume in the stratosphere is estimated to be between 2 and 2.5 years after 15 April 2023. Except during the first week of detection above Réunion Island, the Hunga volcanic plume is not detected above 30 km. The water vapor plume (Fig. 3c) also clearly reveals the unusually high water vapor concentration caused by the volcanic plume. A local peak of 65 ppmv is reached on 13 February 2022. It is almost 15 times higher than the climatological reference value of 4.5 ppmv (see Sect. 2.2). On a monthly basis, the water vapor stratospheric peak in February 2022 is approximately 5 times higher than the climatological reference (4.5 ppmv). This ratio decreases to almost 2 in February 2023. The water vapor plume is thinner than the aerosol one and is located at a higher altitude, 3 to 4 km higher. Such a difference, although not so accentuated, is observed zonally at 15° S during the first 6 months of the year 2022 (Schoeberl et al., 2022). The heights of the peak of the aerosol and water vapor layers (respectively, red and black lines in Fig. 3b and c) have opposite tendencies as of April 2022: the aerosol plume slowly descends, whereas the moist layer ascends slowly until October 2022 and at a higher rate afterwards. Schoeberl et al. (2022) explain that “the water vapor is transported upward with the diabatic circulation that gives rise to the tropical trace gas tape recorders (Schoeberl et al., 2018) whereas the aerosols are gravitationally settling”. Legras et al. (2022), who analyze the same period (first 6 months of the year 2022), explain that the ascent of the moist layer is due to the Brewer–Dobson circulation. The ozone cycle (Fig. 3d), with highs in the austral summer (January–April) and lows in the austral winter (July–October), reflects the higher production of ozone in summer due to the peak of solar radiation compared to winter (Abdoulwahab, 2016). Apart from this natural cycle of stratospheric ozone at subtropical latitudes, no other spatiotemporal variation, potentially caused by the Hunga volcano eruption, is visible to the naked eye in Fig. 3d. Some authors mention that, following the Hunga volcano eruption, the ozone concentration increase in the middle stratosphere and decrease in the lower stratosphere were caused by enhanced sulfate aerosol (Lu et al., 2023), and others claim that the midlatitude and tropic ozone reductions observed by MLS were mainly linked to circulation effects (Wang et al., 2023). Above the Indian Ocean, Millet et al. (2024) reported a transient ozone-depleted area during the first week after the eruption. At this early stage of our understanding of the effects of the Hunga volcano on the stratospheric ozone, the present study does not consider any potential increase/decrease in stratospheric ozone due to the Hunga volcano eruption.

4 Impact of the Hunga volcanic plume on the Earth's radiation budget

In order to analyze the radiative impact of the aerosols and the water vapor, separately and all together, four parameterizations of GAME are performed and summarized in Table 2. The perturbed condition is the full parameterization with observed sAOD and water vapor mixing ratio. For the unperturbed conditions, the impact of aerosols is assessed by assuming an aerosol-free stratosphere, the impact of water vapor is assessed by assuming that the water vapor mixing ratio in the Hunga moist layer is equal to the climatological value of 4.5 ppmv obtained from the MLS 2021 monthly means, and the impact of aerosols and water vapor is assessed by assuming both an aerosol-free stratosphere and a water vapor mixing ratio of 4.5 ppmv in the Hunga moist layer.

Table 2 Aerosol and water vapor parameterizations for the perturbed/unperturbed simulations of GAME.

Figures 4, 5, and 6 show the radiative impact of aerosols only, of the water vapor only, and of both aerosols and water vapor, respectively, in terms of time plots of DRE(BOA) and DRE(TOA) and time–height plots of the heating $/$ cooling (H $/$ C) rate anomaly in the SW, LW, and SW + LW spectral ranges. We analyze three different periods of time (see shaded gray areas in Fig. 4), excluding from now on the first 2 weeks after the eruption to allow for some dispersion to happen:

• The entire period from February 2022 to February 2023 (M2–M14, with M1 being January 2022);

• February 2022–April 2022 (M2–M4), the short-term period;

• May 2022–February 2023 (M5–M14), the mid-term period.

The first period is representative of the radiative impact of the Hunga volcano since the eruption to date, while the second and third periods are representative of the short- and mid-term tendencies, respectively. The radiative effects of the three simulations (aerosols only, WV only, aerosols and WV) associated to the three periods at the three atmospheric levels (BOA, ATM, TOA) are summarized in Table 3. Figure 7 shows the averaged H $/$ C rate anomaly profiles over both the M2–M4 and M5–M14 periods and for the three simulations (aerosols only, WV only, aerosols and WV).

Figure 4 Aerosol direct radiative impact. (a) Radiative effect (W m⁻²) at TOA; time–height evolution of the H $/$ C rate anomaly in the (b) SW, (c) LW, and (d) SW + LW spectral ranges; and (e) radiative effect (W m⁻²) at BOA. In panel (a), the shaded gray areas represent the short- and mid-term periods: February to April 2022 (M2–M4) and May 2022 to February 2023 (M5–M14).

Figure 5 Same as Fig. 4 but for the water vapor. Note that the left-axis scale in panels (a) and (e) is different to that in Fig. 4.

Figure 6 Same as Fig. 4 but for the aerosols and water vapor.

Figure 7 Mean profiles of aerosols only, water vapor only, and aerosol and water vapor daily heating $/$ cooling rate anomaly averaged over the short-term period (M2–M4) and the mid-term period (M5–M14), all excluding weeks 1 and 2 after the eruption.

Table 3 SW + LW DRE at BOA and at TOA produced by aerosols only, WV only, and both aerosols and WV. These values are the average over the entire period M2–M14, the short-term period M2–M4, and the mid-term period M5–M14, all excluding weeks 1 and 2 after the eruption.

Aerosols (Fig. 4). In all moments, the aerosol SW component dominates over the LW one at TOA (Fig. 4a), producing a negative net radiative effect (−0.78 ± 0.35 W m⁻²). The net aerosol effect at TOA is stronger during M2–M4 (−0.93 ± 0.32 W m⁻²) than during M5–M14 (−0.74 ± 0.36 W m⁻²), which is caused by the sAOD decrease. Gupta et al. (2025) estimated a net instantaneous clear-sky radiative energy loss caused by an enhanced sAOD of −0.48 ± 0.04 W m⁻² at TOA in the Southern Hemisphere. At BOA, the aerosol LW component is nearly zero, so the aerosol net DRE is that of the SW component: a cooling is observed, stronger during M2–M4 (−1.26 ± 0.50 W m⁻²) than during M5–M14 (−0.98 ± 0.32 W m⁻²). It is extremely important to contrast those estimations to the uncertainties estimated in Sect. 2.3.2, as the relative variability observed is in the same order of magnitude as this uncertainty. As a consequence of the neutral SW scattering and strong LW absorption, the time–height plots of H $/$ C rate anomaly (Fig. 4b–d) caused by the SW + LW aerosol impact is locally positive with daily heating rates (averaged over both the period considered and the altitude range 17–32 km; see Fig. 7) of +0.04 K d⁻¹ (peaking at 26 km) during M2–M4 and of +0.02 K d⁻¹ (peaking at 23 km) during M5–M14. Our mid-term (M5–M14) heating rate (+0.02 K d⁻¹) is in quite good agreement with the annual mean zonal value at 20° S averaged between 20 and 36 km, ∼ +0.02 K d⁻¹, estimated by Gupta et al. (2025). The time series of the H $/$ C rate profiles (Fig. 4d) follows that of the aerosol extinction profiles, with a decreasing tendency starting in April 2022 (Fig. 3b).

Water vapor (Fig. 5). The water vapor radiative effect is dominated by the cooling effect of water vapor longwave emission and absorption in the moist layer (Fig. 5c), although the SW warming is not negligible in the first 3 months (Fig. 5b). This layer produces a slightly negative effect at TOA of −0.07 ± 0.02 W m⁻² during M2–M4, which decreases to −0.03 ± 0.01 W m⁻² during M5–M14. Sellitto et al. (2022) also estimated a negative DRE(TOA) caused by water vapor for the fresh plume (instantaneous values of −0.7 and −0.4 W m⁻²) but a positive DRE(TOA) caused by water vapor for the aged plume (8 February 2022) of +0.8 W m⁻² and attributed to the descent in altitude of the moist layer. Our analysis supports neither this change in sign nor this direction of the vertical motion of the moist layer over the long term. Zhu et al. (2022), who use a global climate model to simulate the radiative effect in the first 2 months of 2022, find that, when only water vapor is injected in their model (without sulfur dioxide; see their supplementary material), the zonal DRE(TOA) at the latitude of Réunion Island is neutral to positive and much smaller than that caused by aerosols (their global mean over the first 2 months is −0.02 W m⁻²). The discrepancy with our findings stems from an excess of water vapor in the Zhu et al. (2022) simulation, since the reaction of sulfur dioxide (not present) and hydroxide is not happening in their simulation, in different heights of the moist layer, and/or in zonal vs. local computations. The water vapor radiative effect at BOA is negligible. The SW + LW water vapor radiative impact in the stratosphere is mostly negative, with daily cooling rates (averaged over both the period considered and the altitude range 17–32 km; see Fig. 7) of −0.07 K d⁻¹ (peaking at 26 km) during M2–M4 and −0.03 K d⁻¹ (peaking at 26 km) during M5–M14. The same ascending behavior of the water vapor concentration (Fig. 3c) is observed on the profiles of the SW + LW water vapor cooling rate (Fig. 5d). Our mid-term (M5–M14) cooling rate (−0.03 K d⁻¹) is in quite good agreement with the annual mean zonal value at 20° S averaged between 20 and 36 km, ∼ −0.015 K d⁻¹, estimated by Gupta et al. (2025). Schoeberl et al. (2022), who estimated the LW zonal impact of water vapor at 15° S for the first 6 months of the year 2022, show a cooling effect in the stratosphere with a stronger peak at ∼ −0.5 K d⁻¹ at the end of February and decreasing afterwards.

Aerosols and water vapor (Fig. 6). The overall Hunga aerosol and water vapor impact on the Earth's radiation budget is negative (cooling) for the first 13 months after the eruption: DRE(TOA) = −0.82 ± 0.35 W m⁻², and the aerosols (WV) are responsible for ∼ 95 (5) % of this cooling. The breakdown in short- and mid-term tendencies shows a decrease in the net DRE(TOA) from −1.01 ± 0.34 W m⁻² during M2–M4 to −0.76 ± 0.35 W m⁻² during M5–M14. Locally, Sellitto et al. (2022) found the instantaneous aerosol and water vapor DRE(TOA) was negative in the first 2 weeks (ranging from −12.5 to −20.1 W m⁻²) and positive (+0.2 W m⁻²) for what they call the “aged plume” (8 February 2022). At the global scale, Zhu et al. (2022), mentioned earlier, find a “TOA instantaneous aerosol effect” (equivalent to DRE(TOA) here) of −0.19 W m⁻² for the first 2 months of 2022. Gupta et al. (2025) estimated a net instantaneous clear-sky radiative energy loss of −0.48 ± 0.06 W m⁻² at TOA in the Southern Hemisphere, resulting from its effects on stratospheric water vapor, aerosols, and ozone, and showed a certain zonal homogeneity over a 1-year period. Our results, representative of Réunion Island and likely of the zonal hemisphere at 21° S, are a little less than 2 times larger than those of Gupta et al. (2025), manifesting that the aerosol distribution is not latitudinally uniform. At the surface, a marked cooling produced by the Hunga volcanic aerosols is observed (DRE(BOA) = −1.04 ± 0.36 W m⁻² for the entire period), with a decreasing tendency in absolute value with time (DRE(BOA) = −1.26 ± 0.50 W m⁻² for M2–M4 and −0.98 ± 0.32 W m⁻² for M5–M14). At the global scale, Zhu et al. (2022) find an “adjusted radiative effect” (not exactly equivalent to DRE(BOA) here) of −0.21 W m⁻² for the first 2 months of 2022. Also, Zuo et al. (2022) modeled the global surface temperature in the first year after the Hunga volcano eruption and found a negative anomaly of −0.004 K but recognized that it is “within the amplitude of internal variability at the interannual timescale and thus not strong enough to have significant impacts on the global climate”. Our analysis shows that the eruption of the Hunga volcano might have had a clear cooling impact at the surface locally and likely zonally at 21° S. An interesting result at this point is how the aerosol and water vapor H $/$ C rates distribute vertically in the atmosphere. It is clear in Fig. 6d that the negative longwave H $/$ C rate caused by water vapor and the positive one caused by the aerosols coexist at different altitude levels. During M2–M4, the small height difference between the sulfate and the moist layers and the higher rate of cooling of the latter result in a negative aerosol and water vapor H $/$ C rate in most of the altitude range considered (solid red line in Fig. 7). During M5–M14, H $/$ C rate profiles show a clear vertical difference locally in the stratosphere between the aerosol warming impact (18 to 26 km) and the water vapor cooling (22 to 30 km). The resulting aerosol and water vapor H $/$ C rate profile follows an S-shaped curve with peaks slightly larger for the moist layer (−0.09 K d⁻¹ at 26 km) than for the sulfate layer (+0.06 K d⁻¹ at 22 km).

Put in the context of other recent volcanic studies in the literature, the Hunga radiative effect at TOA is quite unique as far as the combined aerosol and WV impact is concerned, and the literature with which to compare our results is rather limited. A general agreement is observed for the shortwave aerosol component, which is systematically found negative in all studies (volcano/year of the eruption: Raikoke/2019, Ulawun/2019 (Kloss et al., 2021), Ambae/2018 (Kloss et al., 2020), Kasatochi/2008, Sarychev/2009, and Nabro/2011 (Andersson et al., 2015)) and through multi-year (2000–2013) global volcanic aerosol forcing estimations (Ridley et al., 2014). For the longwave component, persuasive evidence of the volcanic longwave effect was missing for a long time in the past. However, the longwave effect of stratospheric sulfate aerosols is now well understood and calculated properly (Schmidt et al., 2018). Our work emphasizes the importance of including the aerosol and water vapor interaction with the longwave radiation in stratospheric volcanic studies, especially in the case of submarine eruptions. It also puts forward the complexity of correctly estimating the TOA radiative forcing and thus the ambiguity generated when different datasets, models, spatiotemporal averaging, etc. are used. The use of unambiguous, direct measurements of radiative fluxes from space, such as in Minnis et al. (1993), are highly encouraged.

5 Conclusions

A whole 13 months after the eruption of the Hunga volcano, aerosols and water vapor are still present in the stratosphere of the southern tropical Indian Ocean region. During the first 3 months after the eruption, the stratospheric aerosol optical depth increases and reaches a peak at 0.035 (∼ 13 times the background sAOD) in April 2022, the highest in the last decade and the third highest in the last 40 years (after Pinatubo and El Chichón). From April to November 2022, the sAOD decreases and then stabilizes at a value of 0.012 (∼ 5 times the background sAOD). Unusually high water vapor concentrations are also observed in the stratosphere. On a monthly basis, the water vapor stratospheric peak reaches a maximum in February 2022, which is approximately a factor 5 above the climatological reference. By February 2023, this ratio decreases down to almost 2.

In all moments, the water vapor plume is located at a higher altitude than the aerosol plume. The heights of the peak of the aerosol and water vapor layers have opposite tendencies as of April 2022: the aerosol plume is slowly descending, mostly by gravitational settling, whereas the moist layer ascends slowly until October 2022 and at a higher rate afterwards. The upward transport of the moist layer is due to the Brewer–Dobson circulation. Both aerosol and WV plumes are still present on 15 April 2023. The aerosol plume is located at 18.5–23.5 km height, and the moist layer is above 30 km. As far as aerosols are concerned, the plume peak height decreases from April 2022 at an average steady rate of ∼ −8 m d⁻¹. Assuming this rate is constant in time and omitting other removal processes, the remaining lifetime of the volcanic plume in the stratosphere is estimated to be between 2 and 2.5 years after 15 April 2023.

The radiative impact of both aerosol and water vapor layers is estimated at our site. Averages are made over three different periods of time in order to explain the temporal evolution: the first 13 months and the short-term (M2–M4) and mid-term (M5–M14) periods. During the first 13 months after the Hunga volcano eruption, the overall aerosol and water vapor impact on the Earth's radiation budget is negative (cooling, −0.82 ± 0.35 W m⁻²) and dominated by the aerosol impact (∼ 95 %; the remaining ∼ 5 % is due to WV). At the Earth's surface, aerosols are the main driver and produce a negative (cooling, −1.04 ± 0.36 W m⁻²) radiative impact. Between the short- and mid-term periods, the aerosol and water vapor radiative effect at both the surface and TOA reduces by 22 % and 25 %, respectively. Heating $/$ cooling rate profiles show a clear vertical difference in the stratosphere during the mid-term period between the aerosol warming impact (18 to 26 km) and the water vapor cooling (22 to 30 km). The mean profile during the mid-term period follows an S-shaped curve with peaks slightly larger for the moist layer (−0.09 K d⁻¹) than for the sulfate layer (+0.06 K d⁻¹). This study shows that the eruption of the Hunga volcano has so far had a clear radiative impact on the Earth's radiation budget locally and likely zonally at 21° S.