How well are aerosol–cloud interactions represented in climate models? – Part 2: Isolating the aerosol impact on clouds following the 2014–2015 Holuhraun eruption

Aerosols significantly influence Earth's radiative balance, yet considerable uncertainty exists in the underpinning mechanisms, particularly those involving clouds. Aerosol–cloud interactions (ACIs) are the most uncertain element in anthropogenic radiative forcing, hampering our ability to constrain Earth's climate sensitivity and understand future climate change. The 2014–2015 Holuhraun volcanic eruption in Iceland released sulfur dioxide (SO₂) into the lower troposphere on a level comparable to continental-scale emissions. The resultant volcanic plume across an often near-pristine region of the northern North Atlantic Ocean presents an ideal opportunistic experiment to explore ACI representation within general circulation models (GCMs). We present Part 2 of a two-part AeroCom (Aerosol Comparisons between Observations and Models) Phase III inter-model comparison study that utilises satellite remote sensing observations to assess modelled cloud responses to the Holuhraun attributed volcanic aerosol within eight state-of-the-art GCMs during September and October 2014. We isolate the aerosol effect from meteorological variability and find that the GCMs – particularly their multi-model ensemble response – adeptly capture the observed cloud microphysical changes associated with the ACI first indirect effect (i.e. Twomey effect). Meanwhile, a clear divergence exists in the GCM responses of large-scale cloud properties, namely cloud liquid water content, expected from the precipitation suppression mechanism of the ACI second indirect effect (i.e. rapid adjustments). We attribute this to limitations and differences in their autoconversion schemes under high aerosol loading, specifically in sub-grid-variability representations. Finally, our multi-model ensemble estimates that Holuhraun had a global radiative forcing of −0.11 ± 0.04 Wm⁻² across September and October 2014.

1 Introduction

Aerosols have a major influence on the Earth's energy budget through their interactions with solar and terrestrial radiation via direct and indirect mechanisms. The direct mechanism – termed aerosol–radiation interactions – describes the scattering and absorption of radiation by the aerosol itself (e.g. Bellouin et al., 2020; Myhre et al., 2013), whilst the indirect mechanism – known as aerosol–cloud interactions (ACIs) – centres on changes to cloud properties caused by aerosols via their role as cloud condensation nuclei (CCN) (e.g. Bellouin et al., 2020; Fan et al., 2016). Overall, aerosols exert a negative radiative forcing (RF) on the Earth helping offset a portion of the warming from increased greenhouse gas emissions, yet the magnitude of this key effect continues to be a major source of uncertainty in anthropogenic climate change (Forster et al., 2021; Gryspeerdt et al., 2020; Watson-Parris and Smith, 2022). This uncertainty stems predominantly from ACIs, meaning it is of paramount importance that we improve our knowledge of these cloud-mediated processes to improve future climate estimates.

Aerosols prompt cloud modifications through a causal network of events (e.g. Haywood and Boucher, 2000; Fan et al., 2016). For liquid-only clouds, added aerosol can serve as additional CCN, which increases cloud droplet number concentration (N_d) (Twomey, 1974). Holding cloud liquid water content constant (cloud liquid water path, LWP), an increase in N_d leads to a decrease in cloud droplet size (cloud droplet effective radius, r_e), causing an enhancement in cloud albedo (Twomey, 1977). This chain of events is referred to as the “first indirect effect” or the “Twomey effect”. Furthermore, smaller cloud droplets decrease the efficiency of collision-coalescence processes delaying the formation of precipitation. Consequently, liquid clouds polluted by aerosol may have longer lifetimes and/or greater cloud fraction (CF) (Albrecht, 1989) and increased depth (Pincus and Baker, 1994), all of which act to increase LWP and further enhance cloud albedo. This subsequent chain of events has historically been referred to as the “second indirect effect”, although now further aerosol-induced cloud adjustments are often captured under this term too. Such adjustments include those in non-precipitating clouds whereby the aerosol-induced reduction in r_e increases evaporation and decreases sedimentation, causing feedbacks that help accelerate entrainment and deplete LWP (Ackerman et al., 2004; Bretherton et al., 2007; Hill et al., 2009; Small et al., 2009). For mixed-phase and ice-only clouds, additional cloud modification processes exist (e.g. Bellouin et al., 2020; Fan et al., 2016; Forster et al., 2021). The myriad of mechanisms underpinning ACIs – each with their own dependency on conditions both meteorological (e.g. atmospheric stability, humidity, temperature) and environmental (e.g. aerosol background concentrations, marine versus land region) – is testament to how challenging constraining ACI uncertainty is.

To alleviate this complexity, studies can focus on aerosol perturbations to systems where the meteorology and environment are well understood. Known as “opportunistic experiments”, these instances include industrial plumes, ship tracks, wildfires, regulatory changes, and volcanic eruptions (Christensen et al., 2022). A notable example of the latter is the Holuhraun eruption, an effusive eruption that occurred continuously between 31 August 2014 and 27 February 2015 in the Bárðarbunga volcanic system in Iceland (64.85° N, 16.83° W) (Gislason et al., 2015; Pedersen et al., 2017). Characterised by its non-explosive nature, Holuhraun released an estimated 9.6–11.8 Tg of sulfur dioxide (SO₂) (Gislason et al., 2015; Pfeffer et al., 2018) – approximately one-tenth of current global annual anthropogenic SO₂ emissions (Aas et al., 2019; Szopa et al., 2021) – into the lower troposphere (Carboni et al., 2019; Flower and Kahn, 2020; Pfeffer et al., 2018). These SO₂ emissions subsequently oxidised to sulfate aerosol (SO${}_{\mathrm{4}}^{\mathrm{2}-}$), leading to the formation of a vast aerosol plume. Such widespread pollution to an often near-pristine region of the northern North Atlantic Ocean over a 6-month duration has made Holuhraun a focal point in studying ACIs at the climatic scale.

Previous Holuhraun studies have provided valuable insight into ACIs through a variety of approaches. For example, Malavelle et al. (2017) and Zoëga et al. (2025) use general circulation models (GCMs) to generate climatologies within the North Atlantic Ocean and Arctic Ocean respectively, enabling the volcanic aerosol effect on cloud properties to be disentangled from meteorological variability. Both studies find that GCMs simulate a decrease in r_e during the months following the eruption, yet their LWP responses range from negligible change to a strong increase. Alternatively, Haghighatnasab et al. (2022) and Peace et al. (2024) use an “in-plume versus out-of-plume” approach to isolate the aerosol-induced cloud impacts during September 2014. The studies find increases in N_d and decreases in r_e inside the plume compared to outside, whereas the in-plume changes to LWP are mixed and hard to isolate. Moreover, Zoëga et al. (2025) use a GCM to explore the cloud response sensitivity to Holuhraun with respect to eruption season and size of emissions, noting a stronger response occurs during spring and summer and a plateauing of the response with increasing emissions. McCoy and Hartmann (2015) perform an entirely observation-based study, noting a decrease in r_e post-eruption, yet no appreciable changes in LWP or CF. Additionally, Chen et al. (2022) trained a machine learning model to produce a “counterfactual” satellite remote sensing representation of the region absent of Holuhraun emissions, again finding that N_d increases and r_e decreases due to the eruption, with minimal changes to LWP. Interestingly, Chen et al. (2022) propose that the additional aerosol prompted a 10 % increase in cloud cover, a result not found in other Holuhraun studies exploring this cloud property.

Here we build on this established set of works by presenting Part 2 of a two-part AeroCom (Aerosol Comparisons between Observations and Models) Phase III inter-model comparison study of the Holuhraun plume and its interactions with clouds. In Part 1, the spatial and chemical evolution of the volcanic plume was assessed (Jordan et al., 2024). Differences in the secondary SO${}_{\mathrm{4}}^{\mathrm{2}-}$ aerosol production amongst the GCMs, as well as with observations, were noted, yet overall the modelled representations of the Holuhraun plume were deemed sufficient to explore the impacts of the eruption on ACIs in the region. Here we follow on from Part 1 and assess the ACI representations from eight state-of-the-art GCMs against satellite remote sensing observations. Here we focus on stratocumulus clouds over an often near-pristine marine region (i.e. minimal anthropogenic influence). We narrow our attention to September and October 2014 as these months offer the most favourable conditions for isolating the aerosol signal relative to the later months of the eruption. This is due to (a) more reliable satellite retrievals, (b) peak SO₂ emissions, (c) a well-defined volcanic plume with minimal plume dilution, and (d) a plethora of studies providing insights into the conditions of this period. We compare model analyses and observations to identify differences in ACI representations, seeking to understand the point at which the models depart from the observed ACI casual chain. We conclude with an updated multi-model ensemble forcing estimate of the Holuhraun eruption.

