Multi-model analysis of the impact of water vapor on the radiative forcing of volcanic aerosols after the 2022 Hunga Eruption

On 15 January 2022, the Hunga volcano eruption released unprecedented amounts of water vapor into the atmosphere alongside a modest amount of SO₂. In this work we analyse results from multiple Earth system models as part of the Hunga Tonga-Hunga Ha'apai Volcano Impact Model Observation Comparison Project. Our results show a good model agreement over the climatic outcomes of the eruption, overall indicating a significant negative radiative forcing from the Hunga eruption. The multi-model mean of global instantaneous radiative forcing averaged over 2022–2023 is estimated at −0.19 ± 0.06 W m⁻² at the top-of-atmosphere (TOA), and −0.16 ± 0.06 W m⁻² at the surface. Simulations with free-running meteorology and climatological sea surface temperatures and sea ice yield a global mean TOA forcing of −0.14 ± 0.10 W m⁻² across two models for the first 2 years, decreasing to −0.09 ± 0.10 W m⁻² on average between 2022 and 2027. However, these global values may be underestimated by about 50 %, considering that recent SO₂ injection retrievals suggest nearly twice the amount than the 0.5 Tg-SO₂ used in the protocol. We also find that the contribution from added stratospheric water vapor is minimal and that the injected SO₂ and the resulting formation of stratospheric sulfate dominate the radiative forcing. However, water vapor played a key role in the initial aerosol growth, leading to a stronger negative radiative forcing during the first 6 months after the eruption compared to simulations without water vapor co-injection.

1 Introduction

Explosive volcanic eruptions that release large amounts of SO₂ into the stratosphere have long been identified as a significant contributor to the global energy budget (Schmidt et al., 2018). The release of SO₂ into the stratosphere results in the formation of sub-micron, supercooled sulfate aerosols which efficiently reflect incoming solar radiation, resulting in a negative forcing at the top of the atmosphere and thus a corresponding surface cooling (Kremser et al., 2016). While often studies of volcanic aerosols' impact on climate focus on very large eruptions like Mt. Pinatubo, that erupted in 1991 (Quaglia et al., 2023) releasing anywhere between 10 and 20 Tg-SO₂ in a few days (Baran and Foot, 1994; Fisher et al., 2019), recent decades have seen no eruptions with injections exceeding a few Tg of SO₂, but rather multiple moderate stratospheric eruptions (Carn et al., 2017; Brodowsky et al., 2021). Instead, moderate eruptions, such as Raikoke (2019), Ulawun (2019), Ambae, La Soufrière (2021), and Hunga (2022), have received increasing attention in recent research, offering valuable insights into aerosol microphysics, plume dynamics, and radiative forcing at the top of the atmosphere (Kloss et al., 2021; Wrana et al., 2023). Despite their smaller scale compared to Pinatubo, moderate eruptions can still substantially double the quiescent stratospheric sulfate burden (Schmidt et al., 2018; Andersson et al., 2015; Brodowsky et al., 2024), underscoring the importance of including these events in both climate modeling and aerosol research.

The co-emission of other byproducts of volcanism, together with sulfate, is not infrequent: some ash and water vapor usually reaches the stratosphere (Zhu et al., 2020), and in some occasions chlorine and bromine in significant quantities are also co-injected, but whether they reach the stratosphere are debatable (Staunton-Sykes et al., 2021). Hunga erupted on 15 January 2022 in the South Pacific, releasing an amount of materials into the atmosphere unprecedented in the satellite era (Carr et al., 2022). Current estimates of injected material are 146 ± 5 Tg of water vapor based on Aura Microwave Limb Sounder (MLS) retrievals (Millán et al., 2022), with a broader range for injected SO₂, ranging from 0.41 ± 0.02 Tg of SO₂ for MLS to values greater than 1.0 Tg of SO₂ from the Infrared Atmospheric Sounding Interferometer (IASI) (Sellitto et al., 2024). While the amount of injected sulfur was within the range of many moderate past explosive volcanic eruptions, the amount of water vapor co-injected was unprecedented in the historical record, representing roughly ∼ 10 % of all background stratospheric water vapor. Although comparable increases in moisture can be achieved through the indirect pathway, such as during Pinatubo-like eruptions in the tropics (Kroll and Schmidt, 2024), what made the Hunga eruption truly unprecedented was the height at which the overall volcanic plume was injected, reaching up to 55 km (Carr et al., 2022). Early estimates assumed the water vapor forcing would result in a small net warming (Jenkins et al., 2023). However, other estimates found that the sulfate aerosols produced by the co-injections of SO₂, which normally oxidize to form sulfuric-acid aerosols, would result in a significant negative forcing after the volcanic cloud spread uniformly (Sellitto et al., 2022; Zhu et al., 2022; Schoeberl et al., 2024), larger than the water vapor forcing. Later modeling analyses by Stenchikov et al. (2025), using WRF-chem, also found a negative forcing by considering similar factors. Notably, the conditions of significant stratospheric hydration have been suggested as an important factor contributing to the substantial radiative impact of the Hunga cloud, even in the case of a modest SO₂ injection. This is because enhanced hydration promoted faster aerosol growth and the formation of aerosol that reached optically efficient scattering sizes (see Fig. 7 in Murphy et al., 2021) in a couple of weeks, rather than the 4 months observed after the Pinatubo eruption (Zhu et al., 2022; Li et al., 2024; Sellitto et al., 2025; Asher et al., 2025).

(former name in  Zhu et al., 2025)

Table 1 Summary of experiments.

Table 2 Injection parameters. Adapted from Table 7 in Zhu et al. (2025)

