Process-level contributions to uncertainty in aerosol effective radiative forcing: a perturbed parameter ensemble with ICON–HAM

Changes in aerosols since the preindustrial era have altered the top-of-the-atmosphere radiation balance by scattering and absorbing solar radiation (aerosol–radiation interactions, ARI) and indirectly interacting with clouds (aerosol–cloud interactions, ACI), known as aerosol effective radiative forcing (ERF_aer). ERF_aer persistently remains one of the most uncertain components in climate projections, due to imperfect representations of aerosol and cloud processes in climate models. Here, we construct a perturbed parameter ensemble (PPE) with the aerosol–climate model ICON2.6.4–A–HAM2.3 (hereafter ICON–HAM) to quantify key sources of ERF_aer uncertainty. We perturb 42 aerosol and cloud parameters over 383 PPE members. Parametric uncertainties in aerosol and cloud processes yield a 90 % credible range of −1.42 to −0.65 W m⁻² for ERF_aer (mean: −1.04 W m⁻²) for the period 2024–2025. The parameters related to emissions (anthropogenic sulfur dioxide, natural dimethyl sulfide, and emitted particle size) dominate ACI uncertainty and hence ERF_aer uncertainty (80 %), while absorption-related parameters (anthropogenic black carbon emissions and aerosol refractive indices) drive ARI uncertainty (60 %). Cloud parameters account for 13 % of ERF_aer uncertainty, mainly via convection and entrainment processes. The sensitivity analysis of model diagnostics to parameters reveals that many present-day aerosol and cloud observables share dominant causes of uncertainty with ACI and ARI forcing, highlighting the potential for constraining ERF_aer using existing space- and ground-based measurements. Notably, model biases against SPEXone and MODIS observations coincide spatially with parametric uncertainties, suggesting that much of these biases may be mitigated through appropriate constraint with observations, while the remainder requires structural model developments in combination with improved observations.

1 Introduction

Earth's climate system is modulated by atmospheric aerosols, which perturb the radiative energy balance both directly through scattering and absorption of radiation (aerosol–radiation interactions, ARI) and indirectly by acting as cloud condensation nuclei (CCN) and ice-nucleating particle (INP) that modify cloud microphysical and radiative properties (aerosol–cloud interactions, ACI). Since the industrial era, human activities have significantly increased aerosols and altered their composition, leading to an effective radiative forcing from aerosols (ERF_aer) that arises from both ARI (ERF_ari) and ACI (ERF_aci) (Bellouin et al., 2020; Forster et al., 2021). The latest assessment by the Intergovernmental Panel on Climate Change (IPCC) provided an ERF_aer estimate of −1.1 W m⁻² over 1750–2019, with a 90 % credible range of −1.7 to −0.4 W m⁻² (Forster et al., 2021). This aerosol-induced cooling partly offsets the warming driven by anthropogenic greenhouse gases (Szopa et al., 2021).

Despite sustained advances in observations, modeling, and process understanding over recent decades, the estimate of ERF_aer remains highly uncertain, with particularly large uncertainty in ERF_aci (Seinfeld et al., 2016; Kahn et al., 2023; Rosenfeld et al., 2023; Im et al., 2026). Since IPCC AR5, the observational estimation of ERF_aer has substantially progressed along several axes (Boucher et al., 2013; Forster et al., 2021), leading to improved estimates with reduced uncertainty (Hasekamp et al., 2019b; Bellouin et al., 2020; Jia et al., 2026). In contrast, progress on the modeling side has been more limited: from phase 5 to phase 6 of the Coupled Model Intercomparison Project (CMIP), the spread in simulated ERF_aer has narrowed only marginally (Forster et al., 2021). This implies that advances in our fundamental understanding have not yet been translated meaningfully into reduced model uncertainty and reliable climate projections, despite increases in model complexity and resolution (Carslaw et al., 2026).

Multi-model ensembles (MMEs) have often been used to assess model spread and thereby prediction reliability (i.e., the consistency and robustness of model results across different scenarios and application contexts) in model intercomparison projects (MIP) such as AeroCom (Myhre et al., 2013a; Gliß et al., 2021) and CMIP (Durack et al., 2025). However, the model spread in an MME reflects a complex mixture of structural uncertainty (differences in model equations and process representations), parametric uncertainty (uncertainty in parameter values within a fixed model structure), and internal variability (natural fluctuations of the climate system) (Knutti et al., 2010). This makes it difficult to pinpoint the processes driving model uncertainty, prioritise future model development, and ultimately narrow the uncertainty in climate predictions and projections (Carslaw et al., 2026).

Perturbed parameter ensembles (PPEs) provide a complementary approach that disentangles structural and parametric origins of model-observation discrepancies and attributes model uncertainty to specific physical processes (Allen, 1999; Carslaw et al., 2013; Regayre et al., 2026). A PPE comprises a suite of simulations from a single model in which multiple physical parameters are simultaneously varied to sample parametric uncertainty (Murphy et al., 2004). PPEs have been widely applied to atmospheric science, and more recently adopted to aerosol–climate modeling. Many modeling centres have developed PPE infrastructures, including CAM5 (Zhao et al., 2013), CAM6 (Eidhammer et al., 2024; Gettelman et al., 2024), E3SM (Qian et al., 2018; Nugent et al., 2026), ECHAM–HAM (Bhatti et al., 2026), GISS ModelE (Elsaesser et al., 2025), HadGEM (Regayre et al., 2018; Yoshioka et al., 2019), and UKESM1 (Regayre et al., 2023), through which important sources of ERF_aer uncertainty have been identified and constrained in individual models (Watson-Parris et al., 2020; Johnson et al., 2020; Song et al., 2025; Regayre et al., 2026). Most PPEs typically perturb a relatively limited number of parameters (∼ 25) or focus primarily on either aerosol or cloud processes, potentially limiting their ability to capture complex nonlinear interactions among processes. Another important aspect of PPE design is the simulations-to-parameters ratio (r_sp), which typically ranges from about 5 to 10 and governs the trade-off between emulator accuracy and computational cost. Yet, optimal values of r_sp that balance this trade-off remains uncertain. This poses a practical challenge for future PPE MIPs, given their substantial demands on computational resources and data storage; for example, a single model PPE can generate ∼ 47 TB of data (Eidhammer et al., 2024).

Building on the prior knowledge of parameter importance, we construct a comprehensive aerosol–cloud PPE using the ICON-HAM model to systematically investigate parametric uncertainty in ERF_aer and key aerosol and cloud properties. As a relatively new aerosol–climate model, ICON-HAM particularly benefits from the PPE framework that supports parameter optimisation and targeted model development. Our PPE comprises 383 simulations in which 42 parameters are perturbed, spanning 13 categories that include emissions, aerosol/cloud properties, and aerosol/cloud processes. To move beyond the discrete PPE simulations, we use machine-learning statistical emulation (Gaussian Process) to continuously sample the 42-dimensional parameter space, enabling quantitative assessment of parametric uncertainty in ERF_aer and aerosol-cloud properties. We further evaluate model performance against satellite observations (SPEXone and MODIS) and identify potential causes of both model uncertainty and model-observation bias.

This study aims to address the following questions. (i) Which properties and processes drive uncertainties in ERF_aer and in key aerosol and cloud properties? (ii) What present day (PD) observable quantities can best constrain the preindustrial (PI) to PD forcing? (iii) To what extent can model-observation bias be reduced through parametric optimisation? In addition, we also examine a technical aspect of the PPE experimental design: what is the optimal r_sp value that balances computational cost, data-storage demands, and emulator accuracy?

This paper is organised as follows. We first describe the model configuration, PPE setup, emulator, and satellite observations in Sect. 2. Section 3.1 presents the emulator performance and the determination of the optimal r_sp. Section 3.2 and 3.3 quantify and attribute parametric uncertainties in ERF_aer and key aerosol and cloud properties. Section 3.4 presents the sensitivity of model outputs to perturbed parameters, while Sect. 3.5 evaluates model against satellite observations and relates model biases to parametric uncertainties. In Sect. 4, we discuss the main findings and their implications.

2 Methods and data

2.1 Model description

The global aerosol–climate model ICON2.6.4–A–HAM2.3 is a combination of the icosahedral non-hydrostatic global atmosphere model ICON–A version 2.6.4 (Giorgetta et al., 2018), with the aerosol microphysics module HAM version 2.3 (Stier et al., 2005; Tegen et al., 2019). The coupling of ICON-A and HAM models is described in Salzmann et al. (2022). For brevity, we refer to ICON2.6.4–A–HAM2.3 as ICON–HAM throughout this paper. The model can be configured for different resolutions. In this study, we use the triangular R2B4 horizontal grid (approximately 160 km resolution) with 47 hybrid-sigma vertical layers extending to 80 km.

ICON–HAM employs the Max Planck Institute physics package inherited from the ECHAM6 model (Stevens et al., 2013; Giorgetta et al., 2018; Mauritsen et al., 2019). Radiative transfer calculations are performed using the PSrad scheme (Pincus and Stevens, 2009), an adaptation of the Rapid Radiative Transfer Model (RRTM; Mlawer et al., 1997) optimised for general circulation models (Iacono et al., 2008; Morcrette et al., 2008). Moist convection is parameterised using a mass-flux scheme (Tiedtke, 1989; Nordeng, 1994), allowing either shallow, mid-level, or deep convection to develop within the atmospheric column. Entrainment rates for different convection regimes and for downdrafts are prescribed as model parameters. The turbulent diffusion is parametrised using a total turbulent energy (TTE) scheme (Mauritsen et al., 2007; Pithan et al., 2015), where TTE is the sum of turbulent potential energy and turbulent kinetic energy (TKE), with TKE used to compute the vertical velocity for cloud droplet activation. An overly strong vertical diffusion was identified (Franke and Giorgetta, 2024) and subsequently corrected in our version, through three bug fixes related to the total turbulent energy, the vertical interpolation of variables from full to half levels, and the maximum mixing length scale (Klamt et al., 2025). The land surface model is a simplified version of JSBACH4 (Raddatz et al., 2007) – JSBACH4-lite (Salzmann et al., 2022), providing only static fields of land surface properties. Following Sundqvist et al. (1989), stratiform cloud cover is parameterised as a function of relative humidity (RH). Cloud droplet and ice crystal mass and number concentrations are treated prognostically using a two-moment cloud microphysics scheme (Lohmann et al., 2007; Lohmann and Neubauer, 2018). Droplet activation is computed based on a recent modified version of the Abdul-Razzak and Ghan (2000) activation scheme (Ghosh et al., 2025). To account for the absence of nitrate aerosols and the underestimation of secondary organic aerosols, a lower bound of 10 cm⁻³ for cloud droplet number concentration (CDNC) was imposed; see Lohmann and Neubauer (2018) and Neubauer et al. (2019) for discussion.

The HAM2.3 module simulates a range of aerosol processes, including emission, nucleation, transport, transformation and removal (Stier et al., 2005; Tegen et al., 2019). Aerosol size distributions and composition are represented by the modal M7 aerosol model (Vignati et al., 2004), with major aerosol species (sulfate, organic carbon, black carbon, sea-salt, and mineral dust) treated as internally mixed particles in seven log-normal modes. The seven modes cover four size ranges (nucleation, Aitken, accumulation and coarse) and two hygroscopic states (hydrophilic and hydrophobic). Aerosol size distribution, hygroscopicity, and optical properties evolve throughout their life cycle, hence influencing interactions with radiation and clouds (Salzmann et al., 2022).

Cloud properties derived directly from model outputs can differ markedly from those obtained using a satellite simulator, potentially misleading model evaluation, constraint, and thus development (Saponaro et al., 2020). To enable a better model-satellite comparison, the Cloud Feedback Model Intercomparison Project (CFMIP) Observation Simulator Package (COSP) version 2.0 (Bodas-Salcedo et al., 2011; Swales et al., 2018) was recently implemented into the ICON–HAM model to derive MODIS-style cloud diagnostics (Henkes et al., 2025) and is used here. By accounting for the characteristics and limitations of the MODIS instrument (i.e. variable definitions, retrieval assumptions, sensor limitations, spatial resolution), COSP-simulated fields can be meaningfully compared to MODIS cloud retrievals.