2 Methodology

Here we briefly introduce the experimental set-up and ACI relevant components of the eight GCMs, provide an overview of the four remote sensing products used to assess the GCMs, outline the theoretical framework used to disentangle the aerosol effect from meteorological variability, and describe the identification of regions subject to significant SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations attributed primarily to Holuhraun emissions.

Kirkevåg et al. (2018)Liu et al. (2016)Neale et al. (2012)Michou et al. (2015, 2020)Séférian et al. (2019)Neubauer et al. (2019); Tegen et al. (2019)Dietlicher et al. (2018)Kokkola et al. (2018)Walters et al. (2019)Mulcahy et al. (2020)Mulcahy et al. (2020)

Table 1 Models used in this study. Aerosol module: name of the aerosol module with type given in brackets. Stratiform cloud microphysics: name of the stratiform cloud microphysics scheme (MG1.5 – Gettelman and Morrison, 2015; Morrison and Gettelman, 2008; Lopez – Lopez, 2002; Lohmann – Lohmann et al., 2007; Lohmann and Hoose, 2009; Lohmann and Neubauer, 2018; P3 – Dietlicher et al., 2018; WB –Wilson and Ballard, 1999). Activation: name of the cloud droplet activation scheme (ARG – Abdul-Razzak and Ghan, 2000; Menon – Menon et al., 2002). Autoconversion: name of the autoconversion parameterisation (KK – Khairoutdinov and Kogan, 2000; Smith – Smith, 1990). ACIs: aerosol indirect effects represented. Lat. × long.: atmospheric grid resolution. Levs.: number of vertical levels. References: key references.

¹ Production-tagged: size-resolving through offline lookup tables. ² Refers explicitly to an absence of aerosol-induced precipitation suppression effects on large-scale cloud properties. ³ Aerosol indirect effects simulated in liquid clouds only. ⁴ With significant developments in the warm rain microphysics by Boutle et al. (2014).

2.1 General circulation models

The relevant features of the eight GCMs that participated in Part 2 of this inter-model comparison study are listed in Table 1. Among the models, three types of aerosol modules are employed – modal, sectional, and production-tagged – alongside five distinct stratiform cloud microphysics schemes. Of the eight GCMs, some are based on the same core model. ECHAM6 appears in three configurations, each with a different combination of aerosol module and stratiform cloud microphysics scheme, whilst two versions of UKESM1 are assessed, one with and one without boundary layer nucleation (BLN). In terms of aerosol activation to cloud droplets, CNRM-ESM2-1 uses the empirical-based parameterisation of Menon et al. (2002), whereas all the other models adopt the Köhler-theory-based parameterisation of Abdul-Razzak and Ghan (2000). For CNRM-ESM2-1, HadGEM3, and the UKESM1 variants, N_d is calculated diagnostically, whereas a prognostic approach is taken in the remaining models. In each model, the calculated N_d from these activation parameterisations is passed to the radiation scheme and used in the calculation of r_e, which in turn is used to calculate cloud albedo, thereby enabling the simulation of the ACI first indirect effect. The autoconversion parameterisation of Khairoutdinov and Kogan (2000), hereafter KK2000, is adopted in seven of the GCMs, albeit with different approaches to account for sub-grid variability (see Sect. A). As KK2000 has an inverse power-law dependence on N_d, an increase in the number of (smaller) cloud droplets acts to inhibit the rate of autoconversion. Hence, the models using KK2000 are able to represent the precipitation suppression mechanism of the ACI second indirect effect. The exception is CNRM-ESM2-1, which instead adopts the autoconversion parameterisation described in Smith (1990) – a scheme not dependent on N_d, thus preventing aerosol influencing large-scale cloud properties via precipitation suppression. Finally, the additional ACI mechanism of enhanced entrainment evaporation has been shown to be possible in GCMs (e.g.  Mülmenstädt et al., 2024), yet the effects are found to be negligible and so will not be considered here.

Performed in their atmosphere-only configurations using prescribed sea surface temperature and sea ice fraction (“AMIP-style”), each model provided simulations with and without Holuhraun SO₂ emissions – the former for 2014 only and the latter for 2002–2014 that acts as a long-term control. To reduce model internal variability and to obtain a model meteorology that closely resembles the observed meteorology during the eruption, horizontal winds are constrained (“nudged”) towards ERA-Interim reanalysis (Dee et al., 2011) on a 6 h timescale. The Holuhraun simulations distribute the volcanic SO₂ equally between 0.8 and 3 km within the grid cell containing the eruption vent following the emissions profile shown in Table 2, which is based on empirical estimates by Thordarson and Hartley (2015). Both Holuhraun and control simulations contain additional background SO₂ emissions from anthropogenic and natural sources, including passive degassing volcanoes, in accordance with AeroCom Phase III guidelines. Where possible, in-cloud diagnostics directly outputted from the models are used (i.e. model performs necessary calculations during simulation), rather than dividing grid cell mean values by mean CF post-simulation. All model output is regridded to a regular 1.0° × 1.0° latitude–longitude grid using linear interpolation, aside from precipitation diagnostics which use first-order conservative interpolation to preserve precipitation totals.

Table 2 Sulfur dioxide (SO₂) emissions profile used to represent the Holuhraun eruption. Emissions are prescribed in the grid cell containing the eruption vent (64.85° N, 16.83° W) and follow empirical estimates by Thordarson and Hartley (2015).

2.2 Satellite observations

2.2.1 MODIS

This study uses the Moderate Resolution Imaging Spectroradiometer (MODIS) MCD06COSP version 6.2.0 Level-3 product (Pincus et al., 2023) to quantify the volcanic impact on cloud properties. The MCD06COSP dataset combines observations from MODIS instruments on board the Aqua and Terra satellites obtained using the 3.7 µm Cloud Optical Properties Retrieval Algorithm (Platnick et al., 2017). The Level-3 data are outputted at daily and monthly timescales to a regular 1.0° × 1.0° latitude–longitude grid having been sampled from pixel-scale (Level-2) data. These pixel-scale data estimate cloud properties for sunlight pixels (solar zenith angle <81.3731°) flagged as either “confidently” or “probably cloudy”. Cloud phase – liquid, ice, or undetermined – is decided at 1 km resolution following Marchant et al. (2016). The pixel-scale data are aggregated to daily Level-3 data, which themselves are aggregated to monthly data by weighting each day based on pixel count; this differs from the standard monthly MODIS product (MOD08_M3), which treats each day equally. Aqua and Terra satellites have a 16 d return period, so sampling within months is largely uniform yet reduced in the winter hemispheres due to limited illumination. We use the monthly mean product, except for all-sky LWP and N_d, which are calculated at a daily resolution before averaging to monthly means adopting the pixel count weighting above. All-sky LWP is calculated as the product of the in-cloud LWP (cloudy portions of observed region only) and liquid CF, whilst N_d is derived from liquid phase r_e and cloud optical depth using the Idealised Stratiform Boundary Layer Cloud (ISBL) model (Bennartz and Rausch, 2017; Quaas et al., 2006, 2008). The validity of the assumptions required for our N_d derivation is discussed at length in Grosvenor et al. (2018). To ensure only the most reliable retrievals are considered for estimating N_d, pixels are restricted using r_e and cloud optical depth bounds of 4–30 and 4–70 µm respectively (e.g. Chen et al., 2022; Haghighatnasab et al., 2022; Peace et al., 2024).

2.2.2 GPCP