As part of the Hunga Tonga-Hunga Ha'apai (HTHH) Impacts activity that was established in the World Climate Research Programme (WCRP) Atmosphere Processes And their Role in Climate (APARC), a multi-model-observation intercomparison named the Hunga Tonga-Hunga Ha'apai Volcano Impact Model Observation Comparison (HTHH-MOC) Project was proposed in Zhu et al. (2025) with the aim to better ascertain the radiative and climatic impacts of Hunga in a multi-model context and to try and separate the volcanic impacts from other natural or anthropic perturbations in the same years (Forster et al., 2023), especially in light of the “unprecedented” surface temperature warming in 2023 (Cattiaux et al., 2024; Quaglia and Visioni, 2024). Other parallel works within this same project discuss in depth changes to stratospheric water vapor, aerosols, temperature, and ozone (Zhuo et al., 2025) and potential impacts on climate (Bednarz et al., 2026). In this paper, we examine the radiative forcing resulting from these changes using a multi-model approach, by using several proposed HTHH-MOC experiments. Our goal is to highlight the sources of agreement and differences across models in estimating the magnitude of the forcing. Three different atmospheric setups are considered for the Hunga eruption simulation: (1) temperature and meteorology nudged to observations; (2) a free running meteorology setup with fixed sea surface temperature and sea ice extent; (3) and free-running meteorology setup with interactive ocean. The framework allows to distinguish the contributions from the direct interactions between the forcing agent and radiation, the rapid adjustments to these forcings, and slower climate feedback. Additionally, sensitivity studies performed within the nudged setup compare the effects of SO₂-only and water vapor-only injections with the combined impact of both. The long-term evolution of Hunga water vapor and aerosols in the free-running models provides broader projections of the eruption climate impact. This comprehensive approach differentiates our study from previous work by integrating multiple experimental setups to better capture the complex interactions and feedbacks associated with the eruption.

2 Methods

2.1 Models and simulations

We use a suite of experiments described in depth in Zhu et al. (2025), which includes also a detailed descriptions of the models under analyses. The first set of experiments is a 2-year experiment over 2022–2023 with nudged temperature and meteorology as well as observed sea surface temperatures (SSTs) and sea ice (Exp2a in Zhu et al., 2025, here denoted “Nudged”; note this is distinct from Exp2b, not included in this study, which does not require prognostic aerosols). The second set of experiments is a 10-year long experiment over 2022–2031 with free running meteorology, using either imposed climatological SSTs and sea ice (Exp1_fixedSST in Zhu et al., 2025, here “Fixed-SST”) or with the atmosphere model coupled with the ocean (Exp1_coupled, here “Coupled”). Each experiment includes at least one simulation that involves the combined injection of 0.5 Tg of SO₂ and 150 Tg of H₂O (“SO2andH2O”), alongside a control simulation without any injections (“NoVolc”). The Nudged experiment also includes single-forcing simulations that inject only SO₂ (“SO2only”) or only H₂O (“H2Oonly”). The location, altitude, and amount injected are different between models in order to better match the observed plume in the first couple of days, but they remain consistent across experiments between models. Models account for the interactive coupling between aerosol, water vapor, radiation, and dynamics, allowing the fast descent due to water vapor longwave cooling (Sellitto et al., 2022) to be simulated through this coupling. A summary of the experiments is provided in Table 1, and detailed model-specific injection settings are listed in Table 2.

Nudged was conducted by 5 models: the Whole Atmosphere Community Climate Model version 6 (WACCM6;  Gettelman et al., 2019; Davis et al., 2023) coupled with the four-mode modal aerosol module (MAM4,  Liu et al., 2012, 2016; Mills et al., 2016) and the Community Aerosol and Radiation Model for Atmospheres (CARMA;  Tilmes et al., 2023), WACCM6-MAM and WACCM6-CARMA, respectively; the Model for Interdisciplinary Research On Climate – CHemical Atmospheric general circulation model for Study of atmospheric Environment and Radiative forcing version 6 (MIROC-CHASER;  Sekiya et al., 2016); the atmospheric component of CESM1, the Community Atmosphere Model version 5 (CAM5;  Lamarque et al., 2012) using the sectional aerosol microphysics model CARMA (CAM5-CARMA;  Yu et al., 2015), and the UK Earth System Model version 1.1 (UKESM;  Mulcahy et al., 2023). Fixed-SST was carried out with two models, WACCM6-MAM and MIROC-CHASER, while only WACCM6-MAM participated in Coupled. WACCM6-MAM ran 30-member ensembles for both Fixed-SST and Coupled, whereas MIROC-CHASER ran 10-member ensembles for Fixed-SST. Further information about the participating models can be found in the references above, as well as in Zhu et al. (2025).

2.2 Observational datasets

We use two observational datasets in this study: the Global Space-based Stratospheric Aerosol Climatology version 2.22 (GloSSAC,  NASA/LARC/SD/ASDC, 2023) for zonal monthly-mean stratospheric aerosol optical depth (AOD), and the Stratospheric Water and Ozone Satellite Homogenized dataset version 2.6 (SWOOSH,  Davis et al., 2016) for water vapor.

GloSSAC provides a long-term and global record of stratospheric aerosol properties, including the stratospheric AOD at 525 nm used here. It is primarily based on Stratospheric Aerosol Gas Experiment (SAGE) measurements up to mid-2005, with Optical Spectrograph and Infrared Imaging System and Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations used thereafter, and SAGE III/ISS extending the climatology to the present. Additional data from other satellites, as well as ground-based, airborne, and balloon-borne instruments, are included to fill observational gaps (Thomason et al., 2018; Kovilakam et al., 2020). Although SAGE III/ISS data have limited spatial and temporal coverage, particularly during the first weeks following the Hunga eruption, they are considered robust because they rely on solar occultation and do not require assumptions about aerosol type or particle size distribution. In contrast, limb-scatter datasets such as Ozone Mapping and Profiler Suite Limb Profiler (OMPS-LP), used for comparison in Zhuo et al. (2025), depend on additional assumptions about aerosol properties (Kovilakam et al., 2025). For these reasons, GloSSAC is adopted as the primary observational reference, providing a consistent and reliable basis for evaluating modeled stratospheric aerosol evolution over longer timescales.

SWOOSH is a long-term and global record of stratospheric ozone and water vapor measurements from multiple satellite instruments (SAGE II, Meteor-3M SAGE III, Halogen Occultation Experiment, Upper Atmosphere Research Satellite Microwave Limb Sounder, and Earth Observing System Aura Microwave Limb Sounder) spanning 1984 to the present.

2.3 Radiative forcing estimations

There are important differences among the radiative forcing estimates of volcanic eruptions derived from nudged simulations, which provide the instantaneous radiative forcing (IRF); free-running atmosphere-only simulations, which yield the effective radiative forcing (ERF), following the definition of Forster et al. (2016); and fully coupled simulations with an interactive ocean, which are generally referred to as radiative forcing (RF).