2.2 Experimental setup

2.2.1 Simulations

To diagnose aerosol forcing and associated parametric uncertainty, we conducted two PPEs for PD and PI aerosol conditions. Each PPE consists of 383 simulations with different parameter combinations (Sect. 2.2.2). We selected the earliest 12-month period (1 August 2024 to 31 July 2025) for which new observations from both the PACE and EarthCARE satellites are available. The PD and PI ensembles differ only in prescribed anthropogenic and biomass burning aerosol emissions. Hence, the PI and PD simulations are driven by the same boundary conditions over the simulating period, including the greenhouse gases and atmospheric oxidants concentration from the Shared Socioeconomic Pathways SSP245 scenario as well as monthly sea surface temperature and sea ice cover from the ERA5 reanalysis. Both PI and PD horizontal winds throughout the atmospheric column are nudged towards ERA5 winds in the PD period, with a relaxation timescale of 6 h. To our knowledge, this is the first study applying nudging in ICON-A-(HAM). Evaluation of the nudging setup indicates that it achieves a nudging strength similar to that used in ECHAM–HAM. We also compared the aerosol forcing in twin simulations with full-column and above-2 km nudging and found that the global-mean ERF_aer is insensitive to nudging extent (Fig. S1 in the Supplement), consistent with previous studies (Kooperman et al., 2012), while over regions with strong biomass-burning emissions, allowing the lower atmosphere to respond to aerosol absorption leads to a slightly more positive ERF_ari.

The PD ensemble uses anthropogenic emissions from the Community Emissions Data System (CEDS v2025_04_18) (Hoesly et al., 2018) and biomass burning emissions from the Global Fire Assimilation System (GFAS) for the simulation period, whereas the PI ensemble uses CMIP6 emissions for the year 1850. Since CEDS data are only available through 2023, emissions from 2022–2023 are assumed to represent also the 2024–2025 period in the PD simulations. This approximation is not expected to substantially affect the results, because emissions are unlikely to change significantly over two years based on recent rates of change in anthropogenic emissions (Hoesly et al., 2018), and moreover our PPE framework is able to scale emissions accordingly.

Prior to the main PPE simulations, a 1-year run with the default model setup (CTRL) was performed to spin up the model. The PPE runs were then initialized from the end of the 1-year CTRL run. Each PPE member underwent an additional 3-month spin-up followed by a 12-month run used for the final analysis.

ERF_aer is calculated as the difference in net radiative flux at the top-of-atmosphere (TOA) between the PD and PI simulations. Following Ghan (2013), ERF_aer is further decomposed into ERF_aci and ERF_ari, through double calls to the radiation calculation with and without aerosols.

Table 1 Descriptions of the 26 aerosol parameters and their distributions used in the emulations, grouped by physical process. “Default” denotes the baseline values used in the default run (CTRL). Parameters labeled “Rel” are scaling factors applied to modelled variables produced by the original parameterisation, also known as tuning parameters, whereas those labeled “Abs” represent absolute parameter values. Generalized trapezoidal distributions are parameterised in the form (a, b, c, d, κ_l, κ_r), with a, b, c, and d specifying the lower bound, start of plateau, end of plateau, and upper bound, respectively, and κ_l, κ_r governing the curvature of the left and right boundary slopes, respectively. A value of κ = 1 yields linear slopes, while κ > 1 (κ < 1) produces concave (convex) curvature.

Table 2 Same as Table 1 but for the 16 cloud parameters.

2.2.2 Perturbed parameters

Perturbed parameters, their ranges, and probability distribution functions (PDFs) were initially determined based on the elicitation of model domain experts and prior PPE studies (Regayre et al., 2018; Yoshioka et al., 2019; Eidhammer et al., 2024; Bhatti et al., 2026), and subsequently refined through preliminary screening tests. First, we conducted 3-month one-at-a-time perturbation experiments to screen out parameters with little impact on key aerosol and cloud properties and hence ERF_aer (not shown). Then, we perform a 3 d PPE to make sure that the parameter combinations from the initial ranges would not lead to a large fraction of failed simulations. The screening revealed that failures occur when conv_cmfctop and conv_entrscv are simultaneously high; their upper bounds were therefore reduced. Following this screening procedure, a total of 42 parameters and their ranges were determined for perturbation in the PPEs, including 26 aerosol-related parameters (Table 1) and 16 cloud-related parameters (Table 2).

The aerosol-related parameters represent processes/properties such as aerosol emissions, nucleation, aging, deposition, chemistry, and optical and hygroscopic properties (Table 1). These parameters are largely consistent with those used in the ECHAM–HAM aerosol PPE by Bhatti et al. (2026), but we extend the perturbations to better represent size- and species-dependent effects. Specifically, aerosol emissions are perturbed by sector (anthropogenic and biomass burning) and by species (SO₂, BC, and OC), and sea-salt emissions are perturbed by mode (accumulation and coarse). In addition to scaling overall emission fluxes, we perturb key parameters in the sea-salt and dust emission parameterisations (emi_ss_expo and emi_du_fricv), which control not only the emission flux but also the size of emitted particles. The perturbed aerosol parameters affect aerosol sources and sinks, as well as their ability to absorb and scatter radiation and to serve as CCN. Cloud-related parameters describe processes/properties such as droplet activation, cloud microphysics, convection, cloud cover, cloud optical properties, and turbulence (Table 2). These parameters influence cloud properties and the boundary-layer humidity distribution, and are key controls on the model's radiation balance. Together, the aerosol and cloud parameters are expected to modulate ERF_aer through ARI and ACI pathways, or both.

Maximin Latin Hypercube Sampling (LHS) was employed to create parameter combinations spanning the 42-dimensional parameter space (McKay et al., 1979). In this study, two parameter distributions are input to LHS for different purposes: one to generate the original PPEs to train emulators (Fig. S2) and another to construct the emulated PPEs (Fig. S3). For the original PPEs, most parameters are sampled from uniform distributions, while a subset follows trapezoidal distributions with a uniform central region and sloped edges representing relatively unlikely parameter ranges (Fig. S2). This design preserves dense, uniform-like sampling in likely parameter space while efficiently extending the parameter ranges with minimal additional computational cost. For the emulated PPEs that are used for quantifying the parametric uncertainty, the trapezoidal distributions are adopted (Tables 1 and 2, Fig. S3). Compared to Gaussian distributions, this approach applies less tight expert-elicited constraints, avoids central oversampling, thereby preserving a broader uncertainty range. Note that the ranges displayed in Figs. S2 and S3 are the physical values, whereas the LHS operates on normalized ranges to determine the parameter combinations.

In total 400 parameter combinations were designed using LHS of the 42 perturbed parameters. Of these, 383 full-year ensemble members were successfully produced for both PI and PD periods: one corresponds to the default ICON–HAM configuration (CTRL), while the remaining 382 reflect the sampled parameter combinations. This yields a simulations-to-parameters ratio (r_sp) of approximately 9, which is relatively high compared to previous PPE studies. For instance, Regayre et al. (2018), Yoshioka et al. (2019), and Eidhammer et al. (2024) used r_sp values of around 6, while Bhatti et al. (2026) and Nugent et al. (2026) adopted a higher value of 10. The choice of r_sp is crucial because it directly affects both computational cost and emulator accuracy, yet its optimal values, balancing accuracy, computational cost, storage requirements, remain unclear. This question will be addressed in Sect. 3.1.

2.3 Gaussian process emulator