To assess precipitation across the northern North Atlantic Ocean during the Holuhraun eruption, we use the Global Precipitation Climatology Project (GPCP) version 3.2 product (Huffman et al., 2023). The GPCP dataset assimilates satellite remote sensing data (e.g. low-orbit passive-microwave sensors, geostationary infrared sensors, and sounders) and ground-based gauge analyses to provide global surface precipitation estimates on a regular 0.5° × 0.5° latitude–longitude grid. Here we utilise the monthly product regridded to a 1.0° × 1.0° resolution using first-order conservative interpolation to preserve precipitation totals.

2.2.3 CERES-EBAF

To evaluate the influence of Holuhraun emissions on top-of-atmosphere (ToA) radiative fluxes, we use the Clouds and the Earth's Radiant Energy System (CERES) – Energy Balanced and Filled (EBAF) product (Loeb et al., 2018; Kato et al., 2018) – specifically the ToA Edition 4.2 dataset. The CERES-EBAF product contains monthly mean longwave (LW), shortwave (SW), and net radiative fluxes at ToA under all-sky and clear-sky conditions outputted to a regular 1.0° × 1.0° latitude–longitude grid. The dataset combines observations from narrow field-of-view scanning radiometer instruments and imagers on board polar-orbiting Aqua, Terra, Suomi National Polar-Orbiting Partnership (SNPP), and NOAA-20 satellites, along with additional geostationary imagers that provide data between overpasses. The CERES-EBAF product adjusts ToA SW and LW radiative fluxes within their range of uncertainty to correct the discrepancy between the net energy imbalance observed at ToA and the heat storage within the Earth system (Loeb et al., 2009).

2.3 Separating aerosol and meteorological effects

In this study we adopt the simple theoretical framework used by Malavelle et al. (2017) to separate aerosol and meteorological effects on cloud properties. If the properties of cloud, c, are a function of aerosol, a, and meteorology, m, then – neglecting any interdependency between a and m – a change in c can be expressed as

$$\begin{array}{}\text{(1)}& \mathit{\delta }c=\mathit{\delta }a{\displaystyle \frac{\partial c}{\partial a}}+\mathit{\delta }m{\displaystyle \frac{\partial c}{\partial m}}.\end{array}$$

By combining the 2014 Holuhraun and long-term control simulations, we can use Eq. (1) to find the total change of a cloud property during the eruption, as well as isolating the change's aerosol and meteorological components.

2.3.1 Total effect

We estimate the total effect on a cloud property (i.e. Eq. 1) by subtracting the long-term control (NoHol_clim) from the 2014 simulation with the eruption (Hol₁₄). This anomaly is directly comparable to observations and is expressed succinctly as

$$\begin{array}{}\text{(2)}& \text{total effect}={\text{Hol}}_{\mathrm{14}}-{\text{NoHol}}_{\text{clim}}.\end{array}$$

Note that we remove the year 2014 from NoHol_clim to avoid double-counting/dilution of the meteorological variability.

2.3.2 Aerosol-only effect

As the models are nudged, meteorological differences between Hol₁₄ and the 2014 simulations without the eruption (NoHol₁₄) are negligible (i.e. δm≈0). For this special case, Eq. (1) approximates to

$$\begin{array}{}\text{(3)}& \mathit{\delta }c\approx \mathit{\delta }a{\displaystyle \frac{\partial c}{\partial a}}.\end{array}$$

Hence, we estimate the aerosol-only effect on a cloud property using

$$\begin{array}{}\text{(4)}& \text{aerosol-only effect}={\text{Hol}}_{\mathrm{14}}-{\text{NoHol}}_{\mathrm{14}}.\end{array}$$

2.3.3 Meteorology-only effect

With background aerosol largely the same for each year within a particular model, differences in aerosol between NoHol₁₄ and NoHol_clim are negligible (i.e. δa≈0). In this instance, Eq. (1) approximates to

$$\begin{array}{}\text{(5)}& \mathit{\delta }c\approx \mathit{\delta }m{\displaystyle \frac{\partial c}{\partial m}}.\end{array}$$

Hence, we estimate the meteorology-only effect on a cloud property using

$$\begin{array}{}\text{(6)}& \text{meteorology-only effect}={\text{NoHol}}_{\mathrm{14}}-{\text{NoHol}}_{\text{clim}}.\end{array}$$

Figure 1 The multi-model ensemble mean perturbation in sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$) column load for (a) September and (b) October 2014. Perturbation depicted is the aerosol-only anomaly with meteorological variability excluded (i.e. Hol₁₄ – NoHol₁₄) and is expressed in Dobson units (DU). Predominantly volcanically polluted (PVP) regions are defined over ocean areas where the SO${}_{\mathrm{4}}^{\mathrm{2}-}$ column load anomaly exceeds 0.2 DU and anthropogenic aerosol load is low (see main text). These PVP regions are outlined by dotted lines, with the corresponding spatial mean listed above.

2.4 Predominantly volcanically polluted regions

This study focuses on a northern region of the North Atlantic Ocean, which is distanced from the anthropogenic aerosol sources of continental Europe. Here low aerosol optical depth values ranging from 0.04 to 0.15 – indicative of “near-pristine to clean” conditions – feature consistently in September–November climatologies derived from a range of datasets (e.g. global aerosol reanalyses – Xian et al., 2024; global atmospheric models – Gliß et al., 2021; Li et al., 2022; and remote sensing instruments – Bevan et al., 2012; Remer et al., 2008). Clouds in such areas are likely more susceptible to changes in aerosol concentrations making volcanic impacts on ACIs more apparent and easier to isolate.

In the absence of suitable SO${}_{\mathrm{4}}^{\mathrm{2}-}$ observations, and knowing the models capture the spatial and chemical evolution of the plume with sufficient fidelity (see Part 1, Jordan et al., 2024), we use the modelled SO${}_{\mathrm{4}}^{\mathrm{2}-}$ column load to identify areas subject to high aerosol loading due to the eruption. As the background aerosol levels in comparison are low, we denote these areas predominantly volcanically polluted (PVP) regions. We avoid using SO₂ to distinguish PVP regions due to limitations in assuming the co-existence of SO₂ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ including divergent spatial dispersions, time lag in SO₂-to-SO${}_{\mathrm{4}}^{\mathrm{2}-}$ conversion, and differing deposition rates. The multi-model ensemble SO${}_{\mathrm{4}}^{\mathrm{2}-}$ column load aerosol-only anomaly (i.e. Hol₁₄ – NoHol₁₄) for September and October 2014 is shown in Fig. 1. The additional SO₂ emissions from Holuhraun clearly increase the SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations within the region, more so in September when the prescribed SO₂ emission rate is higher. The added aerosol loading is not uniformly distributed due to each month's differing meteorological conditions. To identify the PVP regions, we mask the grid cells over land, as well as grid cells with SO${}_{\mathrm{4}}^{\mathrm{2}-}$ column load anomalies below 0.2 Dobson units (DU). The former removes areas likely influenced by anthropogenic pollution, whilst the latter helps ensure a sufficient aerosol concentration to prompt ACI responses.

Nevertheless, despite the isolated nature of the perturbed region, continental pollution may be transported into the area under specific meteorological conditions. Peace et al. (2024) use air mass back-trajectory modelling to show that between 15–21 September, anomalous easterly winds brought anthropogenic pollution to the area south of Iceland, mixing with the aerosol load introduced by Holuhraun and hence diluting the volcanic influence there. Malavelle et al. (2017) further show that the anthropogenic pollution is the dominant cause of r_e perturbations in this area of mixed pollution, meaning any volcanic signal is likely weak. To keep the pollution in the September PVP region ∼ 100 % volcanic, we exclude areas below 62° N. Both monthly PVP regions and their associated multi-model ensemble SO${}_{\mathrm{4}}^{\mathrm{2}-}$ loading are outlined in Fig. 1. Unless otherwise stated, all values hereafter refer to these PVP regions and not – as is often the case in other Holuhraun studies (e.g. Chen et al., 2022; Malavelle et al., 2017) – the entire domain. Using PVP regions, coupled with the framework laid out in Sect. 2.3, will help attribute any cloud modifications found in this study to volcanic emissions.