IRF represents the combined effect of direct interactions between the forcing agent and radiation, as well as interactions between the forcing agent and clouds (Smith et al., 2018). ERF, defined by Myhre et al. (2013) as “change in the net TOA downward radiative flux after allowing for atmospheric temperatures, water vapour and clouds to adjust, but with surface temperature or a portion of surface conditions unchanged”, is the sum of IRF and rapid adjustments. These rapid adjustments which occur over weeks to months, before global-mean surface temperatures can respond, are due to changes in tropospheric and stratospheric temperature, water vapor, surface albedo, and clouds, and are distinct from slower feedbacks that are driven by surface temperature changes (Smith et al., 2018; Sherwood et al., 2015). In coupled model simulations, RF includes IRF, rapid adjustments to that forcing, and slower climate feedbacks resulting from the coupled ocean–atmosphere response (Chung and Soden, 2015).

The radiative forcings are calculated under both Clear-Sky (CS, without clouds) and All-Sky (AS, including the effects of clouds) conditions. Unless otherwise specified, all values are assumed to be calculated under Clear-Sky conditions.

IRF is calculated in climate models using a double radiation call that excludes aerosols from online radiative calculations (“Clean-Sky”) following the method proposed in Stenchikov et al. (1998) in order to separate the contribution of aerosols from that of other components. Since nudging reduces variability in meteorological fields, it typically limits any stratospheric temperature adjustments, as radiative forcing is calculated as the difference between the perturbed and unperturbed case (Eq. 1). As a result, the IRF from aerosols and water vapor, either combined or individually, can be approximated using the corresponding Nudged experiments (Eqs. 2 and 3, respectively).

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathrm{IRF}=F_{\text{SO2andH2O,\hspace{0.17em}\hspace{0.17em}Nudged}}-F_{\text{NoVolc,\hspace{0.17em}\hspace{0.17em}Nudged}}\text{(2)}& {\displaystyle }& {\displaystyle }{\mathrm{IRF}}_{\mathrm{aerosol}}=F_{\text{SO2only,\hspace{0.17em}\hspace{0.17em}Nudged}}-F_{\text{NoVolc,\hspace{0.17em}\hspace{0.17em}Nudged}}\text{(3)}& {\displaystyle }& {\displaystyle }{\mathrm{IRF}}_{\mathrm{gas}}=F_{\text{H2Oonly,\hspace{0.17em}\hspace{0.17em}Nudged}}-F_{\text{NoVolc,\hspace{0.17em}\hspace{0.17em}Nudged}}\end{array}$$

In the case of WACCM6-MAM simulations, both methodologies have been applied. Notably, in contrast to previous studies, all forcing estimates presented here treat stratospheric water vapor as a forcing rather than a feedback, due to its direct injection. The second estimate of IRF in WACCM6-MAM for aerosols (sulfate, black carbon, primary organic matter, secondary organic aerosols, sea salt, and dust) and gases (water vapor, dioxygen, carbon dioxide, ozone, nitrous, methane, chlorofluorocarbons) is derived through a double radiation call (Eqs. 4 and 5, respectively). Within this approach, F^clean denotes the Clean-Sky calculation, which excludes aerosols.

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{rl}{\mathrm{IRF}}_{\mathrm{aerosol}}=& (F_{\text{SO2andH2O,\hspace{0.17em}\hspace{0.17em}Nudged}}-F_{\text{NoVolc,\hspace{0.17em}\hspace{0.17em}Nudged}})\\ & -\left(F_{\text{SO2andH2O,\hspace{0.17em}\hspace{0.17em}Nudged}}^{\text{clean}}-F_{\text{NoVolc,\hspace{0.17em}\hspace{0.17em}Nudged}}^{\text{clean}}\right)\end{array}\text{(5)}& {\displaystyle }{\mathrm{IRF}}_{\mathrm{gas}}=F_{\text{SO2andH2O,\hspace{0.17em}\hspace{0.17em}Nudged}}^{\text{clean}}-F_{\text{NoVolc,\hspace{0.17em}\hspace{0.17em}Nudged}}^{\text{clean}}\end{array}$$

ERF and RF are calculated as the difference between the perturbed and unperturbed case (Eqs. 6 and 7, respectively).

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }\mathrm{ERF}=F_{\text{SO2andH2O,\hspace{0.17em}\hspace{0.17em}Fixed-SST}}-F_{\text{NoVolc,\hspace{0.17em}\hspace{0.17em}Fixed-SST}}\text{(7)}& {\displaystyle }& {\displaystyle }\mathrm{RF}=F_{\text{SO2andH2O,\hspace{0.17em}\hspace{0.17em}Coupled}}-F_{\text{NoVolc,\hspace{0.17em}\hspace{0.17em}Coupled}}\end{array}$$

We calculate the different radiative forcings at three key atmospheric levels – top of the atmosphere (TOA), tropopause (TROP), and surface (SURF), as the sign and magnitude of RF can differ by altitude and carry distinct physical implications. The radiative forcing at the three levels is presented for each model according to data availability. Radiative forcing at TOA reflects the overall perturbation to the Earth's energy budget and is commonly used to estimate the potential influence on global average temperature. At the tropopause, RF captures the net energy change affecting the coupled troposphere–surface system, and is considered less affected by upper stratospheric processes and better represents tropospheric heating. In contrast, RF at the surface does not directly correspond to surface temperature responses but is more relevant for understanding impacts on the hydrological cycle, particularly changes in precipitation patterns (Ramaswamy et al., 2018).

Figure 1 Time series of zonal mean radiative forcing at the top-of-atmosphere (TOA) under clear-sky conditions from five models: WACCM6-MAM (a–c), WACCM6-CARMA (d–f), MIROC-CHASER (g–i), CAM5-CARMA (j–l), and UKESM (m). Each column corresponds to a different perturbation scenario from the nudged experiment: the first column shows SO2andH2O, the second column SO2only, and the third column H2Oonly. The data gap is due to missing model output for those months.

3 Results and discussion