The Gaussian Process (GP) emulator provides a flexible, non-parametric statistical framework that can efficiently approximate complex, nonlinear relationships between model parameters and outputs (O'Hagan, 2006), and has been applied extensively in studies of climate model uncertainty (Lee et al., 2011; Yoshioka et al., 2019; Bhatti et al., 2026; Regayre et al., 2026; Prévost et al., 2026). To enable continuous sampling of the 42-dimensional parameter space, we constructed GP emulators trained on the original PPEs with 383 ensemble members using the toolkit Earth System Emulator package (ESEm; Watson-Parris et al., 2021a), with the 42 perturbed parameters as predictors and selected diagnostics (aerosol, cloud and radiation properties) as targets. Here, the emulators are trained for each model grid using monthly mean data. Once trained, the GP emulators can provide rapid predictions at previously unsampled parameter combinations within the defined parameter bounds.

Figure 1 Workflow for evaluating emulator performance for different simulations-to-parameters ratios (r_sp), based on repeated partitioning and resampling of 383 PPE members. Here, n denotes the number of PPE members in each resampled subset and is determined as n = 42 × r_sp, where 42 is the number of perturbed parameters.

Using the trained emulators together with the parameter distributions defined in Tables 1 and 2, we generated the emulated PPEs for PI and PD, each containing 200 000 model variants. We also tested a larger size of 2 million variants and found no change in conclusions on uncertainty quantification and attribution. The resulting ensembles were then used for variance-based decomposition of parametric uncertainty (see Lee et al. (2012) for details of the method) in Sect. 3.3 and sensitivity analysis in Sect. 3.4.

2.4 Workflow for determining the optimal r_sp

The high r_sp (∼ 9) in our PPE provides a testbed to assess emulator performance as a function of r_sp, as illustrated in Fig. 1. The full PPE ensemble (n = 383) is first partitioned into training and test datasets, with 30 members reserved for testing. To ensure representative sampling, this partitioning is randomly repeated 30 times. For each partition, the training dataset is further resampled to generate subsets corresponding to r_sp ranging from 2 (n = 84 members) to 8 (n = 336 members). The resampling procedure preserves the marginal distribution of each parameter so that the resulting subsets remain statistically consistent with the original ensemble. To ensure statistical robustness, each resampling process is also randomly repeated 30 times. The resulting subsets are used to train the emulators, and their predictions are evaluated against the independent test dataset. The 30 × 30 repetitions capture the uncertainty associated with the partitioning and resampling procedures.

2.5 Observations

Aerosol properties are obtained from SPEXone aboard NASA's PACE mission, launched in February 2024 (Werdell et al., 2019). The SPEXone instrument is a multi-angle polarimeter developed by SRON for high-accuracy aerosol retrievals (Hasekamp et al., 2019a; van Amerongen et al., 2019). It measures both polarization and intensity at five viewing angles over the spectral range 385–770 nm. Aerosol retrievals from the Remote sensing of Trace gas and Aerosol Products (RemoTAP) algorithm (Fu et al., 2020; Lu et al., 2022; Hasekamp et al., 2024), produced by the SRON SPEXone scientific processing line (Fu et al., 2025), are employed. Daily level-2 retrieval products are aggregated onto a 1° × 1° resolution for analyses. From these data, we use aerosol optical depth (AOD) and absorption AOD (AAOD) at 550 nm, and Ångström exponent (AE) calculated from AOD at 440 and 670 nm. The same wavelengths are also used for the corresponding model outputs for comparison in Sect. 3.5.1. All quantities have been validated against AErosol RObotic NETwork (AERONET) observations and show good agreements, with root mean square errors of 0.05 (0.04) for AOD, 0.021 (0.015) for AAOD, and 0.366 (0.329) for AE over land (ocean) (Fu et al., 2026).

Cloud properties are derived from the MODIS COSP daily level-3 dataset (MCD06COSP; Pincus et al., 2023). This product relies on the MODIS Collection 6.1 processing stream (Platnick et al., 2017) and is specifically designed to facilitate comparison of MODIS cloud retrievals with climate models, in this study, the cloud diagnostics from ICON–HAM with the MODIS COSP simulator (see Sect. 2.1). From this dataset, we use cloud optical thickness (COT), cloud droplet effective radius (CER), liquid water path (LWP) and liquid cloud fraction (CF_liq), where CER and COT are used to compute CDNC based on the adiabatic approximation (Quaas et al., 2006).

Figure 2 Comparison of globally averaged GP emulator predictions with test data (truth) for (a) AOD, (b) AE, (c) CDNC, (d) LWP, (e) CF, and (f) ERF_aer. Red lines indicate linear regression fits, and the black line denotes the 1:1 reference. The number of data points (N), root mean square error (RMSE), and coefficient of determination (R²) are also shown. Similar comparisons for ERF_aer by component (ACI and ARI) are shown in Fig. S4.

To facilitate meaningful model–observation comparison, the 3 h ICON–HAM outputs are co-located with SPEXone and MODIS COSP observations in space and time, such that each comparison pairs the satellite footprint and overpass time with the nearest model grid cell and closest model time step. The collocation is efficiently handled by Community Intercomparison Suite (CIS) – a tool designed for model-observation comparisons (Watson-Parris et al., 2016). CIS supports various model and observational datasets, but lacks native support for MODIS COSP and SPEXone; we therefore implemented custom plugins to enable their collocation within the CIS workflow.

3 Results

3.1 Emulator results

Comprehensive quantification and attribution of model uncertainties rely on statistical emulators, which rapidly predict outputs for any parameter combinations within the defined ranges, without the need for further computationally expensive model runs. The emulator skill is evaluated prior to subsequent analyses.

3.1.1 Emulator performance

Emulator performance is assessed using a split-sample validation, with 80 % of simulations (307) used for training and 20 % (76) for testing. Figure 2 shows globally averaged comparisons of GP predictions against truth, for key aerosol (AOD and AE) and cloud (CDNC, LWP, and cloud fraction (CF)) properties governing ARI and ACI, as well as ERF_aer. Overall, the GP predictions closely follow the 1:1 reference line. The best performance is found for aerosol properties, with slopes of 0.89 for AOD and 0.85 for AE. Cloud properties show slightly lower slopes (< 0.85), while ERF_aer deviates most strongly from the 1:1 line, with a slope of 0.67. The coefficients of determination (R²) follow a similar behaviour. Root mean square error (RMSE) is small relative to the magnitude of each variable, demonstrating good absolute accuracy. The progressive degradation in emulator performance from aerosol properties to cloud properties and ultimately to ERF_aer likely reflects the increasing influence of nonlinear processes. For ERF_aer in particular, spatial differences between PI and PD emissions further complicate the response because global parameter perturbations cannot capture this structure. Nevertheless, even for ERF_aer, an R² of 0.76 still indicates robust predictive skill.

Figure 3 Comparison of emulated and reference AOD averaged over 76 test members. Shown are 2-D maps of (a) “truth”, (b) absolute bias, (c) relative bias (%), and (d) coefficient of determination (R²) between predictions and test data.

Figure 4 Emulator performance as a function of the simulations-to-parameters ratio (r_sp), measured using regression slope, coefficient of determination (R²), and spatial correlation (r_spatial) between emulator predictions and truth for AOD, AE, CDNC, LWP, CF, ERF_aer, ERF_aci and ERF_ari. Lines show means across the 30 × 30 emulation experiments illustrated in Fig. 1, and shading indicates corresponding standard deviations, representing the uncertainty in partitioning and resampling procedures. Two grey dashed lines denote the uncertainty range at r_sp=8, and the yellow star marks the identified optimal ratio.

Figure 3 shows 2-D maps of comparison between predicted and reference AOD as an example; similar behaviour is found for other variables (not shown). Overall, the emulator faithfully reproduces the spatial patterns of all quantities. Despite visible regional variations of absolute and relative differences, the overall biases remain very small, with global-mean relative biases below 1.5 % and local maxima under 4 %, which are substantially smaller than typical model and observational uncertainties (Gliß et al., 2021). It is interesting to note that regions with lower R² align closely with areas of strong convective precipitation, particularly over the tropical oceans near the Intertropical Convergence Zone, the Amazon, southern Africa, and South Asia over land. The association between reduced emulator performance and strong convective precipitation implies that the emulator may struggle to capture the strong nonlinear interactions among processes, such as emissions, convection, and wet scavenging. Nevertheless, this has only a limited impact on the overall emulator performance, with the minimum R² remaining high (> 0.7).

3.1.2 Optimal simulations-to-parameters ratio

The excellent emulator performance of our PPEs is not surprising given the relatively high simulations-to-parameters ratio in our PPE (r_sp ∼ 9). However, such a large r_sp also implies substantial computational cost and data-storage demands, posing challenges for large PPEs or multi-model PPE (MMPPE) intercomparison project. Identifying the optimal r_sp, defined as the minimum value that yields acceptable performance (within the uncertainty range of metrics at r_sp=8), is therefore essential for guiding future PPE design.

Figure 5 Probability density functions (PDFs) of global mean ERF_aer, ERF_aci and ERF_ari for 1850–2025 from ICON–HAM PPE, together with forcing ranges from this study and other PPE studies, multi-model ensembles (MME), and IPCC assessments. The ICON–HAM estimates (black) are based on 200 000 emulated model variants sampled based on the expert-elicited prior PDFs of parameter values shown in Fig. S3. AER (pink) and CLD (blue) denote ICON–HAM variants with perturbations applied only to aerosol (Table 1) or cloud (Table 2) parameters, respectively, while all other parameters are fixed at their prior medians. The 90 % credible ranges of three ICON–HAM experiments are shown with filled circles marking the means, with numerical values shown alongside. The 90 % credible ranges of other studies are displayed in gray (full range for CAM6), with means (or medians) indicated by open circles, including IPCC AR5 (Myhre et al., 2013b; 1750–2011), AR6 (Forster et al., 2021; 1750–2019), CMIP5 (Zelinka et al., 2014; 1860–2000), CMIP6 (Smith et al., 2020; 1850–2014), ACCMIP (Shindell et al., 2013; 1850–2000), AeroCom IND3 (Gryspeerdt et al., 2020; 1850–2010), and PPE studies from CAM6 (Eidhammer et al., 2024), UKESM1–A (Regayre et al., 2023; 1850–2014), HadGEM–UKCA (Regayre et al., 2018; Yoshioka et al., 2019; 1850–2008), and ECHAM–HAM (Bhatti et al., 2026; 1850–2010). HadGEM–UKCA with suffixes Y and R correspond to Yoshioka et al. (2019) and Regayre et al. (2018), respectively; R_AER and R_ATM denote HadGEM–UKCA_R variants with only aerosol or atmospheric parameters perturbed. Means or medians were not reported in the CAM6 and UKESM1–A studies.

Following the procedure illustrated in Fig. 1, we assess emulator performance as a function of r_sp for key aerosol and cloud properties, as well as aerosol forcing (Fig. 4). Consistent with Fig. 2, the emulator performs best for aerosol quantities (AOD and AE), followed by cloud quantities (CF, LWP, and CDNC), with lower skill for aerosol forcing. This ranking is robust across different r_sp choices and performance metrics. Note that the relatively limited performance of ERF_aer mainly stems from ACI. Emulating ERF_aer is challenging because it highly depends on spatial differences between PI and PD emissions, as well as on the minimum CDNC (CDNC_min) imposed in the model (Hoose et al., 2009). This CDNC_min sets a hard bound on several processes (activation, autoconversion, entrainment mixing, etc.), thereby introducing pronounced nonlinear interactions between model parameters and outputs that are difficult to capture in a PPE (Bhatti et al., 2026).

Overall, the performance metrics (particularly the slope and R²) first increase with r_sp and then begin to level off beyond a certain threshold. This threshold is around 5 for aerosol and cloud properties and ERF_ari, but higher (∼ 8) for ERF_aer and ERF_aci. Nevertheless, even for ERF_aer and ERF_aci, the metric values at r_sp=5 (yellow star) remain within the uncertainty range of those at r_sp=8, and the emulator bias and RMSE decrease by only 0.01 and 0.02 W m⁻² when increasing r_sp from 5 to 8 (Fig. S5). This indicates that r_sp=5 already achieves near-optimal emulator performance while substantially reducing computational cost by 40 %. We therefore consider r_sp=5 as a compromise between accuracy and efficiency, even for ACI-related analyses. Further improvements in emulator performance should therefore focus on refining the machine-learning model (e.g., kernel tuning) rather than simply increasing the ensemble size.

3.2 Uncertainty in aerosol radiative forcing

3.2.1 Global mean

Figure 5 shows probability density functions (PDFs) of global mean ERF_aer, ERF_aci and ERF_ari over 1850–2025, derived from 200 000 emulated variants generated based on the expert-elicited parameter distributions shown in Fig. S3. The mean ERF_aer is −1.04 W m⁻², with a 5–95th percentile spanning −1.42 to −0.65 W m⁻² (Fig. 5a). The ACI component, ERF_aci, shows a very similar estimate of −1.03 [−1.37 to −0.67] W m⁻² with a slightly narrower range (Fig. 5b), indicating that most of the uncertainty in ERF_aer stems from ACI. In contrast, the ARI component ERF_ari is small with a mean of 0.04 W m⁻² but a somewhat wide range of −0.13 to 0.23 W m⁻² spanning both negative and positive values (Fig. 5c). Unlike ERF_aci, which is uncertain only in magnitude, ERF_ari is uncertain in both magnitude and sign.

To isolate the contributions from the aerosol and cloud parameter groups, we use the emulator to vary the parameters in one group over their uncertainty ranges while fixing all parameters in the other group at their prior medians (Fig. 5). Perturbing aerosol parameters only (AER; pink) results in a 5–95th percentile range of [−1.43 to −0.71] W m⁻², close to the total uncertainty range, whilst perturbing cloud parameters only (CLD; blue) yields a narrower range of [−1.18 to −0.9] W m⁻², with the width reduced by 64 % relative to total range, suggesting that aerosol parameters dominate the overall uncertainty. Comparing ACI and ARI (Fig. 5b,c) shows that cloud parameters influence ACI more strongly, although their effect remains weaker than that of aerosol parameters. The specific process-level sources of uncertainty will be identified in Sect. 3.3.

Figure 5 also compares the ICON–HAM PPE results with other PPEs using different climate models, multi-model ensembles (MMEs), and IPCC assessments. Comparison with different PPE studies shows that forcing estimates may vary across models depending on model structure, parameter perturbations considered, experimental setup (e.g., nudged versus free runs and the simulation period). For ECHAM–HAM, which perturbs a similar set of parameters (but only aerosol parameters) within a comparable parameterisation design (both using the HAM aerosol module) and utilises same nudging strategy (the whole column wind nudging), the resulting forcing range is very similar to that of our PPE but slightly narrower (Bhatti et al., 2026). This is partly because cloud-related uncertainty is not sampled, and also because the ECHAM–HAM PPE adopts Gaussian parameter distributions, which sample the central region more densely than the trapezoidal distributions used in our study, leading to a narrower range. The somewhat weaker ERF_ari in our PPE may result from the inclusion of above-cloud aerosols in the forcing calculation, which influences the balance between ACI and ARI (Ghan, 2013).

Models with greater structural differences from ICON–HAM show larger discrepancies in forcing estimates. Four PPE experiments with HadGEM–UKCA (Regayre et al., 2018; Yoshioka et al., 2019) show similar magnitudes and uncertainty ranges but yield more negative ERF_aer than the ICON–HAM PPE. In addition to differences in model structure and parameter perturbations, the stronger forcing may also arise from the different PD reference period used, as a more polluted year (2008) was adopted in those experiments. According to the IPCC AR6 (Forster et al., 2021), ERF_ari was 50 % stronger in 2008 than in 2019, and ERF_aci was 21 % stronger. Although all uncertainties shown here are from PPEs that are not yet constrained using observations, different parameter sampling strategies applied to the emulator implicitly lead to different levels of expert elicitation. ICON–HAM, HadGEM–UKCA, and ECHAM–HAM use expert-elicited parameter distributions (trapezoidal for ICON–HAM and HadGEM–UKCA, Gaussian for ECHAM–HAM), which represent prior expert beliefs about parameter likelihoods. In contrast, CAM6 and UKESM1–A use uniform distributions based only on expert-specified parameter ranges, which sample the full parameter space uniformly and therefore produce much wider uncertainty ranges. Note that the winds in ICON–HAM and ECHAM–HAM PPEs are more tightly nudged compared with the free-running (CAM6) and above-boundary-layer nudging (UKESM1-A and HadGEM–UKCA) PPEs, thereby limiting atmospheric adjustments and resulting in weaker forcing (particularly for ACI). Comparisons among PPEs suggest that different PPEs are not readily comparable. This highlights the need for a coordinated multi-model PPE (MMPPE) with a standardized protocol across models, including common perturbed parameters, parameter distributions, and experimental setups, to better disentangle structural and parametric uncertainties (Carslaw et al., 2026).

The spread in MMEs (CMIP5, CMIP6, ACCMIP and AeroCom) mainly reflects structural diversity and internal variability, because each model uses its own tuned parameter set without sampling parameter uncertainty within the model. Interestingly, the ERF_aer uncertainty ranges from MMEs and PPEs are of similar magnitude, possibly implying a comparable contribution from structural uncertainty and internal variability versus from parametric uncertainty. A more robust attribution, however, would require further coordinated MMPPE efforts. IPCC assessments consider a broader range of evidence than MMEs alone, including individual model studies, uncertainties from missing or poorly represented processes, independent observational, and energy-balance constraints (Myhre et al., 2013b; Forster et al., 2021), thus resulting in broader ERF_aer uncertainty ranges.

Figure 6 (a–c) Mean and (d–f) parametric uncertainty (2σ across the ensemble members) of ERF_aer (left), ERF_aci (middle), and ERF_ari (right) over 1850–2025, derived from 200 000 emulated model variants. Panels (g)–(i) show the percentage contribution of aerosol parameters to the total uncertainty (%), with cloud parameters accounting for the remainder. These model variants are sampled according to the expert-elicited parameter distributions shown in Fig. S3. Top-left numbers indicate global mean values.

3.2.2 Spatial distributions

Figure 6 shows the spatial distributions of ERF_aer, ERF_aci, and ERF_ari and their corresponding uncertainties, along with the relative contributions of aerosol parameters to these uncertainties. Similar to the global-mean results (Fig. 5), ERF_aer and ERF_aci exhibit very similar spatial patterns, indicating ERF_aer is primarily driven by ACI. The strongest forcings (< −5 W m⁻²) are located in industrialized continental regions, particularly Asia, Europe, and North America, followed by marine stratocumulus regions over the Atlantic, North Pacific, and the southeastern Pacific off the west coast of South America. The spatial patterns of ERF_aer and ERF_aci broadly agree with those from the previous ECHAM–HAM PPE (Bhatti et al., 2026), which uses similar aerosol and cloud parameterisations and comparable aerosol parameter perturbations. Our results, however, show weaker forcing over China and stronger forcing over India, likely due to anthropogenic emission reductions in China and increases in India between the respective simulation years (2010 for ECHAM–HAM and 2024–2025 for ICON–HAM). Notably, ERF_aer and ERF_aci turn positive over the UK because its anthropogenic aerosol emissions in 2024–2025 have dropped below 1850 levels (Fig. S6). Overall, regions exhibiting strong ERF_aer and ERF_aci also tend to show larger associated uncertainties (Fig. 6d and e).

The near-zero global mean ERF_ari is a result of the compensation of negative and positive regional contributions (Fig. 6c). Negative ERF_ari is prominent over the northern Indian Ocean, influenced by the transport of light-scattering sulfate aerosols from India. Positive ERF_ari occurs primarily where absorbing carbonaceous aerosols are densely emitted, or where they are transported over high-albedo cloud decks, including the northern India, South America, and North Atlantic. Similar effects of the above-cloud aerosols on ERF_ari were also identified in the HadGEM PPE (Regayre et al., 2018), albeit with different spatial patterns. Additionally, a positive ERF_ari is also seen over the Saharan region but through a different mechanism, namely a reduction in dust aerosols resulting from their interactions with sulfate (Neubauer et al., 2019). Sulfate can coat dust particles, transferring them from insoluble to soluble modes and thereby enhancing their wet deposition, which more strongly reduces PD dust. At the same time, dust provides surfaces for sulfate condensation, suppressing new particle formation. Together, these processes contribute to a positive ERF_ari over dust-dominated regions. Similar positive ERF_ari induced by sulfate–dust interactions has been reported in other models, including GFDL (Fan et al., 2004), GISS (Bauer et al., 2007), and ECHAM–HAM (Neubauer et al., 2019). The weak ERF_ari over China likely results from the partial cancellation between warming from absorbing BC aerosols and cooling from scattering sulfate aerosols, both of which have substantial anthropogenic emissions in this region (Wu et al., 2024). This is further supported by the sensitivity analysis presented later (Fig. S19), where ERF_ari shows comparable but opposite sensitivities to their emission parameters over land (−2.6 to emi_antns_so2 and +2.4 to emi_ant_bc). Substantial uncertainty is found in regions strongly influenced by industrial and biomass-burning emissions, including Asia, the boreal regions of North America and their downwind oceans, as well as parts of South America and Africa (Fig. 6f).

The lower panels of Fig. 6g–i illustrates the relative contribution of aerosol parameters to total uncertainty, with cloud parameters accounting for the remainder. Brownish colors denote aerosol control of uncertainty, and greenish colors cloud control. For ERF_aer and ERF_aci, aerosol parameters dominate over land and across the Northern Hemisphere, where aerosol loadings are higher, whereas the contribution of cloud parameters is limited to tropical regions characterized by intense convection. In contrast, ERF_ari is largely aerosol-dominated across the globe, with only a modest contribution from cloud parameters in convective regions, likely linked to wet scavenging. A detailed process-level attribution is provided in the next section.

Figure 7 Percentage contributions of individual parameters to uncertainties in (a) ERF_aer, (b) ERF_aci, (c) ERF_ari over globe and eighteen climatic regions as shown Fig. A1, derived from 200 000 emulated variants. Regions are ordered from north to south, separately for land (left) and ocean (right): Arctic (ARC), North America (NAM), Europe (EUR), Asia (ASI), Africa (AFR), South America (SAM), Oceania (OCE), Antarctica (ANTA); Arctic Ocean (ARCO), North Pacific Ocean (NPO), North Atlantic Ocean (NAO), Tropical Pacific Ocean (TPO), Tropical Atlantic Ocean (TAO), Tropical Indian Ocean (TIO), South Pacific Ocean (SPO), South Atlantic Ocean (SAO), Southern Indian Ocean (SIO), and Antarctic Ocean (ANTAO). Parameters contributing less than 3 % of the total uncertainty are not shown. The corresponding spatial maps of uncertainty caused by individual parameters are shown in Figs. S7–S9. Parameters in the legend are ordered from aerosol- to cloud-related processes (reading top-down, then left-to-right), and the stacked colors mirror this sequence from bottom to top.

3.3 Sources of uncertainty

3.3.1 Aerosol radiative forcing

Figure 7 presents the global and regional percentage contributions of individual parameters to uncertainties in ERF_aer, ERF_aci, and ERF_ari, with region abbreviations defined in the caption and Fig. A1. The corresponding contributions to uncertainties in global-mean forcing are provided in Table 3.

Table 3 Percentage contributions to global annual mean forcing uncertainties by individual parameters over 200 000 model variants, corresponding to Fig. 7. Parameters contributing more than 2 % are highlighted in bold, while those contributing less than 0.01 % to ERF_aer are omitted.

ERF_aer

Globally, the ERF_aer uncertainty is dominated by emission parameters. Anthropogenic SO₂ (emi_antns_so2), natural DMS (emi_dms), and emitted particle size from biomass burning and fossil fuel sources (emi_cmr_bb and emi_cmr_ff) contribute similarly (25 %, 27 %, and 28 %), together accounting for nearly 80 % of the total uncertainty. This agrees with results from HadGEM (Regayre et al., 2018) and ECHAM–HAM (Bhatti et al., 2026) PPEs, but shows a larger influence of emitted-size parameters, particularly for biomass burning, likely due to stronger boreal fires in our simulation year (Zhong et al., 2024). It should also be noted that a key difference in the representation of ACI between this study and previous studies lies in the droplet activation. The HadGEM (Regayre et al., 2018) and ECHAM–HAM (Bhatti et al., 2026) PPEs employed the default ARG scheme (Abdul-Razzak and Ghan, 2000), whereas we uses a modified version by Ghosh et al. (2025). The higher sensitivity of the revised ARG to Aitken mode (Ghosh et al., 2025) could also help explain the increased sensitivity to particle size found here. The parameters emi_antns_so2 and emi_dms cause the uncertainty by controlling the magnitude of the PI–PD aerosol change. In particular, the role of DMS highlights the importance of the PI aerosol background. Similarly, uncertainties in biogenic volatile organic compound (VOC) emissions and processing may affect the PI aerosol background, too, but are not explored here because full chemistry is not activated in the current model configuration. For a given emitted mass, emi_cmr_bb and emi_cmr_ff determine the emitted particle number, thereby influencing CCN concentrations and cloud properties, and hence ACI. In comparison, the convection parameter conv_cprcon, controlling the conversion of convective cloud water to rain, contributes a smaller fraction of 6 % (Table 3). It does not directly influence the anthropogenic aerosol perturbation but modifies the cloud radiative state by altering cloud amount, which sets the background cloud state for aerosol–cloud interactions.

It is unsurprising that the emission parameters dominate the ERF_aer uncertainty, given their first-order control on PD–PI aerosol differences. Once anthropogenic, biomass burning, and DMS emissions are fixed at their prior medians, cloud parameters dominate the remaining uncertainty (76 %; Fig. S10), primarily through convection (conv_cprcon, conv_entrpen, and conv_cmfctop), cloud cover (cld_cov_crs), and auto-conversion (micro_ccsaut). The remaining uncertainty mainly from aerosol properties (refractive indices) and processes (SO₂ chemistry and wet deposition).

Because of the short lifetime of tropospheric aerosols, the contributions of these parameters to ERF_aer uncertainty vary substantially across regions. In the main industrial regions (NAM, EUR, and ASI) and their outflow regions (ARC, NPO, and NAO), where the influence of biomass burning aerosols is also present, emi_cmr_bb appears to play a dominant role (Figs. 7a and S7), with contributions ranging from about 25 % in NAM to as high as 55 % in EUR. In contrast, in biomass-burning affected regions (AFR, SAM, TAO and SPO) in the Southern Hemisphere, where anthropogenic sulfate is less abundant, emi_cmr_bb's contribution is much smaller (< 10 %). This contrast reflects the role of aerosol aging: abundant soluble materials (e.g., sulfate) can coat insoluble biomass-burning particles, transforming them into more hygroscopic particles with higher CCN efficiency and thereby modulating ACI forcing. Anthropogenic SO₂ emissions (emi_antns_so2) affect ERF_aer uncertainty in most regions, except in the very pristine Antarctic Ocean. Its contribution is generally larger over land (10 %–40 %). Over the oceans, however, other processes dominate: natural DMS emissions (emi_dms) exert a stronger influence over the pristine Southern Ocean (SPO, SAO, and SIO), contributing about 30 %–55 %; while convection and entrainment parameters (conv_entrpen, conv_cprcon, and conv_cmfctop) contribute more substantially over tropical oceans (TPO, TAO, and TIO) where deep convection is active, accounting for 30 %–47 %.

ERF_aci

ERF_aci exhibits a similar uncertainty structure to ERF_aer (Fig. 7b) in most regions, although the relative contributions differ slightly. Compared with ERF_aer, uncertainty in ERF_aci is more strongly influenced by natural emissions (emi_dms; 32 %) and emitted particle size (emi_cmr_bb and emi_cmr_ff; 34 %), while the contribution from anthropogenic SO₂ emissions (emi_antns_so2) is smaller (8 %). Convection parameters (conv_cprcon and conv_entrpen) also contribute more strongly to ERF_aci uncertainty than to ERF_ari. Sea-salt emissions were found to be the dominant natural aerosol source of ERF_aci uncertainty in the HadGEM (Regayre et al., 2018) and UKESM1–A (Regayre et al., 2023) PPEs, whereas they show almost no influence in our PPE and only a similarly weak effect in the ECHAM–HAM PPE (Bhatti et al., 2026). This difference may be partially explained by an underestimation of marine stratocumulus clouds in ICON–HAM and ECHAM–HAM (Neubauer et al., 2019), where sea-salt is expected to exert an influence. These discrepancies highlight the need for MMPPE efforts to better separate parametric and structural uncertainties across models with substantially different structures.

ERF_ari

In contrast, ERF_ari shows a distinctly different uncertainty pattern (Fig. 7c). Besides anthropogenic SO₂ emissions (emi_antns_so2; 18 %), most of the uncertainty stems from parameters controlling aerosol absorption, including anthropogenic BC emissions (emi_ant_bb; 12 %) and refractive indices of BC (rad_bc_ni; 23 %), OC (rad_oc_ni; 26 %), and dust (rad_du_ni; 5 %). The stronger role of absorption compared with other PPE studies likely results from enhanced biomass-burning emissions and lower SO₂ levels in our simulation year 2024–2025 (Jia and Quaas, 2023; Zhong et al., 2024), as well as the inclusion of above-cloud aerosol effects (Ghan, 2013), which are not included in all PPEs, together increasing the relative importance of aerosol absorption.

ERF_ari shows substantial regional variations in parameter contributions (Fig. 7c). The parameter emi_antns_so2 contributes strongly to ERF_ari uncertainty in regions with strongly negative ERF_ari (Figs. 6c and S9), including the industrial regions EUR, ASI, and AFR and their outflow oceans (NPO, TPO, and TIO). Strong sulfate transport from India leads to the largest emi_antns_so2 contribution in TIO, reaching about 60 %. In biomass-burning source regions (NAM, ARC, and SAM) and their downwind regions (EUR, ARCO, and NAO), uncertainty is primarily driven by the refractive index of OC (rad_oc_ni), contributing 23 %–70 %. The refractive index of BC (rad_bc_ni) contributes to ERF_ari uncertainty in nearly all regions, reflecting the efficient long-range transport of small BC to remote areas. Associated with this transport, the coating (so4_coating) and in-cloud wet deposition (wetdep_ic) parameters also contribute to the uncertainty over the Southern Ocean (SAO, SPO, and SIO), by changing the absorption efficiency and by controlling the removal of BC from the atmosphere, accounting for about 10 % and 20 %, respectively. Other parameters related to hygroscopicity (kappa_so4) and chemistry (chem_so4), although identified as major contributors in the ECHAM–HAM PPE, do not appear to be important here.

In addition to the parametric uncertainty, another source of the ERF_aer uncertainty arises from the prescribed PI emission baseline. In existing PPEs, the spatial pattern and reference period of PI emissions are fixed, while the emission scaling parameters are perturbed in the same way as in the PD runs. Different assumptions about PI biomass-burning emissions and the choice of PI period could substantially alter the PI aerosol background and consequently the magnitude and regional pattern of PI-to-PD aerosol forcing (Hamilton et al., 2018; Brown et al., 2026). The PI baseline may further affect the relative contributions of individual parameters presented here. For example, an emission perturbation may exert a stronger effect under a cleaner PI background than under a more polluted background because of the nonlinear response of clouds and radiation to aerosol perturbations (Jia and Quaas, 2023). Although scaling natural aerosol emissions provides some means of constraining the PI aerosol background, uncertainties associated with the prescribed PI emissions are not captured by the current PPE design and cannot be constrained by PD observations alone.

Figure 8 As in Fig. 7, but for (a) AOD, (b) AAOD, and (c) AE. The corresponding spatial maps of uncertainty caused by individual parameters are shown in Figs. S11–S13.

Figure 9 As in Fig. 7, but for (a) CDNC, (b) LWP, and (c) CF_liq. The corresponding spatial maps of uncertainty caused by individual parameters are shown in Figs. S15–S17.

3.3.2 Present-day aerosol properties

Because the PI atmosphere is not directly observable, constraints on PI–PD aerosol forcing must rely on PD observations (Carslaw et al., 2013; McCoy et al., 2020). Identifying the processes controlling key PD aerosol and cloud properties helps elucidate whether and how these quantities can constrain the forcing uncertainty. Here, we focus on three aerosol properties observable from space (Fig. 8), including AOD (total aerosol amount), AE (aerosol size) and AAOD (aerosol absorption), which are closely linked to aerosol extinction and CCN activity. The role of cloud properties (CDNC, LWP, and CF_liq; Fig. 9) is discussed separately in the following Sect. 3.3.3.

AOD

Figure 8a indicates that natural emission parameters dominate global AOD uncertainty, accounting for nearly 50 % of the variance. The main contributors include DMS emissions (emi_dms; 17 %), sea-salt emissions (emi_ss_expo and emi_ss_acc; 21 %), and dust emissions (emi_du_fricv; 8 %). Their influence, however, varies strongly by region (Figs. 8a and S11). DMS and sea-salt parameters primarily affect oceanic regions but also influence nearby continents such as OCE and ANTA, with contributions increasing toward the Southern Ocean. Although dust emissions contribute less than 10 % to global AOD uncertainty, they exert a substantial influence in the regions where its has a significant source (77 % in AFR and 35 % in ASI) and in downwind regions (20 % in TAO and 30 % in EUR). Similarly, despite their small global contributions, biomass burning OC emissions (emi_bb_oc) can explain around 45 % of the uncertainty in NAM and ARC due to boreal fires, and up to 20 % in SAM due to tropical fires, while anthropogenic SO₂ emissions (emi_ant_so2) and sulfate hygroscopicity (kappa_so4) collectively contribute around 30 % of the uncertainties in EUR and ASI due to intense human activities. These strong regional variations suggest that observational constraints should be applied on a regional basis, as different regions are characterised by distinct processes driving the uncertainty (Regayre et al., 2026). In addition, parameters affecting aerosol removal, e.g., conversion of convective cloud water to rain (conv_cprcon) and in-cloud wet deposition (wetdep_ic), contribute relatively uniformly to global AOD uncertainty (40 %), with slightly weaker influence over the Southern Ocean (Figs. 8a and S11).

AAOD

As a component of total AOD, absorbing aerosol optical depth (AAOD) shares three common sources of uncertainty with AOD (emi_du_fricv, wedep_ic, conv_cprcon; Fig. 8a and b). Overall, AAOD uncertainty is dominated by absorption-related parameters, which together explain about 80 % of the variance, in which emission parameters for absorbing aerosols (BC, OC, and dust) account for about 33 %, while their optical property parameters (rad_bc_ni, rad_oc_ni, and rad_du_ni) play an even larger role, contributing ∼ 45 %. These optical property parameters control the imaginary refractive index of BC, OC, and dust, directly governing the absorption efficiency and influencing AAOD, and were also identified as key contributors to single-scattering albedo in the ECHAM–HAM PPE (Bhatti et al., 2026). Although so4_coating has little impact on global AAOD uncertainty, it can affect regions where both coating material (sulfate) and coating cores (biomass-burning or dust particles) are abundant (Fig. S12), by enhancing aerosol absorption through internal mixing of strongly and negligibly absorbing compositions (Bond and Bergstrom, 2006).

AE

Similar to AOD, the Ångström exponent (AE) uncertainty is also dominated by natural emission parameters, explaining an even larger fraction of the variance (about 74 %) than for AOD (Fig. 8c). Parameters associated with coarse-mode emissions (emi_ss_coarse and emi_du_fricv) play a particularly important role. Over land, dust emissions dominate the uncertainty, accounting for about 60 %–70 % in dust-influenced regions (EUR, ASI, and AFR) and ∼ 20 % elsewhere. In oceanic regions, coarse-mode sea-salt emissions dominate the uncertainty (∼ 35 %), with dust additionally contributing in the Northern Hemisphere through long-range transport (Figs. 8b and S13) . Although the wind-exponent parameter in the sea-salt emission scheme (emi_ss_expo) is not directly linked to coarse-mode production, stronger winds preferentially generate coarse sea-salt particles, shifting the emitted size distribution toward larger particles (Gong, 2003; May et al., 2016) and thus influencing AE.

In comparison, fine-mode emission parameters have only a limited influence on AE uncertainty. DMS emissions cause the AE uncertainty through sulfate formation in most regions (∼ 16 % globally), while anthropogenic SO₂ mainly affect their source regions contributing < 5 % globally and < 10 % regionally. Of course, other parameters controlling aerosol sink contribute to AE uncertainty, too, but through different mechanisms. For example, in dust-rich regions, wetdep_bc and so4_coating each contribute about 10 % of the AE variance by modifying precipitation scavenging and dust hygroscopicity through aerosol aging, respectively. Moreover, convection-related parameters (conv_entrpen and conv_cprcon) together cause ∼ 9 % of global uncertainty, by affecting wet removal and aerosol growth through cloud processing (Hoose et al., 2008).

Overall, our uncertainty attribution reveals that present-day aerosol properties share several driving parameters with aerosol forcing, indicating their potential to constrain forcing. For example, emi_dms and conv_cprcon influence AOD, AE, and ERF_aci, while rad_bc_ni, rad_oc_ni, rad_du_ni, and emi_ant_bc affect both AAOD and ERF_ari. The utility of different observations to constrain aerosol forcing is discussed further in Sect. 3.4.

In comparison with the ECHAM–HAM PPE (Bhatti et al., 2026), we identify several common parameters related to aerosol emissions, deposition, and optical properties, but their contributions differ substantially. These differences partly reflect the different expert-elicited parameter distributions used in the two studies, and likely also the inclusion of cloud-related parameters in our PPE, on top of structural differences. By perturbing both aerosol and cloud processes simultaneously, our PPE allows stronger interactions among parameters, which alter the relative contributions of individual parameters. For instance, the convection parameter conv_cprcon alone explains about 20 % of the global AOD uncertainty and causes regional uncertainty of up to ± 0.03 in subtropical regions (Fig. S11), by affecting aerosol removal and regeneration from cloud processing. This highlights the importance of including cloud perturbations even in an aerosol-focused PPE study.

3.3.3 Present-day cloud properties

CDNC

Cloud droplet number concentration (CDNC) is central to aerosol–cloud interactions, linking aerosol perturbations to responses of cloud microphysical (LWP) and macrophysical (CF_liq) properties. In climate models, CDNC is primarily determined by CCN concentrations, updraft velocity, and activation treatment. The activation treatment governs how efficiently aerosol particles activate into cloud droplets for a given updraft. Its uncertainty in our PPE is represented by the parameter activ_aero, whose range reflects the diversity among commonly used activation schemes (Ghosh et al., 2025).

Figure 9a shows that global CDNC uncertainty arises almost equally from parameters controlling CCN concentrations (35 %; emi_dms, emi_ss_expo, drydep_acc, and chem_so2) and activation scheme (33 %; activ_aero), followed by updraft velocity (14 %; activ_w), with activ_aero being the largest single contributor. For CCN-related parameters, DMS emissions dominate globally; however, the relative importance of individual parameters varies strongly by region (Figs. 9a and S15). For example, anthropogenic (emi_antns_so2, emi_ant_oc, and emi_cmr_ff) and biomass-burning (emi_bb_oc and emi_cmr_bb) emission parameters dominate over land and also contribute substantially to ERF_aer uncertainty (Fig. 7), implying strong potential for CDNC observations to constrain ERF_aer, while over oceans the uncertainty is mainly driven by natural emissions (emi_dms, emi_ss_acc, and emi_ss_expo). The parameters drydep_acc and chem_so2 each cause about ∼ 5 % uncertainty in most regions (except for the Southern Ocean) by altering aerosol lifetime through dry deposition and sulfate production through SO₂ oxidation, respectively.

Comparing the three parameter categories (CCN, updraft, and activation), it is interesting to see that CCN-related parameters play a larger role in the cleaner, aerosol-limited Southern Ocean than in the more polluted, updraft-limited Northern Hemisphere. This suggests that ICON–HAM reasonably captures the regime-dependent activation behaviour reported in previous studies using cloud parcel models (Reutter et al., 2009) and satellite observations (Jia and Quaas, 2023). Beyond these three categories, the convection parameter conv_cmfctop contributes about 17 % uncertainty, but only in tropical regions (TPO, TNO, and TIO), by influencing moisture transport and thus cloud formation (Mauritsen et al., 2012).

LWP

The parameters dominating CDNC uncertainty also influence LWP uncertainty, but with much smaller contributions, decreasing from 82 % for CDNC to only 18 % for LWP (Fig. 9b). Instead, LWP uncertainty is dominated by parameters related to cloud microphysics and convection. Globally, auto-conversion rates from cloud water to rain in stratocumulus (micro_ccraut) and convective (conv_cprcon) clouds explain 31 % and 24 % of the LWP uncertainty, respectively. Although convective mass-flux transport (conv_cmfctop) and entrainment mixing (conv_entrpen) account for only 10 % and 6 % of the global uncertainty, they can have a strong impact in tropical regions (AFR, SAM, OCE, TPO, TAO, and TIO), where deep convection is more active. The parameter conv_cmfctop regulates the strength of convective mass-flux transport that redistributes moisture vertically, whereas conv_entrpen controls entrainment mixing between convective updrafts and surrounding air, both affecting the amount of condensed cloud water and thus LWP (also CF_liq; Fig. 9b and c).

CF_liq

The model diagnoses CF_liq solely from relative humidity (Sundqvist et al., 1989); therefore, the parameters dominating CDNC have little influence on CF_liq (Fig. 9c). Uncertainty attribution shows that convection-related parameters (conv_cmfctop, conv_cprcon, and conv_entrscv) account for most of the uncertainty (about 90 %), with the remainder associated with the cloud-cover scheme (cld_cov_crs). Among these, conv_cmfctop is the largest contributor, explaining around 66 % of the CF_liq uncertainty by regulating convective mass-flux upward transport and thereby modifying the moisture available in the boundary layer where liquid clouds form (Mauritsen et al., 2012). conv_entrscv contributes about 8 % by influencing the mass-flux transport more indirectly through entrainment in shallow convection, which modifies the buoyancy of convective updrafts. In addition, conv_cprcon explains about 12 % of the variance by altering convective precipitation efficiency via auto-conversion and thus the persistence of liquid clouds. Only 5 % of the uncertainty arises from the non-convection parameter cld_cov_crs, which determines the critical RH threshold for cloud formation. These parameters produce similar levels of uncertainty across oceanic regions, where liquid clouds are most prevalent (Figs. 9c and S17).

Beyond these parameters, stratiform cloud microphysics parameters controlling phase transitions between cloud water components and precipitation cause the CF_liq uncertainty over land, accounting for about ∼ 10 % from micro_ccraut and about ∼ 30 % from micro_ccsaut. The parameter micro_ccraut affects liquid cloud formation directly by altering rain formation via auto-conversion. In contrast, micro_ccsaut (snow formation via aggregation) affects the CF_liq mainly by changing the partitioning of liquid and ice water in mixed phase clouds through the Wegener-Bergeron-Findeisen process, in which cloud ice grows through the depletion of liquid droplets (Korolev, 2007). This parameter has an even stronger influence on ice and total cloud cover (figure not shown).

To conclude, cloud properties are predominantly controlled by cloud-related parameters (80 % for LWP and 90 % for CF_liq), except for CDNC, for which aerosol and cloud parameters contribute roughly equally. This uncertainty attribution also provides insight into the processes governing aerosol–cloud interactions. The much smaller aerosol contribution to macrophysical properties implies lower LWP and CF susceptibilities to aerosol perturbation compared to CDNC, consistent with the relatively weak adjustments in LWP and CF_liq relative to the Twomey effect reported in climate models (Gryspeerdt et al., 2020; Duran et al., 2025). This, however, contrasts with observational evidence that suggests that cloud adjustments (particularly CF adjustment) can be comparable to or even stronger than the Twomey effect (Rosenfeld et al., 2019; Chen et al., 2022; Yuan et al., 2023; Wall et al., 2023). Such discrepancies may point to structural deficiencies in the model, i.e., missing or poorly represented processes, such as aerosol effects in cloud fraction parameterisation (Sundqvist et al., 1989) and turbulence–microphysics interactions in cloud-top entrainment (Mülmenstädt et al., 2024).

3.4 Sensitivity of model outputs to model parameters

Sensitivity of model outputs to perturbed parameters captures not only the magnitude but also the direction of parameter effects, complementing the variance-based analysis presented above. Following Eidhammer et al. (2024), it is defined as the linear regression slope of standardized outputs against normalized parameter values, allowing direct comparison across parameters and variables. More importantly, the sensitivity analysis helps identify which PD observations with the greatest potential to constrain PI-to-PD aerosol forcing. Specifically, if a PD observable is strongly sensitive to the same parameters that drive uncertainty in PI-to-PD forcing, it has a strong potential to serve as an effective observational constraint.

A wide range of PD observations from both remote sensing and in situ measurements have been used to constrain ERF_aer in PPEs, including AOD (Johnson et al., 2020; Watson-Parris et al., 2020; Prévost et al., 2026), aerosol index (AI) (Regayre et al., 2023), aerosol and precursor gas concentrations (Johnson et al., 2020; Prévost et al., 2026), CCN (Song et al., 2025), CDNC (Regayre et al., 2023, 2026; Mikkelsen et al., 2025; Song et al., 2025), COT (Regayre et al., 2023, 2026), CER (Regayre et al., 2023, 2026), LWP (Regayre et al., 2023, 2026; Eidhammer et al., 2024; Mikkelsen et al., 2025), CF_liq (Regayre et al., 2023, 2026), precipitation rate (Mikkelsen et al., 2025), radiative fluxes (Regayre et al., 2023, 2026; Eidhammer et al., 2024), and cloud radiative effects (Eidhammer et al., 2024), among others. However, how effectively these observations constrain ERF_aer remains unclear and may vary across models due to structural differences.

Figure 10 Linear regression slopes of standardized global-mean outputs against normalized parameter values for 20 variables. The output variables are divided into PI-to-PD forcing (left) and present-day (PD) observables (right), separated by a dashed line. The parameter order follows Tables 1 and 2. Parameters identified as key drivers of ERF_aer uncertainty are highlighted in bold blue. A positive correlation means that the model diagnostic increases with the parameter value, whereas a negative correlation indicates the opposite. The slope values are also shown in each grid cell for quantitative comparison. Similar results for land and ocean are presented in Figs. S19 and S20, respectively.

Here, we systematically evaluate the constraint potential of commonly observed aerosol and cloud properties from both space- and ground-based measurements. Figure 10 shows the sensitivities of global-mean outputs to parameter changes, with PI-to-PD forcing on the left and PD observables on the right. Given the level of detail in the figure, we focus on the global-mean results in the main text, while corresponding land and ocean results are provided in Figs. S19 and S20. Red colours represent increases in the output with increasing parameter values, whereas blue colours represent decreases. In the following, we focus on parameters whose slopes exceed 1, indicating strong influence on the outputs (a slope of 1 corresponds to one standard deviation across the full parameter range), and discuss them in top-to-bottom order in Fig. 10.

As shown in Fig. 10, the sensitivity analysis identifies the same driving parameters as the uncertainty attribution in Sect. 3.3, i.e., emi_antns_so2, emi_dms, emi_cmr_ff, emi_cmr_bb, conv_cprcon, and conv_cmfctop for ACI and emi_antns_so2, emi_ant_bc, rad_bc_ni, rad_oc_ni, rad_du_ni, so4_coating and conv_cprcon for ARI. Encouragingly, for each of these parameters driving aerosol forcing, at least one PD observable shows strong sensitivity, indicating potential for constraining ERF_aer using existing measurements.

Anthropogenic SO₂ emissions (emi_antns_so2) exhibit a strong negative slope with respect to both ERF_aci and ERF_ari, consistent with enhanced PD–PI aerosol contrast under higher anthropogenic emissions. Three PD diagnostics related to fine-mode aerosol concentrations also show strong sensitivities to emi_antns_so2, including AI (AOD multiplied by AE) and CCN number concentrations at the surface (CCN_S0.3) and at the top of boundary layer (CCN_BL0.3). Compared with AOD or AE alone, AI better represents fine-mode aerosol abundance by combining aerosol loading and size information (Nakajima et al., 2001), and is thus more closely related to SO₂ emissions.

Anthropogenic BC emissions (emi_ant_bc) show a positive correlation with ERF_ari, reflecting increased atmospheric absorption by human-emitted BC since the industrial era. This parameter is also positively correlated with PD AAOD and negatively with single scattering albedo (SSA) because of its strong link to atmospheric absorption.

DMS emissions (emi_dms) influence aerosol forcing primarily through ACI by acting as a source of CCN. The parameter emi_dms shows a strong positive correlation with ERF_aci, suggesting that increasing DMS emissions enhances the natural background CCN concentration, thereby reducing the PD–PI contrast and leading to a more positive forcing. Notably, more than half of the PD observables closely linked to CCN concentrations are highly sensitive to emi_dms. These are primarily commonly used CCN proxies, with stronger sensitivities generally observed for better proxies, i.e., CCN_BL0.3 > AI > fine-mode AOD (AOD_f) > AOD (Gryspeerdt et al., 2017; Jia et al., 2026), implying that sulfate from DMS oxidation is an important source of CCN, particularly over the oceans (Fig. S20). Although not directly related to CCN, PD AE, CDNC, LWP, and even TOA net radiative flux (F_net) and outgoing shortwave radiative flux ($F_{\mathrm{sw}}^{\uparrow }$) also show sensitivity to emi_dms, reflecting the broad influence of DMS-derived marine secondary aerosols on global clouds and radiation by acting as CCN.

Emitted particle sizes from fossil fuel (emi_cmr_ff) and biomass burning (emi_cmr_bb) are also correlated with ERF_aciby modulating particle number for a given emitted mass. However, only surface-level PD observables (CCN_S0.3 and BC_S) show notable sensitivity (with a slope close to 1) to these parameters, likely due to their closer association with emission sources. These observations are spatially sparse compared to satellite-retrieved quantities. Nevertheless, ground-based CCN and BC measurements, particularly at sites close to emission source regions, would provide valuable constraints.

Refractive indices related to anthropogenic aerosols (rad_bc_ni and rad_oc_ni) exhibit strong positive correlations with ERF_ari, indicating that the more absorbing the anthropogenic aerosols, the more positive the ARI forcing. Refractive index of dust (rad_du_ni), in contrast, is negatively correlated to ERF_ari, likely because dust provides a strongly absorbing background that dampens the anthropogenic radiative perturbation of other aerosols. Importantly, several PD diagnostics are also strongly linked to these parameters, including AAOD, SSA, surface BC mass concentration (BC_S), PM_2.5 and AE.

The conversion rate from cloud water to rain in convective clouds (conv_cprcon) correlates positively with ERF_aci but negatively with ERF_ari, owing to its different effects on cloud and aerosol processes. A higher conv_cprcon reduces the LWP and CF_liq baseline, weakening the PD-PI contrast in cloud radiative effect (CRE) and thereby shifting ERF_aci towards a more positive value. In contrast, the reduced cloud amount diminishes above-cloud aerosol absorption (Ghan, 2013), resulting in a more negative ERF_ari. The parameter conv_cprcon is very promising to be constrained, as many PD observables show strong sensitivity to it, including F_net, $F_{\mathrm{sw}}^{\uparrow }$, LWP, CF_liq, AOD, AI, AE, AAOD, AOD_f, and coarse-mode insoluble and soluble AOD (AOD_ci and AOD_cs).

Similar to conv_cprcon, conv_cmfctop is also positively correlated with ERF_aci through its influence on cloud properties. A larger conv_cmfctop enhances the convective moisture transport from the boundary layer into the free troposphere, thereby suppressing cloud formation and reducing the baseline LWP, CF_liq, and CRE, which in turn leads to a more positive ERF_aci. This parameter is also readily constrained by observations, as most PD cloud (LWP and CF_liq) and radiation (F_net and $F_{\mathrm{sw}}^{\uparrow }$) diagnostics show strong sensitivity to it.

To summarise, for each of the key parameters driving ERF_aer uncertainty, at least one–and often multiple–PD observables are sensitive to it. This underscores the considerable potential to constrain ERF_aer using existing in situ and remote sensing measurements. For parameters that can be readily constrained by multiple PD observations, a key challenge thus becomes to prioritise observed variables during the constraint process, favouring those with lower measurement errors and, potentially, those less affected by model internal inconsistencies that may lead to failed constraints (Regayre et al., 2023; Prévost et al., 2026). This could be achieved using well-defined constraint measures, e.g., the implausibility metric (Johnson et al., 2020), which weights constraint variables according to their associated uncertainties across different aspects. In addition, several parameters (particularly cloud parameters related to activation, auto-conversion, aggregation, and entrainment) while not directly driving aerosol forcing, exert strong control on cloud states (CDNC, LWP, and CF_liq) and thus would be crucial for cloud feedbacks, too (Duffy et al., 2024; Gettelman et al., 2024). Although our PPE is not designed to explicitly explore climate feedbacks, constraining simulations that resemble PD clouds is expected to improve the representation of cloud feedbacks as a by-product, thereby equipping the model with improved capability for climate projections.

3.5 Comparison to observations

This section evaluates aerosol and cloud properties from the ICON–HAM default run against SPEXone and MODIS-COSP observations. The 3 h model outputs are co-located with the observations in time and space to mitigate spatiotemporal representation uncertainties (Schutgens et al., 2016a, b). Comparison of spatial patterns between model bias and parametric uncertainty is also presented, which would inform whether the model bias is correlated with parametric uncertainty. To aid process-level attribution of the bias, global maps of uncertainty contributions from individual parameters are provided in Figs. S11–S13 and S15–S17.

Figure 11 Global annual-mean AOD (top), AAOD (middle), and AE (bottom) from SPEXone retrievals (left), the ICON–HAM CTRL simulation (middle), and their difference (ICON–HAM minus SPEXone; right). The 3 h ICON–HAM outputs are co-located with SPEXone observations in space and time. Green isolines enclose regions with parametric uncertainty in the highest 10 % globally. Top-left numbers indicate global mean values.

3.5.1 Aerosol properties

Figure 11 shows annual mean AOD, AAOD, and AE from the ICON–HAM default run against SPEXone retrievals. Overall, the model appears to underestimate AOD compared to observations (Fig. 11a–c), with a global mean bias of −0.036. The underestimation of AOD is also common in other AeroCom models (Kinne et al., 2006; Gliß et al., 2021). Spatially, the largest differences are primarily found over land, particularly in biomass burning regions over Southern America and Africa tropical rain forests, where the model exhibits pronounced negative biases (Fig. 11c). A similar, though weaker, negative bias occurs in northern boreal regions. Similar bias patterns have also been reported in ECHAM-HAM and an earlier version of ICON–HAM, and have been attributed to too low biomass burning emissions and too strong wet deposition by convective clouds (Tegen et al., 2019; Salzmann et al., 2022). With the PPE, we suggest the dominant role of biomass-burning emissions, as the spatial pattern of the negative bias closely aligns with the parametric uncertainty caused by emi_bb_oc (Fig. S11). In-cloud wet deposition (wetdep_ic) and convective cloud-to-rain conversion rate (conv_cprcon) also contribute to the bias in these regions, albeit to a lesser extent (Fig. S11). In addition, negative biases of AOD are evident in major industrial regions (e.g., China and India) as well. These biases are likely linked to the high uncertainty associated with emissions of anthropogenic aerosols (emi_antns_so2) and their hygroscopicity (kappa_so4) in these regions (Fig. S11). The positive biases, in contrast, occur primarily in the Sahara and its outflow regions, whose spatial pattern is fairly compatible with that of the uncertainty related to dust emissions (emi_du_fricv; Fig. S11). Consistent with its sister model ECHAM–HAM (Bhatti et al., 2026), ICON–HAM overestimates AOD over subtropical oceans, but less strongly (Salzmann et al., 2022). This bias may be explained by relatively large parametric uncertainties caused by both sources (emi_dms, emi_ss_acc, emi_ss_coarse, emi_ss_expo) and sinks (wetdep_ic and conv_cprcon) of natural aerosols (Fig. S11). Overall, regions with larger biases often coincide with high parametric uncertainty (green isolines; Fig. 11c), implying that part of these biases may be reduced through parameter optimisation.

AAOD largely mirrors the AOD bias map, with a global mean bias of −0.003 (Fig. 11d–f). There are, however, two notable exceptions. First, the positive bias over subtropical oceans seen in AOD disappears in AAOD, as these regions are dominated by scattering aerosols (DMS-derived sulfate and sea-salt). Second, the negative AOD bias over China and India shifts to a positive bias in AAOD, likely reflecting large uncertainties in atmospheric absorption, e.g., anthropogenic emissions of BC (emi_ant_bc) and its refractive index (rad_bc_ni) (Fig. S12). Figure 11f illustrates that virtually all regions with large biases closely correspond to areas of high parametric uncertainty, suggesting that parameter optimisation offers also an effective pathway to reduce the AAOD bias.

AE exhibits an overall positive bias relative to SPEXone observations (Fig. 11g–i), with a mean bias of 0.19, pointing to too small particle sizes in ICON–HAM. Over land, the largest positive biases occur in dust source regions in the western US, where the model fails to reproduce the lower AE compared to the eastern US observed in the data. A similarly strong positive bias is found over the Tibetan Plateau, likely attributed to limitations in simulating aerosol transport over complex terrain at high altitudes, as well as reduced retrieval accuracy under complex terrain and low aerosol. Such issues are also seen in AeroCom models (Kinne et al., 2006; Gliß et al., 2021). Another continental region showing positive bias is South America. As shown in Fig. S13, the origin of the AE bias over land is complex due to the interplay between fine- and coarse-mode aerosols across multiple processes, including emissions (emi_dms and emi_du_fricv), particle aging (so4_coating), and wet deposition (wetdep_ic and wetdep_bc). Although rad_oc_ni also contributes to AE uncertainty over land, this does not reflect changes in particle size. Instead, changes in the OC absorption can alter aerosol extinction differently across wavelengths and thereby affect AE (Schuster et al., 2006).

Figure 12 As in Fig. 11, but for comparisons of CDNC (top), LWP (middle), and CF_liq (bottom) between the ICON–HAM MODIS COSP simulator and the MODIS COSP product.

Strong positive AE biases are broadly found over oceans, particularly over large-scale convergence zones, including the Intertropical Convergence Zone (ITCZ), the South Pacific Convergence Zone (SPCZ), the Indian Ocean Convergence Zone (IOCZ), and the South Atlantic Convergence Zone (SACZ) (Fig. 11i). These regions are characterised by intense moisture transport, active convective cloud formation, and strong precipitation, where precipitation is often biased in models (Giorgetta et al., 2018; Neubauer et al., 2019). Notably, the spatial patten of the positive bias agrees well with that of the parametric uncertainties associated with emissions of DMS (emi_dms) and sea-salt (emi_ss_coarse and emi_ss_expo), as well as convective processes such as cloud-to-rain conversion (conv_cprcon) and entrainment in deep convection (conv_entrpen) (Fig. S13). This coincidence, together with the AOD underestimation over oceans (Fig. 11c), suggests that insufficient coarse-mode sea-salt emissions and overly strong wet scavenging in convective clouds contribute to the AE overestimation. Comparing model bias with parametric uncertainty (Fig. 11i), the green isolines largely enclose regions of high parametric uncertainty, highlighting the potential for reducing AE bias through parameter optimisation in most areas. Nonetheless, part of the convergence zones has a positive bias that cannot be explained by the parametric uncertainty alone, pointing to possible structural deficiencies in the convection parameterisation.

To provide a more direct overview of the parameters driving uncertainty in regions with both large biases and high parametric uncertainty, we repeat the analyses of Fig. 8 considering only locations within the green isolines (Fig. S14). It is shown that the dominant parameters are overall consistent with those inferred from the spatial patterns discussed above, although regional averaging can obscure subregional signals, such as the role of emi_bb_oc over central Africa. These results can inform which parameters may be most useful to constrain observationally for reducing model biases in each region.

3.5.2 Cloud properties

Before comparing the model with MODIS, it is worth noting that MODIS retrievals of cloud microphysics (CDNC and LWP) are less reliable over polar regions and land than over mid-to-low-latitude oceanic regions (Grosvenor et al., 2018).

Figure 12 compares model diagnostics and MODIS observations processed consistently using the COSP simulator. ICON–HAM exhibits positive CDNC biases in most regions, except in polar and some land areas where MODIS cloud retrievals are subject to larger errors (Grosvenor et al., 2018). Over oceans, the largest positive biases occur in high-latitude regions, followed by subtropical eastern ocean stratocumulus decks (Fig. 12a–c). This bias pattern is very similar to the spatial distribution of parametric uncertainties associated with natural emissions of DMS (emi_dms) and sea-salt (emi_ss_acc and emi_ss_expo), as well as activation (activ_aero) and vertical velocity (activ_cwtrub), as shown in Fig. S15. Over land, a pronounced CDNC overestimation is found in China and India, despite an underestimation of AOD (Fig. 11c). This positive bias may be thus attributed to too many fine particles from anthropogenic emissions (emi_cmr_ff) and/or too efficient activation (activ_aero) in the model (Fig. S15). Moreover, the good spatial consistency between mode bias and parametric uncertainty (Fig. 12c), suggests that most of the CDNC bias would be reduced by constraining these uncertain parameters.

For LWP, ICON–HAM shows a large negative bias over high-latitude continents (Fig. 12d–f), yet no parameters can explain this bias patten (Fig. S16). This is most likely attributable to large retrieval errors caused by high solar zenith angles and the presence of mixed-phase clouds (Khanal and Wang, 2018). In contrast, in the remaining regions where LWP retrievals are considered more reliable, the model generally overestimates LWP, particularly in the ITCZ. As shown in Fig. S16, processes of auto-conversion in both stratocumulus and convective clouds and entrainment in deep convection, might be responsible for this bias. Again, most regions with large biases coincide with areas of high parametric uncertainty, whereas the remaining regions are likely influenced by retrieval biases and/or model structural uncertainties. For the locations within the green isolines, the parameter contributions to LWP and CDNC to uncertainties in each region (Fig. S18) are broadly consistent with the spatial correspondence discussed above.

CF_liq is underestimated globally by around 25 % (Fig. 12g–i), which is particularly pronounced in regions with persistent low-level cloud decks over subtropical and polar oceans. This may partly explain the weaker ERF_aci in offshore stratocumulus regions in ICON–HAM than in other AeroCom and CMIP6 models (Gryspeerdt et al., 2020; Kalisoras et al., 2024). Unlike other properties, the CF_liq bias cannot be associated to parametric uncertainty, as evidenced by the little spatial overlap between regions of large bias and those of high uncertainty (Fig. 12i). Instead, structural improvements in cloud parameterisation are required, for example, moving from a RH-based scheme (Sundqvist et al., 1989) to a prognostic cloud cover scheme (Muench and Lohmann, 2020). The underestimation of cloud cover in stratocumulus region has been identified as an inherent issue in ICON 2.6.4 (Grundner et al., 2025), and is expected to be addressed in upcoming versions.

4 Conclusions and outlook

In this study, we constructed a perturbed parameter ensemble (PPE) using the ICON–HAM model, with a focus on aerosol effective radiative forcing. A total of 42 aerosol and cloud parameters were perturbed across 13 parameter categories, including aerosol emissions, optical and hygroscopic properties, deposition, chemistry, activation, cloud microphysics, convection, and turbulence, resulting in an ensemble of 383 members. Gaussian Process emulators were then trained on these 383 members to continuously sample the full parameter space. This work aims to quantify and attribute parametric uncertainties in aerosol effective radiative forcing, to identify observable state variables that have the potential to constrain these uncertainties, and to assess the extent to which model biases can be explained by parametric uncertainty. These insights provide a foundation for subsequent observational constraints.

Prior to the main analyses, we first evaluated the performance of the Gaussian process emulators. The coefficients of determination between predicted and reference global means exceed 0.8 for AOD, AE, CDNC, LWP, and CF, and are only slightly lower for ERF_aer (0.76). We note that the emulator shows a positive bias for strongly negative ERF_aer, which arises primarily from ACI. Spatial patterns are also very well captured, with relative biases of all properties below 4 % across most regions globally. Overall, our emulator provides a robust representation of the reference simulations in terms of both global mean and spatial distribution. Furthermore, we illustrated that the emulator performance improves with increasing simulations-to-parameters ratio (r_sp) but levels off around 8–9. A r_sp of 5 is identified as near-optimal because it provides a great balance between accuracy and computational cost, with only marginal improvements at higher values. This result thus offers important guidance for future PPE design, particularly for MMPPE frameworks that involve substantial computational and data storage demands, and suggests that further improvements in emulation accuracy are more likely to come from advances in emulator design rather than increasing ensemble size. For PPEs sampling comparable parameter spaces, and thus subject to similar levels of non-linearity arising from parameter interactions, we would not expect the near-optimal r_sp to vary substantially across models. Nevertheless, this assumption will be interesting to test with the upcoming MMPPE.

Parametric uncertainty in aerosol and cloud processes within our emulated PPE leads to 5–95th percentile ranges of [−1.42 to −0.65] W m⁻² for ERF_aer, [−1.37 to −0.67] W m⁻² for ERF_aci, and [−0.13 to 0.23] W m⁻² for ERF_ari. ERF_aci is uncertain only in magnitude, whereas ERF_ari is uncertain in both magnitude and sign. Aerosol parameters, particularly those related to emissions, dominate the overall ERF_aer uncertainty (80 %). Anthropogenic SO₂, natural DMS, and emitted particle size from biomass burning and fossil fuel contribute equally of around 25 % for each. This underscores the central role of PI-to-PD aerosol perturbations in shaping aerosol forcing. Although cloud parameters play a minor role (13 %) globally, convection and entrainment parameters together account for up to ∼ 50 % of the uncertainty in tropical oceans, where deep convection is prominent. As ERF_aer is largely driven by ACI, the parameters governing ERF_aci uncertainty are similar, but with relatively stronger contributions from DMS emissions and emitted particle size than from anthropogenic SO₂ emissions. In contrast, ERF_ari is largely controlled by absorption-related parameters (60 %), including anthropogenic BC emissions and the refractive indices of BC, OC, and dust, while aerosol scattering associated with anthropogenic SO₂ emissions contributes only less than 20 %. Although shipping emissions are an important target for reductions, e.g., through the IMO 2020 regulations (Diamond, 2023), they contribute only a small ERF_aer uncertainty (0.01 W m⁻²) here, mainly due to their smaller magnitude compared with other anthropogenic emission sectors and the underestimated marine cloud cover limiting their effect through ACI in our model.

Our PPE suggests that aerosol parameters contribute much more strongly than cloud parameters to the ERF_aer uncertainty. Nevertheless, this conclusion is expected to depend on the choice of perturbed parameters. Regayre et al. (2018) and Yoshioka et al. (2019) likewise found a stronger influence of aerosol parameters on ERF_aer than cloud parameters in HadGEM. In contrast, studies using CAM6 (Gettelman et al., 2024) and E3SM (Nugent et al., 2026) reported that cloud parameters contribute more strongly than aerosol parameters. This is partly because our PPE scales the autoconversion rate itself, whereas the CAM6 and E3SM PPEs scale the exponent. In addition, we included a relatively large number of aerosol parameters associated with anthropogenic and natural emissions, which exert a first-order control on PD–PI aerosol changes. These contrasting results again highlight the need for MMPPE experiments, in which different models perturb a consistent set of parameters.

When comparing ERF_aer uncertainties across different PPE studies (Fig. 5), it is also important to recognise that these uncertainties are often defined in different contexts and are therefore not directly comparable. These contexts include (i) unconstrained ranges derived from uniform parameter distributions within expert-elicited bounds, (ii) expert-constrained ranges based on prior knowledge of parameter probability distributions, and (iii) observation-constrained ranges obtained by sub-selecting plausible members from an initial unconstrained ensemble in (i) using observational data. In this study, we compare our PPE only with estimates not yet constrained by observations.

As for present-day properties, the global AOD uncertainty is driven by natural emissions (DMS, sea-salt, and dust) (50 %) and aerosol removal processes (conversion of convective cloud water to rain and in-cloud wet deposition) (40 %). The effect of anthropogenic emissions is only visible in major industrial regions (EUR and ASI), with little global influence. Natural aerosol emissions account for more than 70 % of the AE uncertainty, with coarse-mode sea-salt being the largest contributor (about 30 %), while absorption-related parameters collectively explain about 80 % of AAOD uncertainty. The dominance of natural aerosols in AOD and AE uncertainty is partly due to the wide spatial extent of sea-salt and DMS emissions. Notably, there are common drivers of PD aerosol properties and aerosol forcing, indicating some potential for constraining ERF_aer. For example, AOD and AE share the parameters related to DMS emissions (emi_dms) and convective cloud water-to-rain conversion (conv_cprcon) with ERF_aci, while AAOD shares the parameters related to refractive indices of BC (rad_bc_ni), OC (rad_oc_ni), and dust (rad_du_ni), and anthropogenic BC emissions (emi_ant_bc) with ERF_ari. Through modulating aerosol removal and regeneration from cloud processing, convection parameters can exert a considerable effect on aerosol properties, highlighting the importance of simultaneously perturbing aerosol and cloud parameters to account for their strong interactions.

Uncertainty in CDNC arises from both aerosol and cloud parameters, including those controlling CCN concentrations (natural aerosol emissions, deposition, and chemistry; 35 %), the activation scheme (33 %), and updraft velocity (14 %). In contrast, LWP and CF_liq are dominated by cloud parameters controlling microphysics, convection, and entrainment, with aerosol parameters accounting for merely 5 % of LWP uncertainty and virtually zero for CF_liq. These three cloud properties collectively share three driving parameters with aerosol forcing (emi_dms, conv_cprcon, and conv_cmfctop), underscoring their value for constraining the ERF_aer uncertainty.

To further clarify which PD observations can best constrain PI-to-PD forcing, we systematically assessed the constraint potential of commonly observed aerosol and cloud properties from both satellite and ground-based measurements. The sensitivity of model outputs to parameters revealed that many PD observables share dominant causes of uncertainty with ACI and ARI forcing, mainly through the processes of natural and anthropogenic emissions, aerosol absorption, and cloud convection, suggesting a strong potential for constraint by existing measurements. Specifically, the uncertain parameters dominating ERF_ari could be constrained using the observations of aerosol absorption (e.g., AAOD, SSA) and fine-mode aerosol amount (e.g., AI, CCN_BL0.3 and CCN_S0.3). For ERF_aci, most controlling parameters can likewise be constrained by a range of aerosol and cloud observations that are readily obtained from space with global coverage (e.g., F_net, $F_{\mathrm{sw}}^{\uparrow }$, LWP, CF_liq, AOD, AI, AE, AAOD, AOD_f, AOD_ci, and AOD_cs). An exception is emitted particle size that particularly requires surface-level measurements of CCN and BC concentrations. With multiple PD observables showing constraint potential, the key challenge lies in prioritising those with lower measurement errors, denser spatial and temporal sampling, and, potentially, those less affected by model internal inconsistencies that may result in failed constraints (Regayre et al., 2023). On the other hand, knowledge on conflicting constraints from diverse PD observables characterising different physical processes, are also valuable in a way that helps expose potential model deficiencies (Prévost et al., 2026).

Comparison against SPEXone and MODIS–COSP observations shows a strong spatial correspondence between model biases and parametric uncertainties in key aerosol and cloud quantities. This suggests that much of the detected biases in our model could be mitigated through appropriate parameter optimisation (observational constraints). However, notable exceptions remain in regions where model biases cannot be attributed to parametric uncertainty. For example, the strong underestimation of LWP over high-latitude continents may partly reflect less reliable satellite retrievals in these regions, while the widespread negative CF_liq biases in stratocumulus regions are more likely linked to deficiencies in cloud cover parameterisation in the current ICON–HAM version. The current RH-based cloud cover scheme is a simplification (Sundqvist et al., 1989) and a prognostic cloud cover scheme may help reduce these biases (Muench and Lohmann, 2020).

Overall, the range of behaviour captured by the ICON–HAM PPE underscores the value of PPE framework for investigating complex interactions among aerosol, cloud and radiation by tracking back their driving processes. Although observational constraints are beyond the scope of this study, the systematic identification of common causes in uncertainties of PD observables and PI-to-PD state changes provides a pathway toward more causally robust constraints on ERF_aer. Our PPE also offers a foundation for future ICON–HAM development, particularly in refining uncertain aerosol and cloud representations, and the outcomes are likely applicable to higher-resolution model configurations with similar microphysics schemes (Weiss et al., 2025).

Improving model skill in simulating aerosol forcing and predicting future climate requires a better understanding of both parametric and structural uncertainties, which calls for coordinated progress across observational and modeling communities:

• Additional measurement types are needed, especially those directly linked to emitted particle size (Regayre et al., 2018; Zhong et al., 2024; Bhatti et al., 2026), which remains weakly constrained by existing observations. Moreover, robust estimates of observational uncertainty are essential for appropriately weighting diverse observations in the constraint process, and for disentangling structural model deficiencies from observational errors (Regayre et al., 2023; Prévost et al., 2026). Such pixel-level uncertainty information is, however, available for only a few observational products, for example SPEXone (Fu et al., 2025), highlighting a clear need for systematic and consistent uncertainty reporting by the observational community.

• Accordingly, more elaborate constraint frameworks are required to explicitly account for these observational uncertainties that differ in different products, which would help reduce inconsistencies as additional measurements are included (Johnson et al., 2020; Regayre et al., 2026) and yield more robust constraints on ERF_aer through maximal use of existing observations.

• ERF_aer uncertainties from different PPEs are not readily comparable due to differences in experimental setup (Fig. 5), making it difficult to separate parametric and structural model uncertainties. To tackle this issue, a coordinated multi-model PPE framework is needed (Carslaw et al., 2026) and is currently being developed within CMIP7 and AeroCom Phase 4 under the ACI-PPEMIP initiative (https://wcrp-cmip.org/mips/, last access: July 2026), largely building on the ICON–HAM protocol presented here. Such a framework would enable significant advances in climate projection skill beyond what can be achieved with single-model PPEs.

Appendix A: Region definition

Figure A1 Definition of the 18 regions adopted from Jia et al. (2021).