Figure 2 Monthly mean anomalies in cloud droplet effective radius (r_e) at cloud top for September 2014 from (a) MODIS instruments on board Aqua and Terra satellites, (b) the multi-model ensemble, and (c–h) individual models. Anomalies depicted are the total effect, so they include both aerosol and meteorological components (i.e. Hol₁₄ – NoHol_clim). The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Stippling highlights grid cells with null hypothesis rejections based on applying the False Discovery Method (FDR) at a 10 % control level (see main text). Hatched areas indicate missing data. Note that the total effect on r_e at cloud top cannot be calculated for ECHAM6-HAM-P3 and ECHAM6-SALSA from the output provided to this experiment.

Figure 3 Disentanglement of the aerosol and meteorological effects on (a) cloud droplet number concentration (N_d) and (b) cloud droplet effective radius (r_e) at cloud top within the predominantly volcanically polluted (PVP) region for September 2014. Total, aerosol-only, and meteorology-only effects are depicted by green – no pattern, red – minor diagonal, and blue – major diagonal box plots respectively. Box plots extend to the 25th–75th percentiles with outer whiskers at 5th–95th. Black squares depict means. Green bounding and dashed lines extend the observed total effects across rows for visual comparison with the model responses. Climatological baselines are given in brackets. Note that solely the aerosol-only effect can be calculated for ECHAM6-HAM-P3 and ECHAM6-SALSA from the output provided to this experiment.

3 ACI first indirect effect

The total anomaly (i.e. Hol₁₄ – NoHol_clim) in cloud top r_e for September 2014 observed by MODIS is shown in Fig. 2a alongside the associated spatial mean of the PVP region. We evaluate the null hypothesis – that the increase in aerosol concentration had no effect on a cloud property – at each grid cell using a two-tailed Student t test. As these “local” hypothesis tests are mutually correlated, to avoid overstating their collective significance we apply the False Discovery Rate (FDR) method (Wilks, 2006, 2016) to control the overall expected proportion of false positives across the spatial domain to 10 %. Stippling highlights grid cells where the null hypothesis is rejected after applying the corrective FDR adjustment (i.e. the cloud property change is unlikely due to random variation). There is a clear decrease in r_e observed across the northern North Atlantic Ocean, particularly south-east of Iceland where anomalies can exceed −3.00 µm. The observed decrease in this area is greater than the PVP region and is likely due to the additional continental anthropogenic aerosol introduced by the meteorological conditions at the time (see Sect. 2.4). The associated modelled total anomalies in cloud top r_e are shown in Fig. 2b–h. All models capture the observed r_e anomalies within the PVP region well, notably CNRM-ESM2-1, which agrees to 2 decimal places. The multi-model ensemble mean is within 4 % of the MODIS mean, differing by only 0.07 Wm⁻². The GCMs do underestimate the observed decrease in r_e around the UK and Ireland where the transported continental anthropogenic aerosol exists, a discrepancy likely due to differences in the magnitude of background anthropogenic emissions between the real-world and simulated systems, rather than in the meteorological conditions given that the model runs are nudged. Evidence for a decrease in cloud top r_e during October is also observed, with the GCMs in good agreement (see Fig. C2). Similar spatial figures for cloud top N_d can be found in the Appendices (see Figs. B1 and C1).

A comprehensive disentanglement of the aerosol and meteorological effects on cloud top N_d and r_e for the PVP regions of September and October 2014 is shown in Figs. 3 and C3 respectively, with summary values provided in Tables D1 and D2. The box plots are generated from the monthly mean anomalies of the individual grid cells within the PVP region. For example, the September PVP region includes 579 grid cells, meaning 579 “local” values form the dataset used to construct the September box plots. The effect described by a box plot (i.e. total, aerosol-only, or meteorology-only) is dependent on the anomaly the values represent (see Sect. 2.3). MODIS retrievals depict an increase in N_d, which, with the aforementioned observed decrease in r_e, shows that an ACI first indirect effect initiated by Holuhraun aerosol features in the remote sensing record. The total effects modelled by the individual GCMs all follow the observed directional change for N_d and r_e. This, coupled with the component analysis showing that these changes are chiefly aerosol-induced, evidences the ability of the GCMs to capture the ACI first indirect effect within the PVP regions following the eruption, albeit with differing magnitudes. It is worth mentioning that, despite the varying strengths of the model responses, the multi-model ensemble is in good agreement with the observed cloud modifications highlighting the advantages of ensemble based techniques. Note that ECHAM6-HAM-P3 and ECHAM6-SALSA output provided to the experiment make it only possible to calculate the aerosol-only effect on N_d and r_e.

The variations in the ACI first indirect effect model representations can largely be explained by their configurations. For example, the strong response in N_d in ECHAM6-SALSA compared to the other two ECHAM6 models is likely due to the type of aerosol module employed. Sectional schemes, such as HAM-SALSA, better capture small particle growth following a pollution event than modal schemes, such as HAM-M7, due to their ability to resolve finer size distributions and nucleation events, generating more CCN and, subsequently, N_d (e.g.  Matsui and Mahowald, 2017; Mann et al., 2012; Saponaro et al., 2020). For highly polluted regions, as is the case here, these differences in microphysics can be exacerbated (Kokkola et al., 2018). In addition, the UKESM1 responses with and without BLN imply that including BLN leads to – somewhat counter-intuitively – lower N_d following the introduction of volcanic emissions. The rationale is that the newly nucleated particles from BLN are lofted vertically into the plume where they compete with the aerosol for condensible vapour, which hinders the growth of individual particles to CCN size, reducing the number available to form cloud droplets (i.e. clouds in the BLN simulations are less susceptible to increases in aerosol). Finally, despite similar increases in N_d, HadGEM3 simulates a considerably larger decrease in r_e than UKESM1 and UKESM1-BLN. This is expected due to improvements added to UKESM1 in aerosol processes, including to the cloud droplet spectral dispersion parameterisation (Mulcahy et al., 2018).

Figure 4 Monthly mean surface precipitation rates for September 2014 from the (a) Global Precipitation Climatology Project (GPCP), (b) the multi-model ensemble, and (d–k) individual models, as well as (c) the climatological September mean (2002–2014) derived from GPCP. The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Modelled precipitation rates are for the simulations including Holuhraun emission (i.e. Hol₁₄).

4 ACI second indirect effect

Delaying precipitation formation lies at the heart of the ACI second indirect effect, so it is useful to first assess precipitation totals as even a substantial aerosol perturbation, such as Holuhraun, cannot suppress precipitation in a non-precipitating cloud. Monthly mean surface precipitation rates for September 2014 are depicted in Fig. 4. Observational data from GPCP show that the PVP region is subject to an average 2.70 mm d⁻¹. This value is only 0.07 mm d⁻¹ less than the 2002–2014 climatological mean, indicating that 2014 is an average year for precipitation and not anomalously dry. Individual GCM precipitation rates taken from their Hol₁₄ simulations capture the observed spatial pattern and magnitude well; only a minute difference exists of 0.01 mm d⁻¹ between the multi-model ensemble and GPCP data across the PVP region. Higher precipitation rates are found within the wider domain, notably south-west of Iceland. Observed precipitation in October is slightly higher yet also average in comparison to the climatology and is well captured by the GCMs (see Fig. C4). As evidence exists of appreciable precipitation in both the GPCP data and GCMs, there should be scope for the added aerosol from Holuhraun to influence precipitation processes – and subsequently bring forth changes related to the second indirect effect – within both the real-world and modelled cloud systems.

Figure 5 Monthly mean anomalies in all-sky liquid water path (LWP) for September 2014 from (a) MODIS instruments on board Aqua and Terra satellites, (b) the multi-model ensemble, and (c–j) individual models. Anomalies depicted are the total effect, so they include both aerosol and meteorological components (i.e. Hol₁₄ – NoHol_clim). The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Stippling highlights grid cells with null hypothesis rejections based on applying the False Discovery Method (FDR) at a 10 % control level (see main text).