We present first the results based on the nudged simulations (Sect. 3.1), as this experiment was done by more models (five) and the use of observed meteorological conditions also allows us to better distinguish the overall radiative impact by removing natural variability. Since the same nudging (of temperature and horizontal wind) was applied in both the control and the volcanic injection experiments, taking the difference between the two simulations, as done in Fig. 1, isolates the direct radiative forcing from the volcanic material, as any temperature response cancels out. Following that, we provide analyses of the two models which provided the 10-year long free-running simulations with prescribed climatological SSTs and sea ice (Sect. 3.2): these analyses allow us to understand the long-term behavior of the forcing as well as include the combined chemical and dynamical impacts and temperature adjustments. Finally, we complete those with an analyses of fully-coupled simulations in WACCM6-MAM (Sect. 3.3). In this case, we compare the results from the three different types of simulations within the same model, coupled with analyses of the double-radiation call described in the methods above, to discuss and quantify the different contribution to the forcings analyses elsewhere.

3.1 Multi-model comparison of instantaneous radiative forcing in the nudged simulations

The IRF includes all optically active components, such as stratospheric aerosols, water vapor, ozone, and polar stratospheric clouds (PSCs). For model intercomparison, we use clear-sky forcing to minimize uncertainties associated with aerosol–cloud interactions that arise from differences in cloud parameterizations among models.

All models show qualitative agreement in the spatial distribution of the IRF under clear-sky conditions at TOA (Fig. 1), at TROP (Fig. A1), and at SURF (Fig. A2). At each level, the response is primarily located in the Southern Hemisphere (SH); therefore, most of the following analyses will focus on the SH only.

A consistent pattern of negative forcing appears across models in both the SO2andH2O and SO2only experiments (first and second columns in Figs. 1, A1, and A2). This forcing peaks in the tropics during the first few months following the eruption and then moves to mid- to high-southern latitudes by 2023, with substantial negative forcing persisting through the end of 2023 in the SH. Likewise, all models simulate a negligible IRF from the H2Oonly experiment (third column of the same figures). Therefore, the models consistently attribute nearly all of the radiative forcing to the aerosol perturbations rather than to the injected water vapor.

Figure 2 Time series of net radiative forcing under clear-sky conditions at three atmospheric levels: TOA (a, d, g), TROP (b, e, h), and SURF (c, f, i). Solid lines indicate the multi-model means for the nudged experiments as indicated in the legend, while the shaded areas represent the range (minimum to maximum) across models. Each row corresponds to a different latitudinal averaging region: the first row covers 90° N–S, the second 60° S to the equator, and the third row 60–90° S.

Global means (Fig. 2a–c) and hemispheric means (60° S-EQ in Fig. 2d–f; 60–90° S in Fig. 2g-i) show the multi-model average and inter-model spread, and highlight the magnitude of the IRF in each injection experiment. These means further clarify that the IRF from H2Oonly is negligible: it is slightly negative at TOA (first column) and slightly positive at TROP and SURF (second and third columns). This behavior results from the vertical distribution of water vapor in the stratosphere and the lack of stratospheric temperature adjustment, which together increase outgoing longwave radiation while enhancing downward longwave re-emission. Especially at the surface, the IRF remains positive but very small, as the radiative forcing is dominated by sensible and latent heat fluxes, which are controlled by the prescribed surface temperature in these simulations.

For the aerosol forcing, which constitutes essentially the entirety of the Hunga radiative forcing, the models generally agree on its temporal evolution but exhibit substantial differences in its magnitude. The multi-model mean TOA IRF from SO2andH2O is −0.35 ± 0.12 W m⁻² over 60° S-EQ (2022–2023 average; Fig. 2d–f), where the inter-model spread is roughly one third of the mean magnitude. The global mean forcing is approximately half of the 60° S-EQ values, reflecting the strong SH dominance of the perturbation. Over the high southern latitudes (60–90° S; Fig. 2g–i), the mean TOA IRF is −0.26 ± 0.23 W m⁻², indicating large uncertainty, with the standard deviation nearly matching the mean (September 2022–December 2023 average). Comparing the injection experiments, the IRF from SO2andH2O is generally comparable to SO2only at TOA, but less negative at TROP and SURF: over the 60° S-EQ, TROP IRF is −0.32 ± 0.13 for SO2andH2O versus −0.43 ± 0.12 W m⁻² for SO2only; at the surface, these values are −0.29 ± 0.12 versus −0.36 ± 0.11 W m⁻².

The shortwave (SW) component of the IRF follows the evolution of stratospheric AOD (Fig. A3) but is additionally modulated by the strong seasonality of solar insolation, producing relative minima during the Austral winter. This is particularly evident at high southern latitudes (Fig. 2g–i), where the small negative SW IRF is partially offset by a positive longwave (LW) contribution. Generally, the LW signal coincides with the peak in stratospheric AOD, which, over the high southern latitudes, occurs during the Austral winter of 2023.

The SW IRF is similar between TOA and TROP but smaller in magnitude at SURF in both SO2andH2O and SO2only experiments. However, at all levels, the SW IRF is more negative in SO2andH2O during the first few months after the eruption and then becomes less negative compared to SO2only, consistent with the evolution of stratospheric AOD in the two experiments (Fig. A3j and k). The differences in the net IRF between the two experiments decrease at TOA but increase at TROP due to the different LW responses to water vapor, which is slightly negative at TOA and positive at TROP.

Figure 3 Time series of net radiative forcing at TOA, averaged between 60° S and the equator, in WACCM6-MAM. Panel (a) shows the total IRF (aerosols + gases), (b) shows the aerosol-only contribution, and panel (c) shows the gas-only contribution, all derived using the double radiation call method (see Sect. 2.2). Different colors indicate various nudged experiments, with solid lines representing clear-sky conditions and dotted lines indicating all-sky conditions. Panels (d) and (e) show the 6-month and 2-year averages, respectively, of the radiative forcing from the three nudged experiments, under clear-sky (CS) and all-sky (AS) conditions. The contributions from aerosols and gases alone are highlighted relative to the total forcing.

Figure 3 provides further insight into the TOA IRF contributions from aerosols (IRF_aerosol) and water vapor (IRF_gas), based on analyses from WACCM6-MAM using both methods described in Sect. 2.2. Among the models, WACCM6-MAM shows the strongest aerosol response (first and second columns of Fig. A4). Multiple interesting features emerge from an analysis of Fig. 3, which also includes a comparison of clear-sky and all-sky IRF.

• The IRF from the SO2andH2O and SO2only simulations (black and orange lines in Fig. 3a) show similar behavior from late 2022 onward. However, substantial differences are evident during the first 6 months (on the order of 20 % in WACCM6-MAM, Fig. 3c), with the exact duration being model dependent (Fig. A4). This is in agreement with previous studies indicating that the co-injection of stratospheric water vapor promotes faster particle growth to optically efficient scattering sizes, reaching an effective radius of about 0.4 µm over less than a month (Fig. A5a–c for SO2andH2O versus Fig. A5d–f for SO2only). Further details on the microphysical aspects, using the same protocol and model simulations, are provided in Khaykin et al. (2025) (comparison of the particle size with observation in Fig. 3.18). As a result, SO2andH2O produces a larger initial stratospheric AOD (Fig. A3j and k) and corresponding negative radiative forcing. Later on, the forcing becomes similar due to a compensating effect: SO2only shows a larger stratospheric AOD, while SO2andH2O shows a stronger gas contribution (negative at TOA). The SO2andH2O forcing becomes slightly weaker only in the last few months of 2023, when the gas IRF in the two simulations becomes comparable (Fig. 3c).

• A comparison between the aerosol IRF (IRF_aerosol; labeled as “Aerosol” in the figure) and the total IRF indicates that IRF_aerosol accounts for approximately 80 % of the overall IRF response, with the remaining contribution attributed to changes in other stratospheric gases (see Fig. 3b vs Fig. 3a; time-averaged values are summarized in Fig. 3d and e).

• H2Oonly simulations do not show any significant forcing in clear-sky, and a more negative, but also highly variable, forcing in all-sky (blue lines in Fig. 3a–c). This further confirms that the water vapor forcing alone does not result in either a cooling or a warming. The small impact of stratospheric water vapor is altitude-dependent, with the largest effect occurring near the tropopause rather than at higher altitudes (Solomon et al., 2010). Indeed, stratospheric water vapor concentrations are primarily confined to the middle stratosphere (40–10 hPa) by the end of 2022, after which they rise and extend into the lower mesosphere (Millán et al., 2024; Zhuo et al., 2025). Interestingly, a comparison of the IRF_gas in the SO2andH2O and H2Oonly experiments (black and blue lines in Fig. 3c) suggests that the all-sky contribution is not coming from water vapor itself, but from changes to other radiatively-active gases. This will be further discussed in Sect. 3.3.

• The estimation of IRF_aerosol and IRF_gas differs depending on whether they are derived from the SO2only and H2Oonly experiments, respectively, or from the double radiation call method applied to SO2andH2O experiment (Fig. 3d and e). Under both clear-sky (CS) and all-sky (AS) conditions, during the first 6 months after the eruption, IRF_aerosol and IRF_gas obtained from the SO2only and H2Oonly experiments, respectively, are smaller than those calculated using the double radiation call in the SO2andH2O experiment (green and red boxes, respectively). This discrepancy occurs because, in the SO2andH2O experiment, water vapor influences aerosol growth, thereby affecting IRF_aerosol, while any simultaneous changes in other gases (e.g. ozone) will also contribute to IRF_gas.

Figure 4 (a–c) Time series of moving average of effective radiative forcing (ERF) under clear-sky conditions, averaged between 60° S and the equator, from the free running simulatons with fixed climatological SSTs and sea ice (Fixed-SST). ERF is shown at TOA, TROP, and SURF for WACCM6-MAM and MIROC-CHASER. Solid lines represent the ensemble mean, and shaded areas represent ± 1 standard deviation across ensemble members. Each shaded area uses the same color as the line representing that model, as shown in the legend. Zonal mean of changes in stratospheric aerosol optical depth (sAOD) for WACCM6-MAM, MIROC-CHASER and GloSSAC (d, g, j), tropical (30° S–N) H₂O mixing ratio for WACCM6-MAM, MIROC-CHASER and SWOOSH (e, h, k) and O₃ mixing ratio for WACCM6-MAM and MIROC-CHASER (f, i). Changes are calculated as the difference between perturbed and unperturbed simulations for the models, and relative to the 2005–2021 climatology period for the observations. Gray areas indicate regions where the differences are not statistically significant at the 5 % level based on Student's t-test.

Table 3 Effective radiative forcing (in W m⁻²) at the Top of Atmosphere under clear-sky conditions. Forcing is averaged between 60° S and the equator.

3.2 Comparison of the effective radiative forcing in the free running simulations

Figure 4 shows the longer-term evolutions of the moving average forcing, calculated as:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{1}}{N_t}}\sum _{i=\mathrm{1}}^{N_t}{\text{RF}}_i\end{array}$$

where N_t represents the number of months elapsed, and RF_i is the radiative forcing at month i. We use the moving average forcing instead of the instantaneous forcing to highlight the time-averaged impact of the eruption over time, which provides a better measure of its overall climatic effects. This is shown for the free-running simulations with fixed climatological SSTs and sea ice (Fixed-SST). For this experiment, only two models, WACCM6-MAM and MIROC-CHASER, performed the simulations. The TOA ERF estimates show a very good agreement between the two models, both in terms of the overall magnitude and the rate of dissipation of the forcing, which is partly true also for the ERF estimates at the surface (see Tables 3 and A3). Larger differences between the models are present for the ERF estimates at the tropopause, whereby WACCM6-MAM shows a more negative average forcing (by 0.1 W m⁻²) in the first year (Table A2) and a different overall trend throughout the simulations.

Although the ERF response is broadly consistent between the two models, notable differences emerge in their simulated stratospheric AOD (Fig. 4d and g). In particular, we observed significant discrepancies in the timing of aerosol formation and the onset of its decline. In MIROC-CHASER, stratospheric AOD starts to decline 2 months after the eruption, right after reaching its peak. In contrast, WACCM6-MAM shows a peak in stratospheric AOD 3 months post-eruption, which remains at peak levels for 6 months before beginning to decrease. However, the models generally agree on the volcanic cloud's southward transport, with the volcanic cloud moving towards 60° S after 4 to 5 months and reaching the pole by December 2023. The persistence of the aerosol cloud at low-to-mid latitudes in the SH is different between the two models, with WACCM6-MAM showing some persistence all the way to the end of 2023, whereas MIROC-CHASER show a full dissipation by late 2023.