We explore the spatial pattern of a possible second indirect effect using LWP – a common proxy for precipitation suppression. The total perturbation in all-sky LWP observed by MODIS during September 2014 is shown in Fig. 5a. As before, stippling indicates grid elements with rejected null hypotheses after applying the FDR method at 10 %. Within the PVP region, our statistical testing identifies a small area where the observed change in LWP is unlikely due to random variation, suggesting the clouds here are retaining more water as a result of the aerosol perturbation introduced by Holuhraun (i.e. precipitation suppression). However, for the vast majority of the domain, any aerosol influence on the observed changes in LWP – both inside and outside the PVP region – are relatively minor compared to meteorological variability. Here, specifically, the positive correlation between LWP and precipitation is shown – areas with higher (lower) cloud liquid water content often support more (less) cloud droplet formation and so increased (decreased) precipitation. Modelled total anomalies in LWP are depicted in Fig. 5b–j. Whilst the GCMs capture the spatial patterns well, there is clear variation in the magnitude of the anomalies, with CAM5.3-Oslo and the ECHAM6 variants showing strong LWP responses relative to MODIS in and near the PVP region. Nonetheless, the response of the multi-model ensemble differs only slightly to the observed (Δ LWP =4.49 gm⁻²). However, unlike the ACI first indirect effect, the multi-model ensemble response in October (see Fig. C5) is only ∼ 50 % that of MODIS, suggesting the September response is likely coincidental. Aside from this, similar observed and modelled behaviour is found for October. Equivalent spatial figures for CF can be found in the Appendices (see Figs. B2 and C6).

Figure 6 Disentanglement of the aerosol and meteorological effects on (a) all-sky liquid water path (LWP) and (b) total cloud fraction (CF) within the predominantly volcanically polluted (PVP) region for September 2014. Total, aerosol-only, and meteorology-only effects are depicted by green – no pattern, red – minor diagonal, and blue – major diagonal box plots respectively. Box plots extend to the 25th–75th percentiles with outer whiskers at 5th–95th. Black squares depict means. Green bounding and dashed lines extend the observed total effects across rows for visual comparison with the model responses. Climatological baselines are given in brackets.

Figure 7 Monthly mean anomalies in the rate of cloud droplet autoconversion for September 2014 from (a–e) select individual models and (f) the multi-model ensemble. Model responses depict aerosol-only anomalies (i.e. Hol₁₄ – NoHol₁₄). The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Note that the aerosol-only effect on cloud droplet autoconversion cannot be calculated for HadGEM3 and ECHAM6-HAM from the output provided to this experiment, whilst CNRM-ESM2-1 is not considered here (see main text).

A breakdown of the aerosol and meteorological components of the modelled LWP and CF responses alongside MODIS observations for the PVP regions of September and October 2014 is given in Figs. 6 and C7 respectively, with summary values provided in Tables D1 and D2. Focusing first on the LWP decomposition, the GCMs clearly diverge in the total effect caused by the eruption, with a roughly equal number of models over- and underestimating the impact noted by MODIS. This discrepancy is due mainly to the variation in the simulated aerosol effects, rather than the meteorological effects. For example, in September the mean meteorological component across the individual GCMs varies by 21.68 gm⁻², whilst for aerosol this spread is 38.95 gm⁻² – almost double. Across the 2 months, the two UKESM1 variants and HadGEM3 simulate a moderate aerosol response (∼ 4–8 gm⁻²), whereas a considerably stronger aerosol response (∼ 20–40 gm⁻²) is simulated in CAM5.3-Oslo and the three ECHAM6 variants. As expected, the aerosol response for CNRM-ESM2-1 is negligible due to the absence of an aerosol–precipitation mechanism within this model (see Sect. 2.1).

To investigate the moderate and strong aerosol responses in LWP, we explore the aerosol-only effect on the monthly mean rate of cloud droplet autoconversion for September and October in Figs. 7 and C8 respectively. For the models with autoconversion data available, the volcanic aerosol acts to decrease the rate of autoconversion. This is likely due to the simulated increase in N_d acting to inhibit cloud droplet growth via the $N_{\mathrm{d}}^{-\mathrm{1.79}}$ inverse dependence in the KK2000 parameterisation that this subset of models all use. Harder to interpret is the variation in the magnitude of the autoconversion decrease. One contributing factor may be how the GCMs represent sub-grid variability in cloud liquid water, which influences the $q_{\mathrm{cl}}^{\mathrm{2.47}}$ term in KK2000 that can counteract the suppression driven by N_d. For instance, implicit descriptions, such as those in the ECHAM6 models, may underestimate this non-linear offset in autoconversion if its scaling term does not accurately reflect the local areas of high cloud liquid water, a feat particularly challenging given the non-uniform conditions of volcanic plumes, whereas explicit descriptions involving probability density functions, such as the log-normal and gamma distributions used in the UKESM1 variants and CAM5.3-Oslo respectively, are potentially able to better capture this non-linear impact on autoconversion from these local high values. Furthermore, it is worth noting that none of the GCMs explicitly account for sub-grid variability in N_d which, whilst unlikely to contribute to the variation in responses, may cause an underestimation of the suppression via $N_{\mathrm{d}}^{-\mathrm{1.79}}$ if sufficiently high N_d localities exist.

Regarding the total CF, no substantial overall change is observed by MODIS within the PVP regions of either month – a finding emulated by the models. The aerosol–meteorology decomposition made possible by the GCMs suggests that the meteorological variability dominates the total effect on CF at the monthly scale, making any conclusion on the aerosol related impact challenging. Nevertheless, a minor increase in total CF due to the added aerosol is simulated by all models except for CAM5.3-Oslo.

Figure 8 Monthly mean anomalies in top-of-atmosphere upwelling shortwave radiation (rsut) for September 2014 from (a) CERES-EBAF, (b) the multi-model ensemble, and (c–j) individual models. Anomalies depicted are the total effect, so they include both aerosol and meteorological components (i.e. Hol₁₄ – NoHol_clim). Here radiative fluxes are positive downward. The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Stippling highlights grid cells with null hypothesis rejections based on applying the False Discovery Method (FDR) at a 10 % control level (see main text).

5 Top-of-atmosphere radiative response

Here we examine the influence of the volcanic aerosol introduced by Holuhraun on the Earth's energy budget. The total effect on ToA upwelling SW radiation (rsut) for September 2014 given by CERES-EBAF is illustrated in Fig. 8a, where increased upward radiative flux is treated as a negative change. Once again, local null hypothesis tests subject to the FDR method at 10 % were conducted. There is mainly an observed increase in rsut across the North Atlantic Ocean following the eruption, with the few areas subject to opposing behaviour largely near land masses in the south (e.g. Celtic Sea, Irish Sea, Baltic Sea, Labrador Sea). Some of the same meteorological features as those depicted in the LWP response are present, suggesting again that meteorological variability is clouding any possible observable aerosol signal on rsut. The associated modelled total effects are shown in Fig. 8b–j. The observed spatial pattern is captured well by the models, yet the magnitude varies with most GCMs overestimating the increase in rsut. This discrepancy is most apparent between 45–55° N. For October, an improvement in the model performance is noted, with only a difference of 0.09 Wm⁻² between CERES-EBAF and the multi-model ensemble (see Fig. C9). Corresponding spatial figures for ToA upwelling LW radiation (rlut) can be found in the Appendix (see Figs. B3 and C10).

Figure 9 Disentanglement of the aerosol and meteorological effects on top-of-atmosphere upwelling (a) shortwave (rsut) and (b) longwave (rlut) radiation within the predominantly volcanically polluted (PVP) region for September 2014. Total perturbations, and their aerosol and meteorological components, are depicted by green – no pattern, red – minor diagonal, and blue – major diagonal box plots respectively. Box plots extend to the 25th–75th percentiles with outer whiskers at 5th–95th. Black squares depict means. Green bounding and dashed lines extend the observed total effects across rows to aid visual comparison with the model responses. Increased upward radiative flux is treated as a negative change. Climatological baselines are given in brackets.