Contrary to the stratospheric AOD, the water vapor upward diffusion (Fig. 4e and h) shows more persistence in MIROC-CHASER than WACCM6-MAM, with the latter showing no residual water vapor anomaly below 1 hPa by late 2026, whereas MIROC-CHASER shows at least 0.5 ppm more in the upper stratosphere remaining all the way to 2028. MIROC-CHASER also shows a stronger positive tropical ozone anomaly between 10 and 20 hPa compared to WACCM6-MAM, and a negative ozone anomaly higher up between 1 and 2 hPa, with the negative anomaly more persistent over time in both models (Fig. 4f and i), explained by a potential enhancement of the HO_x-driven loss cycle due to the water vapor anomaly (Randel et al., 2024; Fleming et al., 2024).

As with the instantaneous radiative forcing, the gas response varies depending on the atmospheric level at which it is calculated and the vertical distribution of the gases. Additionally, temperature adjustments are included here. We speculate that the stronger stratospheric AOD anomaly in WACCM6-MAM, along with the stronger water vapor and ozone anomalies in MIROC-CHASER, may either offset the radiative forcing at the top of the atmosphere (TOA) or amplify it at the tropopause (TROP). This could explain the similar forcing observed in Fig. 4a and the larger differences seen in Fig. 4b. This aspect is only further explored in Sect. 3.3 for WACCM6-MAM, where the availability of separate aerosol and gas radiative contributions allows for a clearer disentangling of each factor.

In Fig. 4 we also provide a comparison with available observations, with more in depth comparisons presented in Zhuo et al. (2025). In general, both models show good qualitative agreement with the observations. While GloSSAC does not capture the pronounced peak immediately following the eruption, which is evident in other observational datasets discussed in Zhuo et al. (2025), this mainly reflects the limited spatiotemporal coverage of the SAGE III/ISS observations on which GloSSAC is based during the earliest post-eruption phase. Despite this limitation, transport toward the Southern Hemisphere occurs on similar timescales in both the observations and the simulations, with WACCM6-MAM showing a better match in terms of the residual aerosol cloud around 60° S. In terms of H₂O, models generally reproduce the upward transport pattern; however, SWOOSH consistently reports higher values. As also noted in Zhuo et al. (2025), the comparison between models and observations is primarily aimed at assessing transport patterns, since anomalies are derived differently in each case, requiring careful consideration for a meaningful quantitative comparisons. In the following section, we will discuss how the results for WACCM6-MAM compare between fully coupled and nudged simulations.

3.3 Further exploration of radiative forcings in WACCM6-MAM

In this section we show a comparison of radiative forcing in the three different experiments, which only WACCM6-MAM conducted in full. This comparison provides useful insights to the different ways to define the radiative impacts of the volcanic cloud between RF, ERF and IRF.

Figure 5 Time series of zonal mean radiative forcing under clear-sky conditions for SO2andH2O experiment in WACCM6-MAM, at TOA (a–c), TROP (d–f), and SURF (g–i). Each column represents a different model setup: first and second columns are free running experiments with either coupled ocean or climatological sea surface temperatures and sea ice (Coupled and Fixed-SST, respectively), third column indicates nudged simulations (Nudged). Results from Fixed-SST and Coupled are ensemble means over all 30 members. Gray areas indicate regions where the differences are not statistically significant at the 5 % level based on Student's t-test.

Figure 6 Time series of net radiative forcing under clear-sky conditions between 60° S and the equator for SO2andH2O experiment in WACCM6-MAM at TOA (a), TROP (b), and SURF (c). Solid lines indicate the results form nudged simulations, dashed lines from the free running simulations with fixed climatological SSTs and sea ice and dotted lines from the coupled experiment. Shading indicates the monthly internal variability, calculated as 1 standard deviation over 10-year simulations using a 30-member ensemble from the unperturbed NoVolc experiment.

Figure 5 shows that, across all experiments and at the three atmospheric levels, the simulated clear-sky radiative forcing is negative and locally statistically significant in the Southern Hemisphere during the first two years following the eruption. The sign, as well as the spatial and temporal evolution, is consistent across the different model configurations and atmospheric levels considered. However, when including the atmospheric temperature adjustments in Fixed-SST and the ocean response in Coupled, the radiative forcing from aerosols and water vapor is substantially smaller than in the nudged simulation, particularly at TOA. This difference is more evident when averaging over the Southern Hemisphere (60° S to the equator), as shown in Fig. 6.

Table 4 Radiative forcing (in W m⁻²) at the Top of Atmosphere under clear-sky conditions in WACCM6-MAM. Values are averaged globally (90° S–N) and between 60° S and the equator, and calculated over 2022–2023.

In Fig. 6a, the clear-sky IRF at TOA during the first 2 years after the eruption is −0.43 W m⁻² in the nudged simulations, whereas in the free running experiments it is −0.27 W m⁻² in both the coupled ocean and atmosphere-only cases. Although smaller than in the nudged simulations, these responses remain outside the range of natural variability, which is estimated in Fig. 6 as one standard deviation from the control ensemble (NoVolc). In general, while most of the results presented here are for the SH only, when considering global mean the values are reduced by approximately half, reaching a value of −0.09 W m⁻² at the TOA for the Coupled experiment in 2022–2023, which falls within the range of natural variability (see Fig. A6 and Table 4).

The reduced net radiative forcing relative to the nudged case arises from the inclusion of temperature and ozone adjustments. Indeed, differences in the simulated stratospheric AOD changes among the experiments are negligible (Fig. A7a–c), consistent with previous studies showing that atmospheric nudging does not necessarily improve the representation of the residual stratospheric circulation (Chrysanthou et al., 2019). Instead, the response is largely dominated by gas-driven variability, which is strongest in the Coupled experiment and includes a late-period negative radiative forcing in 2025 that is not present in the other simulations (Fig. 6a and b).

Figure 7 Time series of net radiative forcing at TOA for the SO2andH2O fully coupled experiment in WACCM6-MAM. RF is averaged between 60° S and the equator and calculated in Clear-sky and All-Sky conditions (a and b, respectively). Different colors indicate the total radiative forcing (gray), the aerosol-only contribution (green), and the gas-only contribution (red). Red line styles distinguish the net forcing (solid), the shortwave (SW) component (dashed), and the longwave (LW) component (dotted).