The disentanglement of the aerosol signal from the meteorological variability for rsut and rlut for the PVP regions of September and October 2014 is shown in Figs. 9 and C11 respectively, with summarising values provided in Tables D1 and D2. All models simulate an overall increase in rsut in the PVP regions, as is observed by CERES-EBAF, yet most models overestimate this change, particularly in September, with notable examples including ECHAM6-SALSA and CNRM-ESM2-1 that respectively simulate responses 106 % and 75 % stronger than observed. The modelled decomposition of the overall increase in rsut shows that the newly introduced aerosol is the predominant cause – likely due to increasing cloud albedo – rather than the meteorological component which often acts to oppose this volcanic influence. In comparison, the aerosol effect on LW radiation leaving the Earth system is minor and dominated by meteorological variability. Nevertheless, for all except UKESM1-BLN, this minor effect is to decrease rlut. This is possibly due to changes in the aerosol direct effect, specifically scattering due to the non-absorbing nature of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, or that low clouds have become thicker due to increased LWP and subsequently trap more upwelling LW radiation. Further analysis with additional diagnostics is needed to confirm this (e.g. using Ghan, 2013, methodology). Overall, for both observed and modelled responses, increases in rsut outweigh decreases in rlut, suggesting the Holuhraun eruption prompted a net cooling effect on the Earth's energy budget.

Table 3 Radiative forcing (RF) estimates from the Holuhraun eruption across the predominantly volcanically polluted (PVP) regions and globe. Global RF estimates are scaled from RF estimates of the entire Northern Hemisphere above 50° N to exclude noise (see main text).

Furthermore, we estimate the strength of this cooling effect using the GCMs. As incoming solar radiation is the same across the Hol₁₄ and NoHol₁₄ simulations, the net change in rsut and rlut between them (i.e. the aerosol-only effect) approximates the RF due to Holuhraun. The local RFs for the September and October PVP regions are listed in Table 3. The model responses vary by ∼ 2 Wm⁻² for both months/PVP regions, with the ECHAM6 variants and UKESM1-BLN generally simulating the strongest and weakest forcings respectively. Overall the RF is stronger in September when the SO₂ emissions are at their highest and solar insolation is greater. In addition, we determine global RF estimates to allow comparison of the influence Holuhraun had on the Earth's energy system versus other events. Global values are scaled up from RF estimates of the entire Northern Hemisphere above 50° N and ignore RF contributions outside this area, a choice made to reduce the influence of noise, namely from equatorial regions, as changes in ToA fluxes there are unlikely due to Holuhraun given the spatial evolution of the plume evidenced in Part 1 (Jordan et al., 2024). Averaged across September and October, we find that all models display a global negative forcing in response to the additional aerosol, with our multi-model ensemble estimating a value of −0.11 ± 0.04 Wm⁻² (± 1σ of the individual model RFs). The ECHAM6-SALSA global RF is nearly twice that of the ensemble, with the additional forcing potentially due to its consistently strong LWP response across September and October relative to the other models. On the other hand, CNRM-ESM2-1 shows the smallest global RF, roughly a third of the ensemble mean, and is likely due to the exclusion of precipitation suppression induced ACI indirect effects within this model. We assume that the difference across the wider Holuhraun domain of ⪅20 % between the September–October and annual solar insolation means is small enough to allow us to extrapolate our September–October global RFs to provide proxy annual global RFs. Our multi-model ensemble suggests that, averaged over a year, the added aerosol from Holuhraun caused a forcing of −0.018 ± 0.007 Wm⁻². Given that Holuhraun released 3.9 Tg of SO₂ in our simulations over this period (see Table 2), we estimate a global mean annual RF efficiency for the eruption of −0.005 ± 0.002 Wm⁻² per Tg of SO₂. In reality, Holuhraun volcanic activity did not cease in October and continued until February, albeit at a lesser extent, and released an estimated total 9.6–11.8 Tg of SO₂; hence our annual forcing estimates should be considered minimums.

6 Summary and conclusion

The continuous degassing of the 2014–2015 Holuhraun eruption into the lower troposphere resulted in a persistent source of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ pollution across the northern North Atlantic Ocean, providing an opportunistic experiment to assess the representation of ACIs in state-of-the-art GCMs. Here we have presented Part 2 of an AeroCom Phase III inter-model comparison two-part study designed to leverage this opportunity and build on previous works utilising GCMs (Gettelman et al., 2015; Jordan et al., 2024; Malavelle et al., 2017). A simple theoretical framework designed to separate the aerosol and meteorological effects on cloud properties is applied to eight GCMs across regions identified with minimal non-Holuhraun sources of aerosol pollution during September and October 2014. By comparing the resulting decomposition of the cloud responses to observations from a range of remote sensing instruments, we review the ACI model representations and highlight those that deviate away from the observed behaviour.

Regarding the ACI first indirect effect (i.e. Twomey effect), MODIS observations suggest notable increases and decreases in cloud top N_d and r_e respectively across the PVP regions of September and October 2014 when compared to their respective long-term averages. All models correctly capture the direction of these observed changes in cloud top N_d and r_e, yet the magnitude of their responses varies. Applying our analysis framework shows that the differences in cloud top N_d and r_e relative to their climatological values are almost entirely due to the aerosol added by the eruption rather than interannual variability driven by meteorological influence, a finding in agreement with previous studies (e.g.  Chen et al., 2022; Malavelle et al., 2017). Despite the differences in the strength of the aerosol induced model responses – which are largely explainable by configuration choices – the multi-model ensemble representation of the ACI first indirect effect agrees well with MODIS observations, increasing our confidence in using ensemble based methods to explore this mechanism elsewhere.

For the ACI second indirect effect (i.e. rapid adjustments), we show that in both the real-world and modelled cloud systems precipitation is close to the GPCP climatological average and not anomalously dry following the eruption, meaning aerosol-invoked precipitation suppression is possible. We use all-sky LWP and total CF as our proxies to assess whether a delay of precipitation formation is causing macrophysical changes in the clouds. Unlike the microphysical changes in N_d and r_e, an aerosol response in LWP and CF is far harder to discern amongst the meteorological variability, a complication previously reported (Malavelle et al., 2017; McCoy and Hartmann, 2015; Peace et al., 2024). Nevertheless, our disentanglement method allows us to isolate the aerosol signal within the PVP regions. The GCMs – excluding CNRM-ESM2-1 due to its autoconversion scheme – show a positive LWP response to the added aerosol suggesting precipitation suppression is simulated. However, a clear divergence in magnitude is evident, which may stem from differences in implementing the KK2000 autoconversion scheme within the models, particularly in the treatment of sub-grid variability. Given the heterogenous nature of the high aerosol load introduced by Holuhraun, inaccurate representations of sub-grid variability are likely to have a greater impact than under more typical conditions. Moreover, aside from CAM5.3-Oslo, all models simulate a positive volcanic influence on CF, yet the magnitude is minor compared to meteorological variability. In comparison, Chen et al. (2022), via machine learning techniques, isolate the aerosol signal within MODIS observations and find a far larger increase in CF. If this is the case, then the model CF responses presented here are underestimated and further work to ascertain why is needed.

We show that the volcanic influence on ToA radiation within the PVP regions is predominantly on SW radiation rather than LW, with the net effect being an increase in radiation leaving the Earth system. Our multi-model ensemble mean estimates that this cooling has a global radiative forcing of −0.11 ± 0.04 Wm⁻² averaged over September and October, revising previous estimates made using individual GCMs (Gettelman et al., 2015; Malavelle et al., 2017). Such forcing is comparable to that caused by weak–moderate explosive eruptions (e.g. Kasatochi, Narbo, Sarychev Peak, Raikoke), with SO₂ emissions an order of magnitude less than Holuhraun yet 10–15 km higher in the atmosphere (Schallock et al., 2023). For Holuhraun, we estimate a global mean annual RF efficiency of −0.005 ± 0.002 Wm⁻² per Tg of SO₂. For comparison, 2014 global anthropogenic SO₂ emissions had an approximate RF efficiency of −0.010 ± 0.004 Wm⁻² per Tg of SO₂ (Aas et al., 2019; Szopa et al., 2021; Thornhill et al., 2021), whereas a recent reduction in shipping SO₂ emissions incited by 2020 regulations yielded a RF efficiency of −0.014 ± 0.002 Wm⁻² per Tg of SO₂ (Jordan and Henry, 2024). Whilst our Holuhraun estimate and these values are in fair agreement, the differences would likely reduce if Holuhraun had occurred during spring–summer and/or in a cloud regime more susceptible to aerosol changes as both would act to increase the cooling effect – a notion shared by other studies (Malavelle et al., 2017; Zoëga et al., 2025). Similarly, as the consensus of the GCMs is that the net effect of the meteorological impact acts to oppose the volcanic influence, a greater cooling effect would also occur if Holuhraun had erupted under more favourable meteorological conditions.

Despite our best efforts, our study is subject to limitations. Observations are subject to the general constraints of satellite remote sensing at high latitudes, whereas modelling caveats include varied cloud system susceptibility due to differing background aerosol concentrations across the models and non-uniformity in the modelled aerosol perturbations/plume representations (e.g. Jordan et al., 2024). Nevertheless, our two-part study of the Holuhraun eruption has used novel techniques to explore GCM representations of ACIs during a high pollution event, confirming their ability to capture the first indirect effect well, whilst highlighting discrepancies in their second indirect effect responses and noting the refinement of their autoconversion schemes as a potential route to improvement.

Appendix A: KK2000 autoconversion parameterisation variations

The majority of GCMs in this study base their autoconversion parameterisation on the scheme by Khairoutdinov and Kogan (2000), hereafter KK2000. KK2000 calculates the rate of autoconversion of cloud water to rain $\left(\frac{\partial q_{\mathrm{cl}}}{\partial t}{|}_{\text{auto}}\right)$ as a function of cloud liquid water content, q_cl, and cloud droplet number concentration, N_d. The general form of KK2000 is given by

$$\begin{array}{}\text{(A1)}& {\left.{\displaystyle \frac{\partial q_{\mathrm{cl}}}{\partial t}}\right|}_{\text{auto}}=\mathrm{1350}{q}_{\mathrm{cl}}^{\mathrm{2.47}}{N}_{\mathrm{d}}^{-\mathrm{1.79}}\end{array}$$

and captures the non-linear dependence of autoconversion on q_cl and N_d. An important characteristic of KK2000 is that it inhibits the rate of autoconversion following an increase in the number of (smaller) cloud droplets. This enables KK2000 to represent the precipitation suppression mechanism of the ACI second indirect effect (see Sect. 1).

KK2000 was intended for spatial resolutions of 𝒪(10 m) – scales able to explicitly resolve cloud-scale processes. To implement KK2000 in large-scale models, such as the GCMs in this study, an additional term is introduced to add a representation of the sub-grid variability of q_cl and/or N_d – an important factor as the non-linear nature of autoconversion means using grid cell average values will likely add substantial biases (e.g. Larson et al., 2001). The simplest way to account for this sub-grid variability is by applying a scaling term, −γ_l, to Eq. (A1) as follows:

$$\begin{array}{}\text{(A2)}& {\left.{\displaystyle \frac{\partial q_{\mathrm{cl}}}{\partial t}}\right|}_{\text{auto}}=-{\mathit{\gamma }}_l\mathrm{1350}{q}_{\mathrm{cl}}^{\mathrm{2.47}}{N}_{\mathrm{d}}^{-\mathrm{1.79}}.\end{array}$$

The tuning of γ_l adjusts the efficiency of the autoconversion rate such that increasing γ_l quickens the rate at which cloud droplets grow by collision–coalescence. Here γ_l implicitly accounts for sub-grid variability. ECHAM6-HAM, ECHAM6-HAM-P3, and ECHAM6-SALSA adopt this form of KK2000 using γ_l values of 10.7, 2.7, and 2.8 respectively. For further information, including the tuning method, see Lohmann and Neubauer (2018) and Dietlicher et al. (2019).

A more detailed approach to characterise sub-grid variability is to use a probability density function (PDF). PDFs can explicitly represent the distribution of q_cl and/or N_d across a grid cell, capturing those important local features. In CAM5.3-Oslo, a gamma distribution, Γ, is used to describe the sub-grid variability of q_cl, altering Eq. (A1) as follows:

$$\begin{array}{}\text{(A3)}& {\left.{\displaystyle \frac{\partial q_{\mathrm{cl}}}{\partial t}}\right|}_{\text{auto}}=-\mathrm{CF}{\displaystyle \frac{\mathrm{\Gamma }(\mathit{\nu }+\mathrm{2.47})}{\mathrm{\Gamma }\left(\mathit{\nu }\right){\mathit{\nu }}^{\mathrm{2.47}}}}\mathrm{1350}{q}_{\mathrm{cl}}^{\mathrm{2.47}}{N}_{\mathrm{d}}^{-\mathrm{1.79}},\end{array}$$

where CF is cloud fraction and ν is a gamma function shape parameter defined by the relative variance of q_cl within the grid cell. Further details are provided in Morrison and Gettelman (2008) and Gettelman and Morrison (2015).

Moreover, sometimes a log-normal distribution is used to describe the PDF of the underlying sub-grid variability. This method is done in the configurations of UKESM1 and HadGEM3 assessed in this study. By assuming q_cl follows a log-normal distribution at the sub-grid scale, along with added correction factors, these models implement a refined version of Eq. (A1) given as follows:

$$\begin{array}{}\text{(A4)}& {\left.{\displaystyle \frac{\partial q_{\mathrm{cl}}}{\partial t}}\right|}_{\text{auto}}=E\left(f_{\mathrm{cl}}\right)\mathrm{1350}{q}_{\mathrm{cl}}^{\mathrm{2.47}}{N}_{\mathrm{d}}^{-\mathrm{1.79}},\end{array}$$

where the sub-grid variability term, E(f_cl), and the parameter f_cl are given by