The separation of aerosol–radiation and gas–radiation contributions to the total radiative forcing (Fig. 7) shows that the negative TOA forcing during the first two years following the eruption is primarily driven by the aerosol contribution from Hunga, partially offset by a positive contribution from changes in gases. This behavior contrasts with the nudged simulations, where the gas–radiation interaction produces a negative forcing (Fig. 3). In the free-running simulations, the gas contribution is strongly influenced by natural variability (red line in Fig. 7), particularly for all-sky radiative forcing, which limits the detectability of the forced response without large ensemble sizes (Fig. A8).

In the second half of 2025, as the tropical stratospheric AOD returns to background levels (Fig. A7a–c), only the Coupled experiment exhibits a significant negative radiative forcing in the tropics at both TOA and TROP, with magnitudes comparable to those observed during the first post-eruption year (Fig. 6a and b), while the forcing at SURF becomes positive. This late-period forcing in the Coupled experiment is primarily driven by gas–radiation interactions, with significant negative clear-sky values emerging in 2025 (Fig. 7a, red lines).

Figure 8 Time series of net, shortwave (SW), and longwave (LW) gas–radiation interaction RF (sAOD; a–c), tropical H₂O concentrations and tropical temperature at 10 hPa (30° S–30° N; d–f, in ppm and K, respectively), and tropical O₃ concentrations and tropical temperature at 70 hPa (30° S–30° N; g–i, in ppm and K, respectively). Changes are computed as the difference between the SO2andH2O and NoVolc experiments for the (a, d, g) Coupled, (b, e, h) Fixed-SST, and (c, f, i) Nudged configurations.

The gas radiative forcing results not only from Hunga-induced changes in stratospheric water vapor (which are similar across the three experiments; Fig. A7d–f) but also from changes in stratospheric temperature and ozone (Fig. A7g–i), as well as associated dynamical adjustments that modulate both, particularly in the lower stratosphere. These contributions are further explored in Fig. 8. The gas–radiation interaction is dominated by the LW component (Fig. 8a–c), with differences among the atmospheric configurations driven primarily by temperature and dynamical responses: Zhuo et al. (2025) has shown an upper-level cooling due to the water vapor and a lower stratospheric level warming from LW absorption from the sulfate aerosols in the free-running simulations.

During the first two years following the eruption, the TOA gas radiative forcing is positive in the Coupled simulation, more variable in Fixed-SST, and negative in the Nudged case (Fig. 8a–c). This behavior reflects the combined effects of increased water vapor in the middle stratosphere and reduced ozone in the lower stratosphere. When temperature and dynamical changes are excluded, as in the Nudged simulation, the water vapor increase and ozone decrease (black curves in Fig. 8f and i) enhance outgoing LW radiation at TOA, resulting in a negative forcing. In contrast, when temperature adjustments are included, as in the Coupled and Fixed-SST simulations, the cooling at higher altitudes where water vapor anomalies peak (black curves for water vapor and red curves for temperature in Fig. 8d and e) reduces outgoing LW radiation, yielding a positive TOA forcing during the first two post-eruption years. Furthermore, the Coupled shows a distinct evolution in the tropical lower stratosphere, characterized by increased ozone and warming during 2022–2023, followed by decreased ozone and cooling in 2025 (Fig. 8g). Because ozone in the lower stratosphere acts as a greenhouse gas, these changes contribute to an additional positive radiative forcing during the first two years and a negative forcing in 2025, when water vapor anomalies have largely returned to background levels. The close correlation of changes in ozone and temperatures in this region, not seen in Fixed-SST (Fig. 8h), is strongly indicative of their dynamical origin, suggesting an associated decrease in tropical upwelling in 2022–2023 and increase in 2025.

As discussed in Bednarz et al. (2026), the coupled ocean WACCM6-MAM simulations shows a significant modulation of the El Niño–Southern Oscillation (ENSO) variability by the eruption, with La Nina like response in 2022–2023 and an El Nino like response in 2025. In general, ENSO is an important driver of interannual variability in tropical upwelling, which in turn modulates lower-stratospheric temperatures and ozone (Randel et al., 2009), and can therefore exert a substantial influence on the overall radiative forcing.

4 Conclusions

In this work we have provided the first multi-model analyses of the radiative impacts of the Hunga volcanic eruption, which co-injected large amounts of water vapor and a small-to-moderate amount of SO₂ into the stratosphere. Our multi-model results confirm previous analyses from Zhu et al. (2022), Stenchikov et al. (2025) and Sellitto et al. (2025) which indicated a potential net negative forcing from the volcanic cloud, due to the formation of a persistent layer of stratospheric sulfate aerosol that rapidly grew to optically efficient sizes, estimated at approximately 0.4 µm for sulfate particles, enhancing shortwave scattering relative to background conditions (English et al., 2013; Asher et al., 2025).

In particular, our analysis indicates a global mean effective radiative forcing at the top of the atmosphere of −0.14 ± 0.10 W m⁻², based on the multi-model mean from free-running simulations with fixed sea surface temperatures. If the averaging is restricted to the Southern Hemisphere, these values nearly double, highlighting the hemispheric asymmetry in the distribution of the aerosol cloud and its radiative impact. When coupled with the ocean, the global forcing is smaller and more noisy (−0.09 ± 0.16 W m⁻²), due to the eruption's significant modulation of ENSO variability, triggering a La Niña–like response in 2022–2023 and an El Niño–like response in 2025 (Bednarz et al., 2026), which, in turn, impacts tropical upwelling and alters lower stratospheric ozone in the tropics. However, the predominant effect arises from the sulfate aerosols (−0.18 ± 0.02 W m⁻², IRF from the multi-model mean in the Nudged SO2only) and only a marginal contribution from the water vapor. The two methods used to estimate the IRF from aerosol–radiation and gas–radiation interactions, one simulating separate injections of SO₂ and water vapor, and the other employing a double radiation call for their co-injection, reveal biases in IRF calculations. Our results provide useful insight that can be used to inform future climate assessments (Forster et al., 2025) that aim to identify the contributions of the single natural and anthropogenic factors to global radiative imbalance and temperatures.

It is important to note that in the simulations analyzed in this work, a value of 0.5 Tg of SO₂ was used (Zhu et al., 2025). Analyses by Sellitto et al. (2024), using the Infrared Atmospheric Sounding Interferometer (IASI), however, suggest that the overall SO₂ burden from the Hunga eruption was larger than what Carn et al. (2022) estimated using UV measurements from the Ozone Monitoring Instrument (OMI) on NASA's Aura satellite, with a lower limit of 1.0 Tg. This suggests that the results presented here might be underestimated, if higher retrieval estimates for the sulfate burden were confirmed, particularly in light of the uncertainties in stratospheric AOD retrievals highlighted by comparisons between the GloSSAC and OMPS-LP datasets in Zhuo et al. (2025).

In general, AOD is often assumed to scale linearly with radiative forcing, especially in the stratosphere, meaning that a doubling of AOD would typically result in a doubling of negative radiative forcing. However, this relationship depends not only on the total aerosol burden but also on the particle size distribution (Sellitto et al., 2025), as AOD is highly sensitive to particle size Consequently, an increased aerosol burden does not necessarily translate into a proportionally larger AOD or radiative forcing.

Our analyses further indicate that stratospheric water vapor plays a significant role in the initial growth of sulfate aerosols and thus in modulating the aerosol–radiation interaction. Results from nudged WACCM6-MAM simulations show that, during the first 6 months following the eruption, the aerosol–radiation interaction IRF in the Southern Hemisphere would be over 10 % smaller (−0.34 W m⁻² vs −0.39 W m⁻²) if the co-injection of water vapor is not included.

As a result, even if the sulfate aerosol burden were doubled, non-linearities in the AOD–forcing relationship would likely be less than double. Ultimately, the enhanced particle scattering efficiency and longer atmospheric lifetime following the Hunga eruption – attributed to its higher injection altitude and increased co-emission of water vapor (Li et al., 2024) – help explain why an SO₂ injection at least 20 times smaller than Pinatubo (10–12 Tg;  Ukhov et al., 2023) resulted in an stratospheric AOD only about 10 times lower, yet still capable of exerting a small but non-negligible impact on Earth’s radiative balance over the following year.

Appendix A: Figures and tables

Figure A1 Time series of zonal mean radiative forcing at the tropopause (TROP) under clear-sky conditions from three models: WACCM6-CARMA (a–c), WACCM6-MAM (d–f), and MIROC-CHASER (g–i). Each column corresponds to a different perturbation scenario from the nudged experiment: the first column shows SO2andH2O, the second column SO2only, and the third column H2Oonly.

Figure A2 Time series of zonal mean radiative forcing at the surface (SURF) under clear-sky conditions from five models: WACCM6-MAM (a–c), WACCM6-CARMA (d–f), MIROC-CHASER (g–i), CAM5-CARMA (j–l), and UKESM (m). Each column corresponds to a different perturbation scenario from the nudged experiment: the first column shows SO2andH2O, the second column SO2only, and the third column H2Oonly.

Figure A3 Time series of zonal mean stratospheric aerosol optical depth anomalies (perturbed minus control) from five models: WACCM6-MAM (a, f), WACCM6-CARMA (b, g), MIROC-CHASER (c, h), CAM5-CARMA (d, i), and UKESM (e). The first two rows correspond to perturbation scenarios from SO2andH2O, and the second row shows SO2only. The third row includes regional averages of the same quantity for two latitude bands: 90° S–N (j), 60° S–equator (k) and 60–90° S (l).

Figure A4 Time series of net radiative forcing under clear-sky conditions at three atmospheric levels: TOA (a, d, g), TROP (b, e, h), and SURF (c, f, i). Each column corresponds to a different perturbation scenario from the nudged experiment: the first column shows SO2andH2O, the second column SO2only, and the third column H2Oonly. The black dashed line represents the multi-model mean, each models is represented with a different color.

Figure A5 Time series of effective radius averaged over 30° S and the equator for (a, d, g) WACCM6-MAM, (b, e, h) WACCM6-CARMA, and (c, f, i) MIROC-CHASER. Each row corresponds to a different scenario from the nudged experiment: the first column shows SO2andH2O, the second column SO2only, and the third column NoVolc.

Figure A6 Time series of global net radiative forcing under clear-sky conditions for SO2andH2O experiment in WACCM6-MAM at TOA (a), TROP (b), and SURF (c). Solid lines indicate the results form nudged simulation (Nudged), dashed lines from the fixed SST simulation (Fixed-SST) and dotted lines from the coupled experiment (Coupled). Shading indicates the monthly internal variability, calculated as 1 standard deviation over 10-year simulations using a 30-member ensemble from the unperturbed NoVolc experiment.

Figure A7 Time series of zonal mean changes in stratospheric aerosol optical depth (sAOD; a–c), tropical stratospheric H₂O concentrations (30° S–30° N; d–f, in ppm), tropical stratospheric O₃ concentrations (30° S–30° N; g–i, in ppm), and tropical stratospheric temperature(30° S–30° N; j–l, in K). Changes are computed as the difference between the SO2andH2O and NoVolc experiments for the (a, d, g, j) Coupled, (b, e, h, k) Fixed-SST, and (c, f, i, l) Nudged configurations. Gray areas indicate regions where the differences are not statistically significant at the 5 % level based on Student's t-test.

Figure A8 Detectability of the 2022–2023 average RF SURF response in the Coupled simulations (left) and the ERF response in the Fixed-SST simulations (right) for clear-sky (top) and all-sky (bottom) conditions. Results obtained by randomly subsampling each ensemble with replacement to obtain 2000 artificial ensembles each of different ensemble size. Black lines denote the mean response, and blue and red lines indicate the ± 2 standard deviation and the max/min ranges, respectively, of the possible responses. The results show that a 2-year mean clear-sky RF/ERF response is detectable already with a few ensemble members (∼ 5) but requires more members to constrain the magnitude with any confidence. In contrast, the uncertainty of the all-sky response is substantially larger, particularly in the Coupled experiment but also in the Fixed-SST experiment, where even 30 ensemble members may be insufficient to determine even the sign of the response with any confidence.

Table A1 Global effective radiative forcing (in W m⁻²) at the Top of Atmosphere under clear-sky conditions.

Table A2 Effective radiative forcing (in W m⁻²) at the Tropopause under clear-sky conditions. Forcing is averaged between 60° S and the equator.

Table A3 Effective radiative forcing (in W m⁻²) at the Surface under clear-sky conditions. Forcing is averaged between 60° S and the equator.