$$\begin{array}{}\text{(A5)}& {\displaystyle }E\left(f_{\text{cl}}\right)=(\mathrm{1}+f_{\mathrm{cl}}^{\mathrm{2}}{)}^{-\mathrm{2.47}/\mathrm{2}}(\mathrm{1}+f_{\mathrm{cl}}^{\mathrm{2}}{)}^{{\mathrm{2.47}}^{\mathrm{2}}/\mathrm{2}}\text{(A6)}& {\displaystyle }{f}_{\mathrm{cl}}=\left\{\begin{array}{l}(\mathrm{0.45}-\mathrm{0.25}\mathrm{CF})\sqrt{(x\mathrm{CF}{)}^{\mathrm{2}/\mathrm{3}}}{\left((\mathrm{0.06}x\mathrm{CF}{)}^{\mathrm{1.5}}+\mathrm{1}\right)}^{-\mathrm{0.17}}\\ \text{if }\mathrm{CF}<\mathrm{1}\\ \mathrm{0.11}\sqrt{(x\mathrm{CF}{)}^{\mathrm{2}/\mathrm{3}}}{\left((\mathrm{0.06}x\mathrm{CF}{)}^{\mathrm{1.5}}+\mathrm{1}\right)}^{-\mathrm{0.17}}\\ \text{if }\mathrm{CF}=\mathrm{1}.\end{array}\right.\end{array}$$

Here x is the grid cell size. For further information, including the correction factors applied, see Boutle et al. (2014). Finally, it is worth noting that none of the models assessed in this study explicitly account for the sub-grid variability of N_d, only q_cl.

Appendix B: Additional September 2014 figures

Figure B1 Monthly mean anomalies in cloud droplet number concentration (N_d) at cloud top for September 2014 from (a) MODIS instruments on board Aqua and Terra satellites, (b) the multi-model ensemble, and (c–h) individual models. Anomalies depicted are the total effect, so they include both aerosol and meteorological components (i.e. Hol₁₄ – NoHol_clim). The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Stippling highlights grid cells with null hypothesis rejections based on applying the False Discovery Method (FDR) at a 10 % control level (see main text). Hatched areas indicate missing data. Note that the total effect on N_d at cloud top cannot be calculated for ECHAM6-HAM-P3 and ECHAM6-SALSA from the output provided to this experiment.

Figure B2 Monthly mean anomalies in total cloud fraction (CF) for September 2014 from (a) MODIS instruments on board Aqua and Terra satellites, (b) the multi-model ensemble, and (c–j) individual models. Anomalies depicted are the total effect, so they include both aerosol and meteorological components (i.e. Hol₁₄ – NoHol_clim). The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Stippling highlights grid cells with null hypothesis rejections based on applying the False Discovery Method (FDR) at a 10 % control level (see main text).

Figure B3 Monthly mean anomalies in top-of-atmosphere upwelling longwave radiation (rlut) for September 2014 from (a) CERES-EBAF, (b) the multi-model ensemble, and (c–j) individual models. Anomalies depicted are the total effect, so they include both aerosol and meteorological components (i.e. Hol₁₄ – NoHol_clim). Here radiative fluxes are positive downward. The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Stippling highlights grid cells with null hypothesis rejections based on applying the False Discovery Method (FDR) at a 10 % control level (see main text).

Appendix C: October 2014 figures

Figure C1 Monthly mean anomalies in cloud droplet number concentration (N_d) at cloud top for October 2014 from (a) MODIS instruments on board Aqua and Terra satellites, (b) the multi-model ensemble, and (c–h) individual models. Anomalies depicted are the total effect, so they include both aerosol and meteorological components (i.e. Hol₁₄ – NoHol_clim). The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Stippling highlights grid cells with null hypothesis rejections based on applying the False Discovery Method (FDR) at a 10 % control level (see main text). Hatched areas indicate missing data. Note that the total effect on N_d at cloud top cannot be calculated for ECHAM6-HAM-P3 and ECHAM6-SALSA from the output provided to this experiment.

Figure C2 Monthly mean anomalies in cloud droplet effective radius (r_e) at cloud top for October 2014 from (a) MODIS instruments on board Aqua and Terra satellites, (b) the multi-model ensemble, and (c–h) individual models. Anomalies depicted are the total effect, so they include both aerosol and meteorological components (i.e. Hol₁₄ – NoHol_clim). The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Stippling highlights grid cells with null hypothesis rejections based on applying the False Discovery Method (FDR) at a 10 % control level (see main text). Hatched areas indicate missing data. Note that the total effect on r_e at cloud top cannot be calculated for ECHAM6-HAM-P3 and ECHAM6-SALSA from the output provided to this experiment.

Figure C3 Disentanglement of the aerosol and meteorological effects on (a) cloud droplet number concentration (N_d) and (b) cloud droplet effective radius (r_e) at cloud top within the predominantly volcanically polluted (PVP) region for October 2014. Total perturbations, and their aerosol and meteorological components, are depicted by green – no pattern, red – minor diagonal, and blue – major diagonal box plots respectively. Box plots extend to the 25th–75th percentiles with outer whiskers at 5th–95th. Black squares depict means. Green bounding and dashed lines visualise the observed total effects across the model responses. Climatological baselines are given in brackets. Note that the aerosol effect is only available for ECHAM6-HAM-P3 and ECHAM6-SALSA (see main text).

Figure C4 Monthly mean surface precipitation rates for October 2014 from (a) the Global Precipitation Climatology Project (GPCP), (b) the multi-model ensemble, and (d–k) individual models, as well as (c) the climatological October mean (2002–2014) derived from GPCP. The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Modelled precipitation rates are for the simulations including Holuhraun emissions (i.e. Hol₁₄).

Figure C5 Monthly mean anomalies in all-sky liquid water path (LWP) for October 2014 from (a) MODIS instruments on board Aqua and Terra satellites, (b) the multi-model ensemble, and (c–j) individual models. Anomalies depicted are the total effect, so they include both aerosol and meteorological components (i.e. Hol₁₄ – NoHol_clim). The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Stippling highlights grid cells with null hypothesis rejections based on applying the False Discovery Method (FDR) at a 10 % control level (see main text). Hatched areas indicate missing data.

Figure C6 Monthly mean anomalies in total cloud fraction (CF) for October 2014 from (a) MODIS instruments on board Aqua and Terra satellites, (b) the multi-model ensemble, and (c–j) individual models. Anomalies depicted are the total effect, so they include both aerosol and meteorological components (i.e. Hol₁₄ – NoHol_clim). The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Stippling highlights grid cells with null hypothesis rejections based on applying the False Discovery Method (FDR) at a 10 % control level (see main text).

Figure C7 Disentanglement of the aerosol and meteorological effects on (a) all-sky liquid water path (LWP) and (b) total cloud fraction (CF) within the predominantly volcanically polluted (PVP) region for October 2014. Total, aerosol-only, and meteorology-only effects are depicted by green – no pattern, red – minor diagonal, and blue – major diagonal box plots respectively. Box plots extend to the 25th–75th percentiles with outer whiskers at 5th–95th. Black squares depict means. Green bounding and dashed lines extend the observed total effects across rows for visual comparison with the model responses. Climatological baselines are given in brackets.

Figure C8 Monthly mean anomalies in the rate of cloud droplet autoconversion for October 2014 from (a–e) select individual models and (f) the multi-model ensemble. Model responses depict aerosol-only anomalies (i.e. Hol₁₄ – NoHol₁₄). The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Note that the aerosol-only effect on cloud droplet autoconversion cannot be calculated for HadGEM3 and ECHAM6-HAM from the output provided to this experiment, whilst CNRM-ESM2-1 is not considered here (see main text).

Figure C9 Monthly mean anomalies in top-of-atmosphere upwelling shortwave radiation (rsut) for October 2014 from (a) CERES-EBAF, (b) the multi-model ensemble, and (c–j) individual models. Anomalies depicted are the total effect, so they include both aerosol and meteorological components (i.e. Hol₁₄ – NoHol_clim). Here radiative fluxes are positive downward. The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Stippling highlights grid cells with null hypothesis rejections based on applying the False Discovery Method (FDR) at a 10 % control level (see main text).

Figure C10 Monthly mean anomalies in top-of-atmosphere upwelling longwave radiation (rlut) for October 2014 from (a) CERES-EBAF, (b) the multi-model ensemble, and (c–j) individual models. Anomalies depicted are the total effect, so they include both aerosol and meteorological components (i.e. Hol₁₄ – NoHol_clim). Here radiative fluxes are positive downward. The predominantly volcanically polluted (PVP) region is outlined by a dashed line, with its spatial mean listed above. Stippling highlights grid cells with null hypothesis rejections based on applying the False Discovery Method (FDR) at a 10 % control level (see main text).

Figure C11 Disentanglement of the aerosol and meteorological effects on top-of-atmosphere upwelling (a) shortwave (rsut) and (b) longwave (rlut) radiation within the predominantly volcanically polluted (PVP) region for October 2014. Total perturbations, and their aerosol and meteorological components, are depicted by green – no pattern, red – minor diagonal, and blue – major diagonal box plots respectively. Box plots extend to the 25th–75th percentiles with outer whiskers at 5th–95th. Black squares depict means. Green bounding and dashed lines extend the observed total effects across rows to aid visual comparison with the model responses. Increased upward radiative flux is treated as a negative change. Climatological baselines are given in brackets.

Appendix D: Aerosol–meteorology disentanglement summary tables

Table D1 September 2014 aerosol–meteorology disentanglement. Shown are the predominantly volcanically polluted (PVP) regional means of the total, aerosol-only, and meteorology-only effects, as well as a climatological baseline, for cloud top cloud droplet number concentration (N_d), cloud top cloud droplet effective radius (r_e), all-sky liquid water path (LWP), total cloud fraction (CF), top-of-atmosphere upwelling shortwave radiation (rsut), and top-of-atmosphere upwelling longwave radiation (rlut). Note that for ECHAM6-HAM-P3 and ECHAM6-HAM-SALSA only the aerosol responses in cloud top N_d and r_e are available (see main text).

^* MODIS observations used for cloud top N_d, cloud top r_e, all-sky LWP, and CF. CERES-EBAF observations used for rsut and rlut.

Table D2 October 2014 aerosol–meteorology disentanglement. Shown are the predominantly volcanically polluted (PVP) regional means of the total, aerosol-only, and meteorology-only effects, as well as a climatological baseline, for cloud top cloud droplet number concentration (N_d), cloud top cloud droplet effective radius (r_e), all-sky liquid water path (LWP), total cloud fraction (CF), top-of-atmosphere upwelling shortwave radiation (rsut), and top-of-atmosphere upwelling longwave radiation (rlut). Note that for ECHAM6-HAM-P3 and ECHAM6-HAM-SALSA only the aerosol responses in cloud top N_d and r_e are available (see main text).

^* MODIS observations used for cloud top N_d, cloud top r_e, all-sky LWP, and CF. CERES-EBAF observations used for rsut and rlut.